Genomic and transcriptional alterations in lung adenocarcinoma in relation to EGFR and KRAS mutation status.

INTRODUCTION: In lung adenocarcinoma, the mutational spectrum is dominated by EGFR and KRAS mutations. Improved knowledge about genomic and transcriptional alterations in and between mutation-defined subgroups may identify genes involved in disease development or progression.
METHODS: Genomic profiles from 457 adenocarcinomas, including 113 EGFR-mutated, 134 KRAS-mutated and 210 EGFR and KRAS-wild type tumors (EGFRwt/KRASwt), and gene expression profiles from 914 adenocarcinomas, including 309 EGFR-mutated, 192 KRAS-mutated, and 413 EGFRwt/KRASwt tumors, were assembled from different repositories. Genomic and transcriptional differences between the three mutational groups were analyzed by both supervised and unsupervised methods.
RESULTS: EGFR-mutated adenocarcinomas displayed a larger number of copy number alterations and recurrent amplifications, a higher fraction of total loss-of-heterozygosity, higher genomic complexity, and a more distinct expression pattern than EGFR-wild type adenocarcinomas. Several of these differences were also consistent when the three mutational groups were stratified by stage, gender and smoking status. Specific copy number alterations were associated with mutation status, predominantly including regions of gain with the highest frequency in EGFR-mutated tumors. Differential regions included both large and small regions of gain on 1p, 5q34-q35.3, 7p, 7q11.21, 12p12.1, 16p, and 21q, and losses on 6q16.3-q21, 8p, and 9p, with 20-40% frequency differences between the mutational groups. Supervised gene expression analyses identified 96 consistently differentially expressed genes between the mutational groups, and together with unsupervised analyses these analyses highlighted the difficulty in broadly resolving the three mutational groups into distinct transcriptional entities.
CONCLUSIONS: We provide a comprehensive overview of the genomic and transcriptional landscape in lung adenocarcinoma stratified by EGFR and KRAS mutations. Our analyses suggest that the overall genomic and transcriptional landscape of lung adenocarcinoma is affected, but only to a minor extent, by EGFR and KRAS mutation status.

Introduction

Lung cancer is a heterogeneous malignancy with poor survival due to diagnosis at an often advanced stage [1]. 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 [2]. In adenocarcinoma, the mutational spectrum is dominated by EGFR and KRAS mutations, where the former is an established predictor of response to EGFR inhibitors [3,4]. EGFR and KRAS mutations are nearly always mutually exclusive and associated with differences in patient gender and smoking history [5]. Together, this suggests that these genetic alterations may be drivers of pathogenesis for specific adenocarcinoma subgroups [5] (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 [5-7]. Several studies have reported genomic or transcriptional alterations between EGFR-mutated and/or KRAS-mutated tumors and corresponding wild-type adenocarcinomas [8-18]. 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. [19]. 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 [19,20]. 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 [19,20].

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). 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.

Figure 1

Schematic diagram of genomic and transcriptional analyses.

(A) Genomic analyses. (B) Transcriptional analyses. Individual cohorts are portrayed.

Materials and Methods

Tumor material

Genomic profiles from 1272 adenocarcinoma tumors and cell lines were obtained from a previous study (n=1210) [21], with addition of adenocarcinomas from Wilkerson et al. [20] (n=62, GSE36363). All genomic profiles were analyzed in an unmatched fashion and sample uniqueness was assured as described [21]. 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 [8,10,20,22-26]. Samples from Chitale et al. [10] 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 and 2, with additional description in File S1.

Table 1

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

Genomic cohort	No. of ACA	Microarray platform	Included in CN / GAP analysisB	No. of tumors / cell lines	No. of EGFR-mutated / KRAS-mutated / EGFRwt/KRASwtC	No. of stage I/II/III/IV tumorsC	Gender female / maleC	Smoking status NS/SC,D
Zhao [49]	36	Affymetrix 100K	Yes/No	36/0	-	-	2/1	-
GSE18252 [50]	4	Affymetrix 100K	Yes/Yes	4/0	0/4/0	-	-	-
Weir [51]	112	Affymetrix 250K Sty	Yes/No	112/0	-	-	-	-
Weir [51]	196	Affymetrix 250K Sty	Yes/No	196/0	15/49/51	41/14/15/3	110/83	17/116
GSE19399 [52]	19	Affymetrix 250K Sty	Yes/Yes	19/0	1/1/8	10/1/4/0	10/9	2/16
GSE17247 [53]	49	Affymetrix 250K Sty	Yes/No	0/49	-	-	-	-
GSE28572 [27]	44	Affymetrix 250K Nsp	Yes/Yes	44/0	9/18/17	21/10/8/2	27/17	4/37
GSE34140 [54]	141	Affymetrix 250K Nsp	Yes/No	141/0	-	79/26/9/7	-	10/107
GSK [55]	10	Affymetrix 250K Nsp	Yes/No	0/10	-	-	-	-
GSE19804 [56]	12	Affymetrix 6.0	Yes/No	12/0	-	6/3/2/1	12/0	12/0
GSE25016 [34]	58	Affymetrix 6.0	Yes/No	58/0	-	-	-	-
GSE33848 [9]	146	Affymetrix 6.0	Yes/No	146/0	-	95/0/0/0	-	-
TCGA-AC [57]	135	Affymetrix 6.0	Yes/No	135/0	-	73/23/27/9	78/57	19/109
GSE36363 [20]	62	Affymetrix 6.0 & 250K Nsp	Yes/Yes	62/0	5/14/24	39/8/11/0	37/25	5/56
Chitale [10]	184	Agilent 44K	Yes/No	184/0	42/45/97	121/26/32/5	107/77	39/145
GSE20393 [58]	2	Agilent 244K	Yes/No	2/0	-	-	-	-
E-TABM-926 [59]	17	Agilent 244K	Yes/No	17/0	13/1/3	8/1/8/0	14/3	17/0
E-TABM-1169 [59]	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

Table 2

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

	Discovery cohorts					Validation cohorts
	GSE31210 [26]	Chitale U133A [10]	Chitale U133 2plus [10]	E-MTAB-923 [24]*	GSE37745 [25]	GSE13213 [22]	GSE32863 [23]	GSE26939 [20]*	GSE8569 [8]
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

EGFR and KRAS mutation status was determined as described in either File S1 (for GSE37745 [25] and GSE28572 [27]) or in each of the original articles.

Genomic analyses

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. [21] and File S1. 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 ([21] and File S1).

A modified version of Genomic Identification of Significant Targets in Cancer (GISTIC) [28], referred to as mGISTIC herein, was used for identification of focal copy number alterations and recurrent amplifications from the 1272-sample cohort (see [21] and File S1). Robustness of identified regions was assessed by permutation analysis (Figure S1 and File S1). 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 [29], integrated with copy number data, and subjected to Genome Alteration Print (GAP) [30] analysis for estimation of allele-specific copy numbers and in silico tumor ploidy (referred to as GAP-ploidy herein) as described [21]. Loss-of-heterozygosity (LOH), copy-neutral LOH, and copy-neutral allelic imbalance were estimated from GAP results as described [21]. The fractions of the genome altered by copy number alterations, LOH, copy-neutral LOH, and copy-neutral allelic imbalance were calculated as described [21]. Data processing steps are further described in File S1 and [21].

Gene expression analyses

Affymetrix cohorts were individually normalized using GC Robust Multi-array Averaging (GCRMA) [31]. For non-Affymetrix cohorts, normalized expression data were obtained from Gene Expression Omnibus [32]. 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). 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.

Results

Copy number alterations in lung adenocarcinoma

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, Table 1). 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, Table S1). 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.

Figure 2

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 log2ratio ± 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.

Copy number alterations in EGFR/KRAS mutation groups

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). 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). 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).

Figure 3

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 log2ratio ± 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.

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 and Table 3). 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.

Table 3

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 Less 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 and S2).

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).

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

Patterns of tumor ploidy and allelic imbalances between mutational groups were evaluated by GAP [30] 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). 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 and Figure S2).

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

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). 96 genes showed consistent differential expression in ≥4 cohorts, while only 21 genes were differentially expressed across all five cohorts (Tables 4 and S3). We validated the 96 identified genes in four independent adenocarcinoma cohorts analyzed by different microarray platforms (n=290 tumors, Table 2). 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).

Table 4

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 and Table S3). Classification of the four independent cohorts showed an overall accuracy of 40-90% in classification across a range of classification cut-offs (Figure S3A). Sensitivity was highest in classification of EGFR-mutated tumors across the different cohorts (80-100%), followed by KRAS-mutated tumors (Figure 4A). However, specificities of the 96-gene classifier were lower (60-90%) for the EGFR and KRAS-mutated groups (Figure 4A). In contrast, for EGFRwt/KRASwt tumors, sensitivity was poor (10-60%) but specificity higher (80-100%).

Figure 4

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.

To further analyze the transcriptional patterns between the three mutation groups we performed unsupervised hierarchical clustering of four Affymetrix discovery cohorts (Table 2). Each cohort was individually clustered using a) a signature of genes overexpressed in EGFRwt/KRASwt tumors [26], b) a KRAS dependency gene signature [33], 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). 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 [33] (Figure S4).

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

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). Overall, EGFR-mutated tumors displayed more CNAs, more amplifications, and higher genomic complexity than non-EGFR-mutated tumors consistent with previous reports [11,15,18]. 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 [34]. Moreover, FGF-FGFR pathway activation has been suggested to be one mediator of resistance to EGFR inhibitors, together with, e.g., MET amplification (see [35] for review and [36]). 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 [9,11-13,17,18], while others such as 5q34-q35.3, appear novel (see Table S4 for literature comparison of 34 previously reported regions from five independent studies [9,11-13,18]). 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). 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 [18], 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 [37]. 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 [21], and also recent reports from breast cancer [38,39], 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 [8,10,12,14,16,26,33,40]. However, the overlap between these public signatures is very poor when directly compared (Figures S3B and C). 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 [8,12,14,16]. 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 [8,14,16,26,33,40] 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 [8,12,14,16,26] (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) [10,41], 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) [42-45].

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. [10] 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 [16,19,20]. Although the bronchioid molecular subtype [19] has been strongly associated with EGFR-mutated tumors, this subtype also includes notable fractions of KRAS-mutated and EGFRwt/KRASwt tumors (see, e.g., [19,20]). Moreover, ~30% or more of EGFR-mutated have been classified as non-bronchioid (magnoid or squamoid) in discovery cohorts in previous studies [19,20,46]. 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). These results also suggest that EGFR-mutated tumors could be divided into additional subgroups, which we have recently demonstrated [46]. 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 [26], 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 [17,20,47,48]. 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 [17]. 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 [46]. 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.

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 -log10(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 [28] 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 [33,40].

(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 [26]. (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. [33]. 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.