Integrative analysis of lung adenocarcinoma across diverse ethnicities and exposures

Abstract

Lung adenocarcinomas (LUAD) are a pressing global health problem with enduring lethality and rapidly shifting epidemiology. Proteogenomic studies integrating proteomics and post-translational modifications with genomics can identify clinical strata and oncogenic mechanisms, but have been underpowered to examine effects of ethnicity, smoking and environmental exposures, or sex on this heterogeneous disease. This comprehensive proteogenomic analysis of LUAD tumors and matched normal adjacent tissues from 406 patients across diverse geographic and demographic backgrounds explores the impact of understudied driver mutations, prognostic role of chromosomal instability, patterns of immune signaling, differential and sex-specific effects of endogenous mutagens and environmental carcinogens, and pathobiology of early-stage tumors with "late-like" characteristics. Candidate protein biomarkers are proposed for unstable tumors with highly fragmented genomes and for carcinogen exposures, and a LUAD subtype-specific atlas of therapeutic vulnerabilities is presented. These observations and the associated data resource advance the objective of precision management strategies for this devastating disease.

Keywords: CPTAC; ICPC; PTMs; biomarkers; genomics; lung adenocarcinoma; mass spectrometry; never-smokers; proteogenomics; smokers.

Figure 1.. Proteogenomic landscape of the CPTAC and ICPC cohorts

A. Left: Summary of included samples and -omes. Sample numbers are stratified into tumor and normal adjacent tissue (NAT). The -omes available are listed to the right of each cohort (ICPC, top; CPTAC, bottom). * denotes that glycoproteome data is available only for the CPTAC.B cohort. Right: Distribution of selected annotations across 406 patients profiled in this study. All annotations are self-reported except for EGFR mutation, which was determined using whole exome mutation calls. “Unknown” denotes that the data were not reported, or in the case of EGFR mutation, that mutation calls were not available for those samples. B. Oncoplot of mutations of interest across the cohort of 406 patients. TMB per MB: Tumor Mutation Burden per Megabase. Selected mutation calls are grouped together in the following categories: RAS Family (HRAS, KRAS, NRAS); RAF Family (ARAF, BRAF, RAF1); TKI Fusion (ALK, RET, ROS1). C. Cell-of-origin analysis across the 402 tumor samples for which RNA-seq passed QC. The 402 samples are grouped based on mutation status. An additional 14 squamous samples collected and analyzed under the auspices of this study were included exclusively in this analysis. D. Multi-omic NMF clustering across 383 tumor samples with available RNA-seq, copy number variation, proteomic, and phosphoproteomic data. E. Clinical annotations significantly enriched in each multi-omic cluster. Color represents nominal p-value from Fisher's exact test on the NMF consensus core samples. Annotations shown have nominal p-value < 0.01 in at least one of the clusters. F. Gene Set Enrichment Analysis (GSEA) on the four multi-omic clusters. Heatmap shows the Normalized Enrichment Score (NES), which is derived from the feature weights for each cluster. * indicates significant enrichment at adj. p<0.01. See also Figure S1 and Tables S1 and S2.

Figure 2.. Impact of mutation and PTMs on downstream signaling and structure

A. GSEA result of protein-level differential expression analysis comparing RBM10 truncated mutation tumors versus RBM10 WT in the context of EGFR mutation. B. Post-translational modification-Signature Enrichment Analysis (PTM-SEA) showed splicing kinase- and pathway-associated phosphorylation enrichment. C. Differentially expressed features between RBM10 mutant and WT tumors in the context of EGFR mutation at the levels of proteome and phosphoproteome. D. Differential features between ALK fusion and WT tumors at the level of proteome and phosphoproteome. E. Candidate pY sites and their abundance in a range of cell lines treated with driver-matched kinase inhibitors including ALK-fusion LUAD cell lines treated with the ALK inhibitor ceritinib. pY sites were downregulated upon kinase inhibition specifically in the ALK lines, and in a time dependent fashion (6 and 24 hours). F. Rank plots showing protein abundance (indicated by signed −log10 p-value) in candidate pY vs. unmodified peptide pull-downs. The red dots show known ALK interactors in the BioGRID database. Key interactors are indicated in the panel. G. Bubble plot representing CLUMPS-PTM results. Phosphoproteome shows clustering results for both ICPC and CPTAC cohorts; acetylproteome and ubiquitylproteome results are shown for the CPTAC cohort. Red circle: Significant results: permutation test, adj. p<0.1. H. HSPB1 phosphorylation cluster on 3D crystal structure (cyan, PDB ID: 6DV5-A, left). Violin plots showing relative protein abundances of AKT1 between EGFR mutant and wild type (right). In boxplots inset in violin plots, the centerline represents median; box bounds represent 25th and 75th percentiles; whiskers represent minimum and maximum; and filled points represent individual data points. I. MDH1 acetylation cluster on 3D crystal structure (cyan, PDB ID: 4WLE-A). Highly differentially expressed acetylation sites are denoted by bright pink. J. PSMA5 ubiquitylation cluster on 3D crystal structure (cyan, PDB ID: 5GJQ-F, left). Violin plots showing protein abundances of PSMA5 across the proteomic NMF clusters (right). In boxplots inset in violin plots, the centerline represents median; box bounds represent 25th and 75th percentiles; whiskers represent minimum and maximum; and filled points represent individual data points. See also Figure S2 and Table S3.

Figure 3.. Genomic instability in LUAD and its impact

A. Distribution of CNV segment lengths for each sample. Each column on the x-axis represents one sample and the colors represent the density of the segment length for that sample. Samples are sorted based on 75th percentile of segment lengths. Samples on the left with a high density of large segments represent “Full-stable” and samples on the right with a high density of short segments represent “Unstable”. B. Kaplan-Meier (KM) survival curves showing the overall survival of patients, stratified by SegLen-Q3. P-value obtained using log-rank test. C. Comparison of recurrently amplified regions in unstable (top) and fully stable (bottom) tumors. D. Classification of tumors according to their breakpoint intensity (x-axis) and clustering score (y-axis) to derive BIC groups. Tumors are categorized into three groups: tumors with highly clustered breakpoints (Contiguous), tumors with low clustering and intensity of breakpoints (Fragmented), and tumors with low clustering and high intensity of breakpoints (Intense). E. KM survival curves showing the overall survival of patients, stratified by BIC. Tumors with high clustering show exceptionally good survival independent of their breakpoint intensity. In groups with low clustering, high intensity of breakpoints results in worse outcomes. P-value obtained using log-rank test. F. The proportion of tumors with NKX2-1 amplification (red) and null expression (blue). Tumors with amplified NKX2-1 show significant enrichment in the BIC-intense group. G. Proliferative activity of tumors stratified by NKX2-1 amplification (left) and expression (right) status. Proliferative activity is measured by ssGSEA analysis of RNA expression for the HALLMARK_G2M_CHECKPOINT pathway. P values are based on the Wilcoxon Rank Sum test. In boxplots, the centerline represents median; box bounds represent 25th and 75th percentiles; whiskers represent minimum and maximum; and filled points represent individual data points. H. Transcriptional state of tumors stratified by BIC groups. The BIC-intense group shows a reduced contribution of the AT2 phenotype. I. Log2 fold changes of genes from differential expression analysis comparing BIC-intense to other tumors at both RNA (x-axis) and protein (y-axis) levels. J. IGF2BP3 IHC and H&E staining shows heterogeneous protein expression of IGF2BP3 in case-3 (partial stable LUAD), uniformly strong expression in case-1 (unstable), and complete abrogation of IGF2BP3 expression in case-4 (full stable). (IHC; H&E). Scale bar = 50μm. K. Independent validation in 3 patient TMAs reveals higher IGF2BP3 expression in metastatic LUAD compared to primary LUAD and primary LSCC. Also noteworthy is the absence of IGF2BP3 in benign lung parenchyma cores. In boxplots, the centerline represents median; box bounds represent 25th and 75th percentiles; whiskers represent minimum and maximum; and filled points represent individual data points. L. Expression of IGF2BP3 (RNA) in LUAD patients from ICGC/PCAWG, stratified by instability. In boxplots, the centerline represents median; box bounds represent 25th and 75th percentiles; whiskers represent minimum and maximum; and filled points represent individual data points. M. KM survival curves showing the overall survival of non-squamous patients treated with the anti-PDL1 antibody MPDL3280A, stratified by IGF2BP3 RNA expression in the OAK trial. P-value obtained using log-rank test. N. PD-L1 (CD274) protein expression in the CPTAC cohort stratified by BIC status. In boxplots, the centerline represents median; box bounds represent 25th and 75th percentiles; whiskers represent minimum and maximum; and filled points represent individual data points. O. Correlation of IGF2BP3 protein expression and average IGF2BP3 promoter methylation. See also Figure S3 and Table S3.

Figure 4.. Immune landscape of LUAD

A. Cell type fractions estimated via multi-omics-based deconvolution analysis (top), protein (middle), and RNA expression (bottom) of cell type-specific markers. The annotation tracks show immune subtypes derived from proteogenomic data, immune score derived via ESTIMATE from RNA-seq data, tumor purity derived via TSNet from RNA-seq data, tumor purity estimated from CNV data, and various clinical and mutation parameters. B. Bubble plot showing summary statistics of association analyses between immune subtypes and biological pathways (from Hallmark, KEGG and Reactome). Bubble size corresponds to adjusted p-value (Wilcoxon test, −log10 scale), while bubble color corresponds to the difference between the averaged pathway score for samples in a particular immune subtype and that of samples allocated to other immune subtypes. C. Bubble-plot showing the association analysis between immune subtypes and kinase activity score for tumor samples. The color of the bubble corresponds to the average kinase activity score in a particular immune subtype, while the size of the bubble corresponds to the adjusted p-value (Wilcoxon test, −log2 transformed). Associations significant at 10% FDR are denoted with an outer black circle. D. Focal adhesion-related kinase cluster associated with immune-related TF cluster in Hot tumors. The bar plots show experimental validation via the LINCS L1000 database. The leftmost blue bar chart shows the enrichment of Innate Immune System Reactome pathway in the set of genes upregulated after kinase CRISPR-Cas knockouts from the LINCS L1000 dataset (p-values from Fisher’s exact test). The centered bar chart shows the enrichment of STAT4 targets in the set of up-regulated genes after L1000 kinase knockouts. The rightmost chart displays the number of overlapping genes between STAT4 targets and the set of up-regulated genes after L1000 kinase knockouts. The red vertical line in the bar plots denotes significance (p = 0.05). Cell lines are in parentheses and their primary disease associations are: A549: Lung Cancer, AGS: Gastric Cancer, YAPC: Pancreatic Cancer, BICR6: Head and Neck Cancer, A375: Skin Cancer, ES2: Ovarian Cancer, HT29: Colon/Colorectal Cancer, MCF7: Breast Cancer, U251MG: Brain Cancer. E. Cell-type specific pathway analysis to characterize pathway activation in tumor cells of Hot and Cold tumors based on RNA (left) and proteome data (right). The heatmap shows the average pathway score for each cell type (i.e., tumor cells of Hot tumors, tumor cells of Cold tumors and immune/stromal cells. F. Selected glycosites stratified by immune subtype. P-values from the Wilcoxon test are reported. Asterisks represent p-values<0.001. In boxplots, the centerline represents median; box bounds represent 25th and 75th percentiles; whiskers represent minimum and maximum; and filled points represent outliers. G. Neutrophil infiltration stratified by immune subtype. P-values from the Wilcoxon test are reported. Asterisks represent p-values<0.001. In boxplots, the centerline represents median; box bounds represent 25th and 75th percentiles; and whiskers represent minimum and maximum. H. KM curves of progression-free survival stratified by immune subtypes. P-value from a log-rank test is reported. I. Association between STK11 mutation status and fraction of different cell types in tumor samples. P-values from the Wilcoxon-rank-sum test are reported. In boxplots, the centerline represents median; box bounds represent 25th and 75th percentiles; and whiskers represent minimum and maximum. J. CD47 stratified by STK11 mutation status and NMF proteomic subtype. P-values from the Wilcoxon test are reported. Asterisks represent p-values<0.001. In boxplots, the centerline represents median; box bounds represent 25th and 75th percentiles; whiskers represent minimum and maximum; and filled points represent outliers. K. Scatter plot showing molecular aberrations between C1/C3 STK11-mutant vs C2 STK11-mutant samples. See also Figure S4 and Table S4.

Figure 5.. Environmental Carcinogens Synergistically Induced LUAD Development.

A. Ten mutation-based subclusters (Figure S5A) are condensed into four major mutagen/carcinogen groups based primarily on the relative contribution of nitro-PAHs (dark blue), PAHs (light blue), nitrosamines (green), APOBEC (yellow) and methylcytosine (gray) signatures. A mixed cluster is indicated in white. B. Summary of association of environmental carcinogens with clinical features and gene mutations. Tumor samples are categorized into four groups based on the relative contribution of exogenous carcinogen signatures: PAHs/nitro-PAHs high, PAHs/nitro-PAHs low, Nitrosamines high, and Nitrosamines low. Clinical features with significantly different enrichment between these groups are annotated and marked with an asterisk. (Fisher’s exact test, adj. p<0.05). C. Summary plot showing different endogenous mutagen and exogenous carcinogen signatures (x-axis) in tumors and their enriched EGFR and KRAS mutations segregated by smoking status (y-axis). D. KM plot showing significant differences in relapse-free survival (log-rank test, p<0.001) between the four carcinogen groups defined in Figure 5A. E. Forest plot of the impact of carcinogen groups on relapse-free survival using the Hazard ratio (HR) with 95% confidence interval (CI). F. KM plot showing significantly worse relapse-free survival (log-rank test, p=0.018) in patients in the high compared to low nitrosamine signature group. G. Canonical pathway activation (activating z-score >1, adj. p<0.05) based on Ingenuity Pathway Analysis of proteins differentially expressed (Fisher's exact test, adj. p<0.05) in NATs from patients in the PAHs/nitro-PAHs and nitrosamine carcinogen signature groups. H. Signaling pathway enrichments in high smoke-exposure score (HSS) NATs based on a protein-protein interaction network derived from proteins differentially expressed between HSS and low smoke-exposure score (LSS) NATs and partitioned into 19 strongly intra-connected modules (hypergeometric test, adj. p <0.05, STAR methods). I. Upregulated proteins and phosphosites associated with metabolic imbalance and chronic inflammation from NATs with HSS, mostly contributed by samples with PAHs/nitro-PAHs signatures (adj. p<0.05). The protein interaction networks and functional categories are represented as connecting lines according to the STRING database. See also Figure S5 and Table S5.

Figure 6.. Nitrosamines and Nitro-PAHs/PAHs carcinogens promote different cancer-relevant pathways

A. Canonical pathway activation (activating z-score >1, adj. p<0.05) based on Ingenuity Pathway Analysis of proteins differentially expressed (Fisher's exact test, adj. p<0.05) in tumors from patients in the PAHs/nitro-PAHs and nitrosamine carcinogen signature groups. B. Differential protein expression profiles demonstrate the impact of specific carcinogens on LUAD tumors. Comparison between tumors manifesting high or absence of nitrosamine signatures (upper panel) or between high, low, or absence of PAHs/nitro-PAHs signatures (lower panel; Table S6) reveal activation of organonitrogen/aromatic compound metabolic processes, immune response, tumor progression and metastasis (activating z-score >1; Fisher's exact test, adj. p<0.05). Specific enriched proteins and phosphosites (Fisher's exact test, adj. p<0.05) underpinning shared carcinogenic processes are color-coded for the associated carcinogen signature (nitrosamine, green; PAHs, light blue; nitro-PAHs, dark blue). Vertical tracks annotate smoking status and key driver mutations (EGFR and KRAS) for the carcinogen groups. Asterisks represent p-values (Wilcoxon rank sum test, * p< 0.05, ** p<0.01, ***p<0.001 and ****p<0.0001). C. Protein biomarker candidates differentially upregulated in tumors with HSS or LSS from self-reported smokers or never-smokers (Kruskal Wallis tests, adj. p<0.05) further filtered for the ability to distinguish between any two groups (Wilcoxon rank sum test, adj. p<0.05). The vertical tracks annotate proteins as kinases, glycosylation-regulated enzymes, oncogenes, tumor suppressors, or immune response markers (ONC, oncogene; TSG, tumor suppressor gene.) D. Metacore pathway analysis of differential proteomic and phosphoproteomic profiles of 1-NP- (5 uM) or BaP- (5 uM) treated A549 and PC9 cells. The differential protein numbers in the enriched pathways (Fisher's exact test, adj. p<0.05) are presented by the bubble plot. E. Selected proteins or phosphosites from differential expression profiles of 1-NP- and BaP-treated A549 and PC9 cells, with tracks showing enriched functions in metabolism, EGFR signaling, cell cycle and proliferation, metastasis or EMT from Metacore pathway analysis. Asterisks indicate proteins consistently regulated in nitro-PAHs and PAHs tumor signatures. F. Migration ability of lung cancer cells treated with carcinogens 1-NP and BaP. For panels F, G, and H, at least two independent experiments were performed, data are presented as mean ± SD, and p-values were determined by one-way ANOVA followed by Dunnett's multiple comparisons test. Asterisks indicate p< 0.05. G. Invasion ability of lung cancer cells treated with carcinogens 1-NP and BaP. Asterisks indicate p< 0.05. H. Cell proliferation assay for lung cancer cells treated with carcinogens 1-NP and BaP. Asterisks indicate p< 0.05. See also Figure S6 and Table S6.

Figure 7.. Proteomics Clustering Identifies Advanced Molecular Features Driving Late-like Early-stage LUAD

A. KM plots for relapse-free survival in tumors of stages I–III (left) and stage I alone (right) (log-rank test, p<0.001). B. Summary of activated canonical pathways (activating z-score >1; Fisher's exact test, adj. p<0.05) based on Ingenuity Pathway Analysis of proteins differentially expressed in tumors of all stages (I–IV, left), C1-C3 clusters including tumors of all stages (center), and C1-C3 clusters subsetted for stage I tumors only (right). C. Summary of activated canonical pathways (activating z-score >1; Fisher's exact test; asterisks indicate Benjamini-Hochberg adj. p<0.05) based on Ingenuity Pathway Analysis of proteins differentially expressed in C2 vs non-C2 tumors of all stages (I–IV, left) and stage I alone (right), separated by smoking status. D. Top 30 candidate protein biomarkers significantly enriched in “late-like” early-stage tumors in the proteomic C2 subtype (Fisher's exact test, adj. p<0.01). The log₂ and relative abundances in stage I-IV or stage I in the C2/non-C2 subtype, −log p-value, protein function and annotation of these candidate biomarkers are presented. E. Sex-specific differences in the incidence frequency (x-axis) and median values of mutagen/carcinogen exposure signatures (y-axis) in proteomic clusters C1–3 (p-values from Fisher’s exact test). F. Sex-dependent mutagen/carcinogen signature distributions within the proteomic C2 subtype. Lower tracks show clinical and mutational annotations. G. Summary of shared and sex-specific enrichment of activated pathways in the C2 proteomic subtype (left). EGFR and BCR signaling-related proteins as well as the immunoproteasome are highly activated in C2 females, while C2 male patients have upregulated KRAS pathways and standard proteasome proteins. Shared activating proteins in the C2 subtype are involved in immune response, cell death and survival, metabolism, and metastasis (right). See also Figure S7 and Tables S1 and S7.

Figure 8.. Candidate drug targets and biomarkers with potential translational utility

A. Plots depicting the top 15 potentially druggable targets identified in LUAD subtypes. B. Relative abundance of EGFR protein (top panels) and EGFR_K716 ubiquitylation sites (middle panels) in tumors and corresponding NATs in EGFR-WT and EGFR-mutant tumors (p-values derived from Wilcoxon rank-sum test) and EGFR shRNA dependency scores from DepMap (lower panels; lower scores indicate greater dependency). In boxplots, the centerline represents median; box bounds represent 25th and 75th percentiles; whiskers represent minimum and maximum; and filled points represent outliers. C. Same as in (B) except depicting CSF1R protein and CSF1R_K812 ubiquitination sites. Statistical analyses and boxplots are as described in (B). D. Boxplots comparing response to a HSP90AA1 inhibitor, tanepsimycin (PRISM) in STK11-mutated or -WT LUAD cell lines. P-value derived from t-test. In boxplots, the centerline represents median; box bounds represent 25th and 75th percentiles; whiskers represent minimum and maximum; and filled points represent outliers. E. Potentially druggable targets for NMF proteomic cluster C1 tumors as defined by candidate upregulation in tumors vs paired NATs using multiple proteomic and PTM data types and a CRISPR effect score below zero in cluster-matched LUAD cell lines. The top 15 overexpressed, targetable dependencies are labeled. Numbers in the bottom right corner indicate the total number of overexpressed dependencies and targets from each drug target tier. F. Same as in (E) except for NMF proteomic cluster C2 tumors. G. Same as in (E) except for NMF proteomic cluster C3 tumors. H. Summary of major translational findings from this study. See also Figure S8 and Table S8.