Introduction {#s1}
============

Lung cancer is a heterogeneous malignancy with poor survival due to diagnosis at an often advanced stage \[[@B1]\]. Lung cancer is broadly divided into small cell lung cancer (\~15% of all lung cancers) and non-small cell lung cancer with adenocarcinoma as the most frequent histological type \[[@B2]\]. In adenocarcinoma, the mutational spectrum is dominated by *EGFR* and *KRAS* mutations, where the former is an established predictor of response to EGFR inhibitors \[[@B3],[@B4]\]. *EGFR* and *KRAS* mutations are nearly always mutually exclusive and associated with differences in patient gender and smoking history \[[@B5]\]. Together, this suggests that these genetic alterations may be drivers of pathogenesis for specific adenocarcinoma subgroups \[[@B5]\] (and references therein). In *EGFR* and *KRAS*-wild type adenocarcinomas (EGFRwt/KRASwt), different potential drivers of pathogenesis exist, including *ALK*, *RET*, and *ROS1* gene fusions, with *ALK* rearrangements being therapeutically relevant \[[@B5]-[@B7]\]. Several studies have reported genomic or transcriptional alterations between *EGFR*-mutated and/or *KRAS*-mutated tumors and corresponding wild-type adenocarcinomas \[[@B8]-[@B18]\]. However, the majority of previous studies are based on relatively small patient cohorts and do not always stratify tumors into all three mutational groups, which may explain conflicting results. *EGFR*-mutated adenocarcinomas have repeatedly been associated with the bronchioid gene expression subtype originally defined by Hayes et al. \[[@B19]\]. Bronchioid classified tumors are generally of lower grade, have a higher expression of excretion, asthma and surfactant genes, occur predominantly in women and never-smokers, and have better overall survival compared with the other two expression subtypes, magnoid and squamoid \[[@B19],[@B20]\]. The magnoid and squamoid subtypes harbor more *KRAS* mutations, seem to be more closely related in gene expression, occur more often in men and smokers, and have poorer overall survival \[[@B19],[@B20]\].

To resolve conflicting reports and provide a comprehensive survey of copy number alterations, allelic imbalances and transcriptional alterations in lung adenocarcinomas stratified by *EGFR* and *KRAS* mutation status, we analyzed 457 genomic and 914 gene expression profiles for differences between the three mutational groups ([Figure 1](#pone-0078614-g001){ref-type="fig"}). We show that a few consistent genomic differences exist between the mutational groups, however with moderate frequencies. Transcriptional analyses identified only a small set of differentially expressed genes across multiple cohorts, and highlighted the difficulty to resolve the three mutational groups as distinct transcriptional entities. Together, our results suggest that the genomic and transcriptional landscape of lung adenocarcinoma is only to a minor extent determined by the mutational status of *EGFR* and *KRAS*.

![Schematic diagram of genomic and transcriptional analyses.\
(**A**) Genomic analyses. (**B**) Transcriptional analyses. Individual cohorts are portrayed.](pone.0078614.g001){#pone-0078614-g001}

Materials and Methods {#s2}
=====================

Tumor material {#s2.1}
--------------

Genomic profiles from 1272 adenocarcinoma tumors and cell lines were obtained from a previous study (n=1210) \[[@B21]\], with addition of adenocarcinomas from Wilkerson et al. \[[@B20]\] (n=62, GSE36363). All genomic profiles were analyzed in an unmatched fashion and sample uniqueness was assured as described \[[@B21]\]. All included tumors represented primary disease based on information from original studies. *EGFR* and *KRAS* mutational status was available for 457 adenocarcinoma, including 113 *EGFR*-mutated, 134 *KRAS*-mutated, and 210 EGFRwt/KRASwt tumors.

Gene expression profiles from 914 adenocarcinoma tumors, including 309 *EGFR*-mutated, 192 *KRAS*-mutated and 413 EGFRwt/KRASwt cases, were collected from eight studies analyzed by different microarray platforms \[[@B8],[@B10],[@B20],[@B22]-[@B26]\]. Samples from Chitale et al. \[[@B10]\] were further divided into two cohorts according to their different Affymetrix platforms (U133A and U133 2plus).

Explicit information on patient ethnicity or specific mutation type was not available for the majority of the included studies, and these parameters were therefore omitted from the analyses. However, the included studies were performed in both western and Asian countries. Patient and tumor characteristics are summarized in [Tables 1](#pone-0078614-t001){ref-type="table"} and [2](#pone-0078614-t002){ref-type="table"}, with additional description in [File S1](#pone.0078614.s005){ref-type="supplementary-material"}.

10.1371/journal.pone.0078614.t001

###### Characteristics of individual aCGH and SNP genomic adenocarcinoma cohorts stratified by microarray platform.

  Genomic cohort           No. of AC^A^   Microarray platform         Included in CN / GAP analysis^B^   No. of tumors / cell lines   No. of *EGFR*-mutated / *KRAS*-mutated / EGFRwt/KRASwt^C^   No. of stage I/II/III/IV tumors^C^   Gender female / male^C^   Smoking status NS/S^C,D^
  ------------------------ -------------- --------------------------- ---------------------------------- ---------------------------- ----------------------------------------------------------- ------------------------------------ ------------------------- --------------------------
  Zhao \[[@B49]\]          36             Affymetrix 100K             Yes/No                             36/0                         \-                                                          \-                                   2/1                       \-
  GSE18252 \[[@B50]\]      4              Affymetrix 100K             Yes/Yes                            4/0                          0/4/0                                                       \-                                   \-                        \-
  Weir \[[@B51]\]          112            Affymetrix 250K Sty         Yes/No                             112/0                        \-                                                          \-                                   \-                        \-
  Weir \[[@B51]\]          196            Affymetrix 250K Sty         Yes/No                             196/0                        15/49/51                                                    41/14/15/3                           110/83                    17/116
  GSE19399 \[[@B52]\]      19             Affymetrix 250K Sty         Yes/Yes                            19/0                         1/1/8                                                       10/1/4/0                             10/9                      2/16
  GSE17247 \[[@B53]\]      49             Affymetrix 250K Sty         Yes/No                             0/49                         \-                                                          \-                                   \-                        \-
  GSE28572 \[[@B27]\]      44             Affymetrix 250K Nsp         Yes/Yes                            44/0                         9/18/17                                                     21/10/8/2                            27/17                     4/37
  GSE34140 \[[@B54]\]      141            Affymetrix 250K Nsp         Yes/No                             141/0                        \-                                                          79/26/9/7                            \-                        10/107
  GSK \[[@B55]\]           10             Affymetrix 250K Nsp         Yes/No                             0/10                         \-                                                          \-                                   \-                        \-
  GSE19804 \[[@B56]\]      12             Affymetrix 6.0              Yes/No                             12/0                         \-                                                          6/3/2/1                              12/0                      12/0
  GSE25016 \[[@B34]\]      58             Affymetrix 6.0              Yes/No                             58/0                         \-                                                          \-                                   \-                        \-
  GSE33848 \[[@B9]\]       146            Affymetrix 6.0              Yes/No                             146/0                        \-                                                          95/0/0/0                             \-                        \-
  TCGA-AC \[[@B57]\]       135            Affymetrix 6.0              Yes/No                             135/0                        \-                                                          73/23/27/9                           78/57                     19/109
  GSE36363 \[[@B20]\]      62             Affymetrix 6.0 & 250K Nsp   Yes/Yes                            62/0                         5/14/24                                                     39/8/11/0                            37/25                     5/56
  Chitale \[[@B10]\]       184            Agilent 44K                 Yes/No                             184/0                        42/45/97                                                    121/26/32/5                          107/77                    39/145
  GSE20393 \[[@B58]\]      2              Agilent 244K                Yes/No                             2/0                          \-                                                          \-                                   \-                        \-
  E-TABM-926 \[[@B59]\]    17             Agilent 244K                Yes/No                             17/0                         13/1/3                                                      8/1/8/0                              14/3                      17/0
  E-TABM-1169 \[[@B59]\]   40             Illumina 370K               Yes/Yes                            40/0                         28/2/10                                                     22/5/13/0                            36/4                      40/0
  GSE31586                 5              ROMA 85K                    Yes/No                             0/5                          \-                                                          \-                                   \-                        \-
  TOTAL                    1272           \-                          1272/733                           1208/64                      113/134/210                                                 515/117/129/27                       433/276                   165/586

A: Number of used adenocarcinoma cases per cohort.B: Included in overall GISTIC analysis (CN) and/or GAP-analysis for EGFR/KRAS mutation groups.C: For tumors only.

D: NS = never-smoker, S = smoker

10.1371/journal.pone.0078614.t002

###### Clinical characteristics of patients with available mutation status in gene expression cohorts.

                                 **Discovery cohorts**   **Validation cohorts**                                                                                                                         
  ------------------------------ ----------------------- ------------------------ ----------------------- ----------------------- ----------------------- ------------- ----------------- ------------- -------------
  **Total number of patients**   226                     91                       102                     99                      106                     117           58                85            30
  **Gender**                                                                                                                                                                                            
  Male                           105                     41                       42                      15                      46                      60            13                37            20
  Female                         121                     50                       60                      84                      60                      57            45                48            10
  **Mutation status**                                                                                                                                                                                   
  *EGFR*-mutated                 127                     15                       24                      49                      18                      45            17                11            3
  *KRAS*-mutated                 20                      11                       36                      17                      43                      15            22                20            8
  EGFRwt/KRASwt                  79                      65                       42                      33                      45                      57            19                54            19
  **Stage**                                                                                                                                                                                             
  I                              168                     53                       70                      57                      70                      79            \-                48            \-
  II                             58                      20                       10                      10                      19                      13            \-                18            \-
  III                            0                       18                       17                      32                      13                      25            \-                11            \-
  IV                             0                       0                        5                       0                       4                       0             \-                3             \-
  **Usage**                                                                                                                                                                                             
  Differential gene expression   x                       x                        x                       x                       x                       x             x                 x             x
  Unsupervised analysis          x                       x                        x                       x                                                                                             
  **Platform**                   Affymetrix U133 2plus   Affymetrix U133A         Affymetrix U133 2plus   Affymetrix U133 2plus   Affymetrix U133 2plus   Agilent 44K   Illumina WG6 V3   Agilent 44K   Custom cDNA

These cohorts contain additional samples with unknown *EGFR* and *KRAS* mutation status.

EGFR and KRAS mutation analysis {#s2.2}
-------------------------------

*EGFR* and *KRAS* mutation status was determined as described in either [File S1](#pone.0078614.s005){ref-type="supplementary-material"} (for GSE37745 \[[@B25]\] and GSE28572 \[[@B27]\]) or in each of the original articles.

Genomic analyses {#s2.3}
----------------

Normalized copy number and B allele frequency estimates for Affymetrix microarrays and Illumina SNP beadchips, and normalized copy number estimates for Agilent 44K, Agilent 244K and ROMA 85K cohorts were generated and/or assembled as described in Staaf et al. \[[@B21]\] and [File S1](#pone.0078614.s005){ref-type="supplementary-material"}. Probe annotations for all array platforms were updated to the hg18/NCBI36 genome build. Genomic profiles were partitioned, centralized, and merged to a common probe set as described (\[[@B21]\] and [File S1](#pone.0078614.s005){ref-type="supplementary-material"}).

A modified version of Genomic Identification of Significant Targets in Cancer (GISTIC) \[[@B28]\], referred to as mGISTIC herein, was used for identification of focal copy number alterations and recurrent amplifications from the 1272-sample cohort (see \[[@B21]\] and [File S1](#pone.0078614.s005){ref-type="supplementary-material"}). Robustness of identified regions was assessed by permutation analysis ([Figure S1](#pone.0078614.s001){ref-type="supplementary-material"} and [File S1](#pone.0078614.s005){ref-type="supplementary-material"}). A genome-wide screen of differential copy number gain and loss between the three *EGFR* and *KRAS* defined mutation groups was performed by division of genomic profiles into 12,698 sequential segments of \~200 Kbp size, excluding reported regions of copy number variation. Each segment was subsequently tested for differences in frequency of copy number gain or loss individually. Fisher's exact test or the Chi-square test was used to identify genomic regions and recurrent amplifications with different frequency between mutation groups.

For tumors analyzed by SNP microarrays (n=141), B allele frequency estimates were partitioned \[[@B29]\], integrated with copy number data, and subjected to Genome Alteration Print (GAP) \[[@B30]\] analysis for estimation of allele-specific copy numbers and in silico tumor ploidy (referred to as GAP-ploidy herein) as described \[[@B21]\]. Loss-of-heterozygosity (LOH), copy-neutral LOH, and copy-neutral allelic imbalance were estimated from GAP results as described \[[@B21]\]. The fractions of the genome altered by copy number alterations, LOH, copy-neutral LOH, and copy-neutral allelic imbalance were calculated as described \[[@B21]\]. Data processing steps are further described in [File S1](#pone.0078614.s005){ref-type="supplementary-material"} and \[[@B21]\].

Gene expression analyses {#s2.4}
------------------------

Affymetrix cohorts were individually normalized using GC Robust Multi-array Averaging (GCRMA) \[[@B31]\]. For non-Affymetrix cohorts, normalized expression data were obtained from Gene Expression Omnibus \[[@B32]\]. In total, nine cohorts were analyzed individually for transcriptional differences between mutation groups as either discovery (n=5) or validation cohorts (n=4) ([Table 2](#pone-0078614-t002){ref-type="table"}). Differentially expressed genes between *EGFR*-mutated, *KRAS*-mutated and EGFRwt/KRASwt tumors were identified by ANOVA with false discovery rate adjustment using a 5% threshold for statistical significance. Hierarchical clustering was performed using Pearson correlation and complete linkage. Data processing steps are further described in [File S1](#pone.0078614.s005){ref-type="supplementary-material"}.

Results {#s3}
=======

Copy number alterations in lung adenocarcinoma {#s3.1}
----------------------------------------------

To identify copy number alterations (CNAs) of general importance in lung adenocarcinoma, which may serve as basis for supervised comparisons between *EGFR*-mutated, *KRAS*-mutated and EGFRwt/KRASwt tumors, we analyzed 1272 tumors and cell lines profiled by SNP or aCGH microarrays ([Figure 2A](#pone-0078614-g002){ref-type="fig"}, [Table 1](#pone-0078614-t001){ref-type="table"}). To pinpoint recurrent CNAs in lung adenocarcinoma, we performed an mGISTIC analysis of the entire 1272-sample set identifying 59 gains and 31 losses distributed across all autosomes ([Figure 2A](#pone-0078614-g002){ref-type="fig"}, [Table S1](#pone.0078614.s006){ref-type="supplementary-material"}). Several of the identified mGISTIC regions harbored known or putative adenocarcinoma driver candidates, such as *EGFR*, *MDM2*, *KRAS*, *MYC*, *TERT*, *MET*, *CCND1*, *NKX2-1/TITF1*, *CDK4*, *ERBB2*, *ID1*, *RB1*, *CDKN2A* and *PTEN*.

![Copy number alterations in lung adenocarcinoma.\
Frequency of copy number gain (red) and loss (green) for adenocarcinoma stratified by *EGFR* and *KRAS* mutational status using log~2~ratio ± 0.12 as threshold for identification of copy number gain and loss. Probes matched to known copy number variations are excluded. Black regions indicate genomic position of significant mGISTIC regions, which were identified from analysis of the entire 1272-sample cohort across chromosomes. Arrows indicate genomic regions with apparently different copy number alteration frequency between EGFR/KRAS mutation groups (6q, 7p, and 16p). (**A**) All 1272 adenocarcinomas. (**B**) 113 *EGFR*-mutated adenocarcinoma tumors. (**C**) 134 *KRAS*-mutated adenocarcinoma tumors. (**D**) 210 EGFRwt/KRASwt adenocarcinoma tumors.](pone.0078614.g002){#pone-0078614-g002}

Copy number alterations in EGFR/KRAS mutation groups {#s3.2}
----------------------------------------------------

Stratification of the 457 adenocarcinoma tumors with known *EGFR* and *KRAS* mutation status into *EGFR*-mutated (n=113), *KRAS*-mutated (n=134) and EGFRwt/KRASwt (n=210) tumors revealed both common alterations across mutation groups, such as gains of chromosome 1q and 8q, and loss of 3p, and regions with apparently different prevalence between mutation groups, including gains on chromosome 7p (*EGFR*-mutated) and 16p (*EGFR*-mutated), and losses on 6q (*KRAS*-mutated) ([Figures 2B-D](#pone-0078614-g002){ref-type="fig"}). In general, *EGFR*-mutated tumors displayed more copy number alterations (estimated by the fraction of the genome altered by CNA, CN-FGA) than non-*EGFR*-mutated tumors ([Figure 3A](#pone-0078614-g003){ref-type="fig"}). This pattern was consistent also in five out of six individual cohorts that included both *EGFR*-mutated and non-*EGFR*-mutated tumors. When the three mutation groups were stratified by clinicopathological variables, *EGFR*-mutated tumors continued to display higher CN-FGA fractions in stage I tumors, female patients, and never-smokers ([Figure 3A](#pone-0078614-g003){ref-type="fig"}).

Analysis of the 90 focal mGISTIC regions (derived from analysis of the 1272 adenocarcinoma sample cohort) identified 17 regions discriminating between the three mutation groups. 15 of these 17 regions showed the highest alteration frequency in *EGFR*-mutated tumors, while the remaining two regions showed highest frequency in *KRAS*-mutated tumors (Bonferroni adjusted Fisher's exact test p\< 0.05 and frequency difference \>20%, [Figure 3B](#pone-0078614-g003){ref-type="fig"} and [Table 3](#pone-0078614-t003){ref-type="table"}). Specifically, *EGFR*-mutated tumors showed higher frequencies of copy number gain on chromosomes 1p34.2 (including *MYCL*), 5p15.33, 5q35.1, 7p22.3-p22.2, 7p21.1, 7p11.2 (including *EGFR*), 7q11.21, 14q21.2, and 16p13.13, and copy number loss in regions at 8p (including *DUSP4*), 9p (including *CDKN2A*), and 10q23.2-q23.31 (*PTEN*). *KRAS*-mutated tumors showed higher frequencies of gain on 12p12.1 (*KRAS*) and loss at 6q16.3-q21.

![Copy number alterations and tumor ploidy in EGFR/KRAS mutation groups.\
(**A**) Pattern of gross copy number alterations measured as fraction of the genome altered by copy number gain or loss in adenocarcinoma tumors stratified by *EGFR* and *KRAS* mutation status (EGFR:red, KRAS:light blue, EGFRwt/KRASwt:gray), stage, gender and patient smoking status. Copy number alterations were called using log~2~ratio ± 0.12 as thresholds for identification of copy number gain and loss. P-values were calculated using ANOVA for indicated groups, \*\*\*: P\< 0.001, \*\*: P\< 0.01, \*: P\< 0.05. Top axis indicates number of cases per group. (**B**) mGISTIC regions discriminating between *EGFR*-mutated (red), *KRAS*-mutated (light blue) and EGFRwt/KRASwt (gray) adenocarcinoma tumors. mGISTIC regions identified by Fisher's exact test (Bonferroni adjusted p-value \< 0.05) with an additional requirement of \> 20% frequency difference between the lowest and highest groups. The y-axis describes the frequency of copy number gain or loss in respective group. (**C**) Distribution of GAP-ploidy across the adenocarcinoma EGFR/KRAS mutation groups for 141 tumors analyzed by GAP. A GAP-ploidy of two equals a diploid, three a triploid genome and four a tetraploid genome. Curves were generated by an Epanechnikov smoothing kernel with 0.1 smoothing bandwidth.](pone.0078614.g003){#pone-0078614-g003}

10.1371/journal.pone.0078614.t003

###### Differences and similarities in genomic alterations and allelic imbalances between EGFR/KRAS mutation groups.

  **Investigated property**                                                                     ***EGFR*-mutated**                                                                                                                                                 ***KRAS*-mutated**               **EGFRwt/KRASwt**
  --------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------------------------------------------------------------------------ -------------------------------- -------------------
  Fraction of the genome altered by copy number gain and loss                                   More ^a^                                                                                                                                                           Less                             Less
  Fraction of the genome altered by LOH                                                         More                                                                                                                                                               Less                             Less
  Fraction of the genome altered by copy number neutral LOH or copy-neutral allelic imbalance   Less                                                                                                                                                               Less                             Less
  Overall frequency of recurrent amplifications and genomic complexity                          More                                                                                                                                                               Less                             Less
  Amplifications differing between mutation groups (mGISTIC regions)                            7p11.2 (*EGFR*), 12q14.2-q14.3, 12q15 (*MDM2*)                                                                                                                                                      8p12 (*FGFR1*)
  Copy number alterations differing between mutation groups (mGISTIC regions)                   +1p34.2, +5p15.33, +5q35.1, +7p22.3-p22.2, +7p21.1, +7p11.2, +7q11.21, +14q21.2, +16p13.13, -8p23.2-p23.1, -8p21.2-p12, -9p24.3, -9p23, -9p21.3, -10q23.2-q23.31   +12p12.1 (*KRAS*), -6q16.3-q21   
  Copy number alterations differing between mutation groups (genome-wide screen)                +1p36.33-p31.1, +5q34-q35.3, +7p22.3-p11.1, +7q11.21, +16p13.3-p12.1, +16p11.2-q11.1, +21q22.11-q22.3, -8p22-p11.21, -8q11.23                                                                       
  Characteristic total LOH regions (≥50% frequency)^b^                                          8p, 9, 13q, 17p                                                                                                                                                    19p                              
  Predominant tumor ploidy (GAP-ploidy)                                                         2N (highest), 3N                                                                                                                                                   2N (highest), 3N                 2N (highest), 3N

a *Les*s indicates relatively lower estimates or frequencies between groups. *More* indicates relatively higher estimates or frequencies.

b Includes LOH caused by copy number loss as well as copy-neutral LOH.

A genome-wide analysis of differences in copy number frequency between the three mutation groups identified nine large coherent genomic regions (seven gains and two losses), all with higher alteration frequency in *EGFR*-mutated tumors. Regions were located on 1p, 5q, 7p, 7q, 8p, 8q, 16p and 21q, and involved 8% (7% gain, 1% loss) of the analyzed genome (Hochberg adjusted Fisher's exact test p\<0.01 and minimum frequency difference \>20%, [Tables 3](#pone-0078614-t003){ref-type="table"} and [S2](#pone.0078614.s007){ref-type="supplementary-material"}).

Similar to copy number gain and loss in general, *EGFR*-mutated tumors also displayed more recurrent amplifications in the 59 mGISTIC regions of gain compared with the non-*EGFR*-mutated tumors (p=0.004, Chi-square test). This finding was consistent also in patients with stage I disease (p=0.02, Chi-square test) or female gender (p=0.004, Chi-square test). In higher stage (≥II) tumors and in male patients the *EGFR*-mutated group also showed more recurrent amplifications, however not reaching statistical significance due to the lower number of tumors in these comparisons. In exploratory analysis, individual recurrent amplifications at 7p11.2 (*EGFR*), 8p12 (*WHSC1L1*, *FGFR1*), and 12q14-q15 (including *MDM2*) discriminated between mutation groups (p\<0.05, Fisher's exact test, [Table 3](#pone-0078614-t003){ref-type="table"}).

Taken together, these results suggest a higher genomic complexity in *EGFR*-mutated adenocarcinomas compared with *KRAS*-mutated and EGFRwt/KRASwt tumors.

Patterns of tumor ploidy and allelic imbalance in EGFR/KRAS mutation groups {#s3.3}
---------------------------------------------------------------------------

Patterns of tumor ploidy and allelic imbalances between mutational groups were evaluated by GAP \[[@B30]\] analysis of 141 tumors (n=43 *EGFR*-mutated, 39 *KRAS*-mutated, and 59 EGFRwt/KRASwt) profiled by SNP microarrays. Primarily, no differences in distribution of tumor ploidy (estimated by GAP-ploidy) were observed between mutation groups (p=0.96, ANOVA, [Figure 3C](#pone-0078614-g003){ref-type="fig"}). Secondly, *EGFR*-mutated adenocarcinomas were weakly associated with higher fractions of total LOH compared with *KRAS*-mutated and EGFRwt/KRASwt tumors both overall and in stage I disease (p=0.05, ANOVA), but not in tumors of higher stages (≥II), or tumors stratified by gender. In contrast, no significant differences in copy-neutral LOH or copy-neutral allelic imbalance fractions were observed between mutation groups overall or when sub-stratified by stage or gender. The highest frequencies of total LOH (\>50%) were most often found in regions of copy number loss, while copy-neutral LOH and copy-neutral allelic imbalance showed an overall lower prevalence across chromosomes in all mutation groups (generally ≤10-15% frequency for copy-neutral LOH, and \<25% for copy-neutral allelic imbalance across chromosomes) ([Table 3](#pone-0078614-t003){ref-type="table"} and [Figure S2](#pone.0078614.s002){ref-type="supplementary-material"}).

Taken together, this implies that the weak associations of differences in allelic imbalances between the mutation groups are predominantly related to LOH caused by copy number loss in *EGFR*-mutated tumors.

Supervised and unsupervised analysis of transcriptional differences between EGFR/KRAS mutation groups {#s3.4}
-----------------------------------------------------------------------------------------------------

To identify a robust set of differentially expressed genes between the three mutation groups we performed supervised analysis of five adenocarcinoma Affymetrix cohorts (n=624 tumors, discovery cohorts, [Table 2](#pone-0078614-t002){ref-type="table"}). 96 genes showed consistent differential expression in ≥4 cohorts, while only 21 genes were differentially expressed across all five cohorts ([Tables 4](#pone-0078614-t004){ref-type="table"} and [S3](#pone.0078614.s008){ref-type="supplementary-material"}). We validated the 96 identified genes in four independent adenocarcinoma cohorts analyzed by different microarray platforms (n=290 tumors, [Table 2](#pone-0078614-t002){ref-type="table"}). In the independent cohorts, 41-96% of the 96 genes were present and thus available for further comparisons. Of the available genes 46-67% showed differential expression between the mutation groups in the independent cohorts (p\<0.05 ANOVA, [Table S3](#pone.0078614.s008){ref-type="supplementary-material"}).

10.1371/journal.pone.0078614.t004

###### Differentially expressed genes between *EGFR*/*KRAS* mutation groups in ≥4 of five Affymetrix adenocarcinoma cohorts.

  Gene       Name                                               Gene           Name
  ---------- -------------------------------------------------- -------------- --------------------------------------------------------
  ACSF2      acyl-CoA synthetase family member 2                ISG20          interferon stimulated exonuclease gene 20kDa
  ADCY9      adenylate cyclase 9                                ITPR3          inositol 1,4,5-triphosphate receptor, type 3
  AGFG1      *ArfGAP with FG repeats 1*                         KCNK5 \*       potassium channel, subfamily K, member 5
  AHR        aryl hydrocarbon receptor                          KIAA0319L \*   polycystic kidney disease 1-like
  APOH       apolipoprotein H precursor                         KIAA0494       hypothetical protein LOC9813
  ARMCX6     armadillo repeat containing, X-linked 6            KIAA0495       hypothetical protein LOC57212
  ARSD       arylsulfatase D                                    KIAA1033 \*    hypothetical protein LOC23325
  BAG1       BCL2-associated athanogene 1                       KRAS           c-K-ras2 protein isoform b
  BLVRA      biliverdin reductase A                             LDLRAP1 \*     low density lipoprotein receptor adaptor protein
  C16orf58   hypothetical protein LOC64755                      LRRC31         leucine rich repeat containing 31
  C7orf23    chromosome 7 open reading frame 23                 MANBA          mannosidase, beta A, lysosomal
  CADPS2     Ca2+-dependent activator protein for secretion 2   MEAF6 \*       MYST/Esa1-associated factor 6
  CAMTA1     calmodulin-binding transcription activator 1       MMP15          matrix metalloproteinase 15 preproprotein
  CLDN10     claudin 10                                         MTPAP          mitochondrial poly(A) polymerase
  COL21A1    collagen, type XXI, alpha 1 precursor              MYST1          MYST histone acetyltransferase 1
  CTNNBIP1   catenin, beta interacting protein 1                NAT15          N(alpha)-acetyltransferase 60, NatF catalytic subunit
  DDAH1 \*   dimethylarginine dimethylaminohydrolase 1          NBPF10         hypothetical protein LOC440673
  DDX21      DEAD (Asp-Glu-Ala-Asp) box polypeptide 21          NFYC \*        nuclear transcription factor Y, gamma
  DNAJC9     DnaJ homolog, subfamily C, member 9                NIPAL3         NIPA-like domain containing 3
  DUSP4 \*   dual specificity phosphatase 4                     PDSS1          prenyl diphosphate synthase, subunit 1
  EFHC2      EF-hand domain (C-terminal) containing 2           PEF1           penta-EF-hand domain containing 1
  EGFR \*    epidermal growth factor receptor                   PER3           period 3
  ELN        elastin                                            PIGV \*        phosphatidylinositol glycan class V
  ENC1       ectodermal-neural cortex (with BTB-like domain)    PIK3IP1        HGFL protein
  ENTPD4     ectonucleoside triphosphate diphosphohydrolase     PPCS           phosphopantothenoylcysteine synthetase isoform
  ETV5       ets variant gene 5 (ets-related molecule)          PPFIBP2        PTPRF interacting protein, binding protein 2
  FAAH       fatty acid amide hydrolase                         PPIF           peptidylprolyl isomerase F precursor
  FAM184A    family with sequence similarity 184, member A      PRDM4          PR domain containing 4
  FGF13      fibroblast growth factor 13                        PYROXD1        pyridine nucleotide-disulphide oxidoreductase domain 1
  FGG        fibrinogen, gamma chain                            RAPGEF5        Rap guanine nucleotide exchange factor (GEF) 5
  FGGY       FGGY carbohydrate kinase domain containing         RFK            riboflavin kinase

A three-group centroid classifier was used to explore the predictive power of the 96 genes in calling true mutation status ([File S1](#pone.0078614.s005){ref-type="supplementary-material"} and [Table S3](#pone.0078614.s008){ref-type="supplementary-material"}). Classification of the four independent cohorts showed an overall accuracy of 40-90% in classification across a range of classification cut-offs ([Figure S3A](#pone.0078614.s003){ref-type="supplementary-material"}). Sensitivity was highest in classification of *EGFR*-mutated tumors across the different cohorts (80-100%), followed by *KRAS*-mutated tumors ([Figure 4A](#pone-0078614-g004){ref-type="fig"}). However, specificities of the 96-gene classifier were lower (60-90%) for the *EGFR* and *KRAS*-mutated groups ([Figure 4A](#pone-0078614-g004){ref-type="fig"}). In contrast, for EGFRwt/KRASwt tumors, sensitivity was poor (10-60%) but specificity higher (80-100%).

![Supervised classification of adenocarcinoma gene expression cohorts with respect to *EGFR* and *KRAS* mutation status.\
Sensitivity (solid line) and specificity (dashed line) by EGFR/KRAS mutation group for classification of four independent validation cohorts using a 96-gene centroid classifier. The x-axis shows Pearson correlation cut-off for assigning a sample to the centroid with the highest correlation. Increasing correlation cut-offs correspond to increased stringency in classification, but introduces growing numbers of unclassified samples excluded in the calculation of sensitivity and specificity.](pone.0078614.g004){#pone-0078614-g004}

To further analyze the transcriptional patterns between the three mutation groups we performed unsupervised hierarchical clustering of four Affymetrix discovery cohorts ([Table 2](#pone-0078614-t002){ref-type="table"}). Each cohort was individually clustered using a) a signature of genes overexpressed in EGFRwt/KRASwt tumors \[[@B26]\], b) a *KRAS* dependency gene signature \[[@B33]\], and c) probe sets derived from three different expression variance filters reflecting at different stringency the variation in expression across all tumors in a cohort (Affymetrix probe set range n=1356-24052). In none of these analyses did clustering resolve the three mutation groups into discrete transcriptional groups without notable inclusion of tumors from other mutation groups ([Figure S4](#pone.0078614.s004){ref-type="supplementary-material"}). However, supportive of results from the supervised analyses we found that *EGFR*-mutated adenocarcinomas in general appeared to display a more distinctive expression pattern with enrichment of *EGFR*-mutated tumors (\~60% of all mutations) in specific clusters. In contrast, *KRAS*-mutated and EGFRwt/KRASwt tumors often appeared more intermixed, even when clustered using the *KRAS* dependency gene signature \[[@B33]\] ([Figure S4](#pone.0078614.s004){ref-type="supplementary-material"}).

Taken together, results from the supervised and unsupervised gene expression analyses suggest that mutation status is not translated into a clearly distinctive and prominent expression signature.

Discussion {#s4}
==========

In the current study we delineate genomic and transcriptional alterations in lung adenocarcinoma stratified by *EGFR* and *KRAS* mutation status. We show that a few specific copy number and transcriptional alterations exist between the three mutational groups, but also a considerable similarity caused by high intra-group heterogeneity and/or less distinctive inter-group differences. Together, this suggests that the overall genomic and transcriptional landscape of adenocarcinoma is affected, but only to a minor extent, by the mutational status of *EGFR* and *KRAS*.

Stratification of genomic profiles from 457 tumors with available *EGFR* and *KRAS* mutation status into three mutation groups revealed differences in the overall pattern of CNAs, amplifications and genomic architecture, as well as specific regions and amplifications differing in frequency between the groups (summarized in [Table 3](#pone-0078614-t003){ref-type="table"}). Overall, *EGFR*-mutated tumors displayed more CNAs, more amplifications, and higher genomic complexity than non-*EGFR*-mutated tumors consistent with previous reports \[[@B11],[@B15],[@B18]\]. Specific patterns of recurrent amplifications in between the mutation groups, such as 8p12 (harboring *FGFR1*) in EGFRwt/KRASwt and 12q amplifications (including the p53 repressor *MDM2*) in *EGFR*-mutated tumors were observed. *FGFR1* mutations are rarely observed in NSCLC, while *FGFR1* amplification is frequent in, e.g., squamous cell lung carcinoma and associated with increased protein levels and a FGFR1 proliferation dependency \[[@B34]\]. Moreover, FGF-FGFR pathway activation has been suggested to be one mediator of resistance to EGFR inhibitors, together with, e.g., *MET* amplification (see \[[@B35]\] for review and \[[@B36]\]). In the current study, *FGFR1* and *MET* amplifications were restricted to the EGFRwt/KRASwt tumor group, and were mutually exclusive (*MET* amplification was borderline non-significant for difference in frequency between mutation groups, p=0.09, Fisher's exact test). Together, this could indicate presence of specific genomic circuits acting as driving forces in pathogenesis in the different mutation groups.

Taken together, the analyses of differential genomic regions point to only a few, variably sized, regions with moderate frequency differences (20-40%) between the mutational groups. These regions predominantly include regions of copy number gain with higher frequency in *EGFR*-mutated tumors. Several of the regions have been reported previously, but as larger and less defined regions \[[@B9],[@B11]-[@B13],[@B17],[@B18]\], while others such as 5q34-q35.3, appear novel (see [Table S4](#pone.0078614.s009){ref-type="supplementary-material"} for literature comparison of 34 previously reported regions from five independent studies \[[@B9],[@B11]-[@B13],[@B18]\]). For instance, 17 of our mGISTIC regions were present in 34 previously reported regions differing between *EGFR*-mutated and *EGFR*-wild type tumors, or *KRAS*-mutated tumors and *KRAS*-wild type tumors ([Table S4](#pone.0078614.s009){ref-type="supplementary-material"}). Although 24 of the 34 reported regions showed statistical significance for the original comparisons in our cohort, only 16 of these 24 regions also showed \>20% frequency difference between the three groups. The absolute majority of these regions (88%) were located on chromosome 1p, 7p, and 16p (gains) and 8p (losses). Together, this emphasizes the need for adequately sized cohorts in order to draw reproducible conclusions when only moderate differences exist between investigated groups.

Few genome-wide analyses of differential allelic imbalance between *EGFR*-mutated, *KRAS*-mutated, and EGFRwt/KRASwt tumors exist in the literature. Blons et al. reported that *EGFR*-mutated tumors in general display more fractional allelic loss than *KRAS*-mutated tumors \[[@B18]\], consistent with our observation of higher fractions of total LOH in *EGFR*-mutated tumors. Moreover, Nakanishi et al. reported that two regions, 4q13 and 4q22, differ in allelic imbalance between the mutation groups \[[@B37]\]. However, in the current study we could not verify this finding using total LOH, copy-neutral LOH or copy-neutral allelic imbalance as measurements. Instead, we found that the frequency of total LOH was strongly correlated with regions of copy number loss. In contrast, the less frequent copy-neutral LOH and copy-neutral allelic imbalance events were overall more evenly distributed across chromosomes in the mutation groups. These findings are consistent with results for lung cancer histology groups in general \[[@B21]\], and also recent reports from breast cancer \[[@B38],[@B39]\], suggesting that LOH is predominantly caused by copy number loss in these tumor types. Moreover, the similarity in the observed tumor ploidy patterns appears consistent with that the mutation groups do not exhibit gross differences in either CNAs or allelic imbalances. However, it should be noted that the analyses of allelic imbalances in the current study are based on a smaller subset of samples (n=141), which could be a source of variability.

Several studies have reported *KRAS* mutant signatures or differentially expressed genes between adenocarcinomas with *EGFR* and/or *KRAS* mutations and respective wild type cases \[[@B8],[@B10],[@B12],[@B14],[@B16],[@B26],[@B33],[@B40]\]. However, the overlap between these public signatures is very poor when directly compared ([Figures S3B and C](#pone.0078614.s003){ref-type="supplementary-material"}). We identified 96 differentially expressed genes by supervised gene expression analyses between mutation groups across multiple discovery cohorts, of which several could be validated in independent cohorts. Reasons for the lower number of significant genes in the independent cohorts likely include smaller sample sizes and different microarray platforms compared with the discovery cohorts. The low number of differentially expressed genes between the mutation groups (only 21 genes consistently differentially expressed in all five discovery cohorts) is similar to results from other studies \[[@B8],[@B12],[@B14],[@B16]\]. This low number of differentially expressed genes argues against that the mutational subgroups represent distinct transcriptional groups. Moreover, the overlaps between our 96 genes and previous studies \[[@B8],[@B14],[@B16],[@B26],[@B33],[@B40]\] were poor (1-5% individual overlap between signatures). These results underline the need for a multicohort approach for identification of robust transcriptional differences between the mutation groups. Notably, our 96 genes mapped to a higher extent (43% of genes) to genomic regions showing differences in frequency of copy number gain or loss between the mutation groups compared to gene signatures from five reported studies \[[@B8],[@B12],[@B14],[@B16],[@B26]\] (6-15% of reported genes). However, the influence of the modest differences in CNA frequency (20-40%) between the mutation groups on transcriptional levels is difficult to assess. In addition to *EGFR* and *KRAS*, differentially expressed genes between the mutation groups included several other genes reported to be involved in tumorigenesis (*DUSP4*, *RPS6KA1*, *ID1*, *TNFRSF10B*, *CAMTA1*) \[[@B10],[@B41]\], and, consistent with the enrichment of never-smokers in the *EGFR*-mutated patient group, genes reported as deregulated by smoking (*AHR*, *CLDN10*, *FGG*, *GGA2*, *GUSB*, *TXNRD1*) \[[@B42]-[@B45]\].

In supervised classification, the 96 differentially expressed genes identified *EGFR*-mutated adenocarcinomas with high sensitivity, but poorer specificity, while opposite results was found for EGFRwt/KRASwt tumors. Together with the results from unsupervised hierarchical clustering of multiple gene expression cohorts using different gene or probe sets these analyses demonstrate the difficulty in separating the mutation groups, especially *KRAS*-mutated and EGFRwt/KRASwt tumors, into more discrete transcriptional entities. Chitale et al. \[[@B10]\] proposed that the more distinctive expression pattern of *EGFR*-mutated tumors compared to *KRAS*-mutated tumors may depend on either a less prominent effect of *KRAS* mutations on expression, a biological or etiological heterogeneity among *KRAS*-mutated tumors, or that *EGFR* mutations arise in a more homogeneous and restricted cell type. Our results may be interpreted as support for potentially all three hypotheses, given the differences observed between and within mutation groups. Together, the results from our supervised and unsupervised gene expression analyses suggest that only modest, reproducible, transcriptional differences exist between the mutation groups. This conclusion appears consistent with the somewhat mixed inclusion of *EGFR*-mutated, *KRAS*-mutated and EGFRwt/KRASwt adenocarcinomas in different reported molecular subtypes of adenocarcinomas \[[@B16],[@B19],[@B20]\]. Although the bronchioid molecular subtype \[[@B19]\] has been strongly associated with *EGFR*-mutated tumors, this subtype also includes notable fractions of *KRAS*-mutated and EGFRwt/KRASwt tumors (see, e.g., \[[@B19],[@B20]\]). Moreover, \~30% or more of *EGFR*-mutated have been classified as non-bronchioid (magnoid or squamoid) in discovery cohorts in previous studies \[[@B19],[@B20],[@B46]\]. In the absence of bronchioid classified tumors we found no significant association between the magnoid and squamoid subtypes and EGFR/KRAS mutation status in any of the five discovery cohorts in the current study (data not shown). These findings appear consistent with our unsupervised analysis showing a more distinct expression pattern of a subset of *EGFR*-mutated tumors across multiple cohorts, while the *KRAS*-mutated and EGFRwt/KRASwt groups are more intermixed ([Figure S4](#pone.0078614.s004){ref-type="supplementary-material"}). These results also suggest that *EGFR*-mutated tumors could be divided into additional subgroups, which we have recently demonstrated \[[@B46]\]. Taken together, *EGFR* and *KRAS* mutational status do not appear to be translated into a clearly distinctive and prominent expression signature in lung adenocarcinoma.

To further delineate the observed heterogeneous patterns of CNAs, allelic imbalances and gene expression patterns in the three mutational groups identification and/or definition of new molecular subgroups within the *EGFR*-mutated, *KRAS*-mutated and EGFRwt/KRASwt tumor groups are needed. For instance, although EGFRwt/KRASwt adenocarcinomas with *ALK* rearrangements are reported to display distinct expression profiles compared with *ALK*-negative tumors \[[@B26]\], it remains unclear whether this is also true for CNAs and allelic imbalances. Recent studies of lung adenocarcinoma have suggested that molecular profiling could be of value in future clinical decision making by providing clues about, e.g., treatment response to EGFR inhibitors \[[@B17],[@B20],[@B47],[@B48]\]. For instance, Yuan et al. recently reported that clustered CNAs (copy number gains) on chromosome 7p were associated with poorer survival and less favorable response to EGFR tyrosine kinase inhibitors in *EGFR*-mutated adenocarcinomas specifically \[[@B17]\]. In support of Yuan et al., we recently identified a gene signature associated with poorer survival for patients with *EGFR*-mutated adenocarcinomas, where the high-risk patient group showed more copy number gains and amplifications on chromosome 7p \[[@B46]\]. As regions on chromosome 7p display some of the largest frequency differences between the mutational groups (\~40%) these findings highlight the need for a more detailed characterization of this chromosome arm. The growing number of detected tyrosine kinase fusions in predominantly EGFRwt/KRASwt adenocarcinomas (including *ALK*, *RET*, and *ROS1*) are also becoming increasingly important in the therapeutic setting, as these alterations are/may become targets for specialized molecular agents. However, it remains to be investigated whether there exist similar regions and/or gene signatures associated with treatment response also for these adenocarcinoma subgroups. Clearly, further molecular stratification within the *EGFR* and *KRAS* mutation-defined lung adenocarcinoma groups has the potential to reveal new targets for synergistic treatment and provide insights into resistance mechanisms.

In summary, our multicohort analyses of genomic and transcriptional alterations demonstrate both differences and strong similarities between the *EGFR* and *KRAS* mutation defined adenocarcinoma groups. Moreover, our results suggest that the overall genomic and transcriptional landscape of adenocarcinoma is only to a minor extent affected by the mutational status of *EGFR* and *KRAS*.

Supporting Information
======================

###### 

**Permutation analysis of mGISTIC regions.** Close support from permutation analysis for an mGISTIC region (n=90) is defined as the % of times the region was enclosed or overlapped by a permuted region based on a 75% sample subset of the 1272 samples (n=100 permutations). (**A**) Cumulative fraction of regions (all, gain, loss) stratified into bins of 10% close support. (**B**) Hexagonal binning of mGISTIC regions (all, gain, loss) for % of close support versus -log~10~(p-value) of detected regions. A general trend of higher p-values connected to lower % close support is observed. Colors of bins indicate number of regions. Taken together, regions showing the lowest permutation detection rates also showed the lowest g-scores \[[@B28]\] and p-values. This is consistent with that these regions are present in only a small subset of the 1272 cases, which makes the identification and delineation of these regions to sensitive to sample composition.

(PDF)

###### 

Click here for additional data file.

###### 

**Pattern of CNAs, LOH, CNN-LOH, and CNN-AI in EGFR/KRAS mutation groups.** Panels show in decreasing order from the top pattern (frequency) of copy number gain (red) and loss (green) relative to GAP-ploidy with mGISTIC regions identified from the 1272 sample cohort indicated by blue dots, LOH, copy-neutral LOH (CNN-LOH), copy-neutral allelic imbalance (CNN-AI), and variation of FGA values versus GAP-ploidy for copy number (black), CNN-AI (red), LOH (blue), and CNN-LOH (light blue) in the bottom panel. For the bottom panel GAP-ploidy estimates were binned in bins of size 0.3, which is represented by tick marks on the x-axis. For each bin the median FGA value of the included samples is plotted (points) for copy number, LOH, CNN-LOH and CNN-AI. Bins contain different numbers of samples (top axis). The 141 tumors with mutation status analyzed by GAP were stratified into (**A**) *EGFR*-mutated (n=43), (**B**) *KRAS*-mutated (n=39), and (**C**) EGFRwt/KRASwt (n=59) tumors.

(PDF)

###### 

Click here for additional data file.

###### 

**Comparison of public EGFR/KRAS signatures and classification by a set of genes differentially expressed between EGFR/KRAS mutation groups across multiple cohorts.** (**A**) Overall accuracy for classification of four independent adenocarcinoma cohorts using a 96-gene centroid classifier. The number of genes in the centroid matching to the different cohorts varies. The x-axis shows Pearson correlation cut-off for assigning a sample to the centroid with the highest correlation. Increasing correlation cut-offs introduces growing numbers of unclassified samples, which are excluded in calculation of accuracy. (**B**) Venn-diagram of the gene overlap between four reported gene lists of differentially expressed genes between *EGFR*-mutated and *EGFR*-wild type adenocarcinoma tumors. (**C**) Venn-diagram of the gene overlap between two reported gene lists of differentially expressed genes between *KRAS*-mutated and *KRAS*-wild type adenocarcinoma tumors, and two reported *KRAS* mutant signatures \[[@B33],[@B40]\].

(PDF)

###### 

Click here for additional data file.

###### 

**Unsupervised analyses of four Affymetrix adenocarcinoma gene expression cohorts using different probe sets.** Unsupervised hierarchical clustering was performed using Pearson correlation and complete linkage using five different probe sets in four adenocarcinoma cohorts that were analyzed by Affymetrix gene expression microarrays. Dendrograms for each cluster tree were cut into the top two or three clusters, and the number of probe sets used in the clustering is shown for each cohort. For each cohort the distribution of *EGFR*-mutated (red), *KRAS*-mutated (blue), and EGFRwt/KRASwt (black) tumors are shown across clusters as bars. Percentages in bar plots correspond to, e.g., how many *EGFR*-mutated tumors of the total number of *EGFR*-mutated cases that reside in a particular cluster. (**A**) Clustering based on probe sets from a list of 190 probe sets reported to be upregulated in EGFRwt/KRASwt adenocarcinomas \[[@B26]\]. (**B**) Clustering based on probe sets with log2ratio standard deviation \>0.3 across tumors in a cohort. (**C**) Clustering based on probe sets with log2ratio standard deviation \>0.5 across tumors in a cohort. (**D**) Clustering based on probe sets with log2ratio standard deviation \>1 across tumors in a cohort. (**E**) Clustering based on matching genes from the list of top 250 genes reported by Singh et al. \[[@B33]\]. Each dendrogram is cut into the top two clusters. Division of dendrograms into three groups did not identify *KRAS*-mutants as a single group without notable inclusion of EGFRwt/KRASwt tumors in all cohorts.

(PDF)

###### 

Click here for additional data file.

###### 

**Document with details concerning used analysis methods.**

(DOC)

###### 

Click here for additional data file.

###### 

**Genomic mGISTIC regions identified from analysis of 1272 lung adenocarcinomas.**

(XLSX)

###### 

Click here for additional data file.

###### 

**Differential regions of copy number gain and loss between mutation groups obtained from genome-wide analysis.**

(XLSX)

###### 

Click here for additional data file.

###### 

**Differentially expressed genes between mutation groups across at least four gene expression cohorts.**

(XLSX)

###### 

Click here for additional data file.

###### 

**Analysis of genomic regions reported in the literature to stratify mutation groups in the current cohort.**

(DOC)

###### 

Click here for additional data file.

[^1]: **Competing Interests:**The authors have declared that no competing interests exist.

[^2]: Conceived and designed the experiments: JS MP SI. Performed the experiments: KE. Analyzed the data: JS. Contributed reagents/materials/analysis tools: JS. Wrote the manuscript: JS MP. Contributed patient data: KE PM JB.
