Development of somatic mutation signatures for risk stratification and prognosis in lung and colorectal adenocarcinomas

Background Prognostic signatures are vital to precision medicine. However, development of somatic mutation prognostic signatures for cancers remains a challenge. In this study we developed a novel method for discovering somatic mutation based prognostic signatures. Results Somatic mutation and clinical data for lung adenocarcinoma (LUAD) and colorectal adenocarcinoma (COAD) from The Cancer Genome Atlas (TCGA) were randomly divided into training (n = 328 for LUAD and 286 for COAD) and validation (n = 167 for LUAD and 141 for COAD) datasets. A novel method of using the log2 ratio of the tumor mutation frequency to the paired normal mutation frequency is computed for each patient and missense mutation. The missense mutation ratios were mean aggregated into gene-level somatic mutation profiles. The somatic mutations were assessed using univariate Cox analysis on the LUAD and COAD training sets separately. Stepwise multivariate Cox analysis resulted in a final gene prognostic signature for LUAD and COAD. Performance was compared to gene prognostic signatures generated using the same pipeline but with different somatic mutation profile representations based on tumor mutation frequency, binary calls, and gene-gene network normalization. Signature high-risk LUAD and COAD cases had worse overall survival compared to the signature low-risk cases in the validation set (log-rank test p-value = 0.0101 for LUAD and 0.0314 for COAD) using mutation tumor frequency ratio (MFR) profiles, while all other methods, including gene-gene network normalization, have statistically insignificant stratification (log-rank test p-value ≥0.05). Most of the genes in the final gene signatures using MFR profiles are cancer-related based on network and literature analysis. Conclusions We demonstrated the robustness of MFR profiles and its potential to be a powerful prognostic tool in cancer. The results are robust according to validation testing and the selected genes are biologically relevant.


Background
Lung and colon cancer are the leading cause of death over all cancers in the United States in 2017, with 155,870 and 50,260 deaths, respectively [1]. Prognostic signatures and risk stratification are vital to clinical decision making of treatment options in cancer precision medicine. As patient prognosis remains poor [2], researchers are seeking to develop improved prognostic signatures using molecular information, such as incorporating long non-coding RNA expression [3,4].
However, incorporating somatic mutation profiles into prognostic signatures has remained a challenge and is often overlooked due to the sparse and binary nature of somatic mutation data [5]. The sparsity of the data arises from the observation that the vast majority of mutated genes are not shared among patients [6]. Save for a few frequently mutated driver genes, most somatically mutated genes are likely to be composed of only passenger mutations that do not provide growth advantage [7].
To investigate the prognostic value of somatic mutations, studies have chosen to tackle the challenge by confronting the sparsity problem. Le Morvan et al. [8] uses gene-gene networks as prior knowledge to de-sparsify the data. A patient's binary somatic mutation profile is transformed by removing non-essential mutations and adding proxy mutations based on gene-gene network topology to normalize tumor mutational burden within a sample of patients. However, gene-gene networks vary from tissue to tissue and a single set of canonical gene-gene networks as prior knowledge may omit or overemphasize some interactions [9]. To address this issue, other studies have elected to use cancer-specific co-expression networks based on RNA expression data [10] or canonical pathways [11].
In this study, we confront the challenge of the binary nature of somatic mutation data rather than the sparsity problem. We propose the usage of the quantitative mutation frequency ratio of tumor vs. normal tissue from whole exome sequencing in building somatic mutation profiles. Using somatic mutation data for lung adenocarcinoma (LUAD) and colorectal adenocarcinoma (COAD) from The Cancer Genome Atlas (TCGA) [12,13], we evaluate the risk stratification and prognostic performance of somatic mutation signatures generated by using two types of continuous somatic mutation profiles: mutation frequency ratio (MFR) profiles and tumor mutation frequency (TMF) profiles. We compare to two existing types of binary mutation profiles, raw binary mutation (BM) profiles and gene-gene network normalized profiles provided by NetNorM [8]. We show that the somatic mutation signatures generated by MFR profiles consistently provides statistically significant risk stratification while the other types of profiles do not.

Identification of prognostic somatically mutated genes
To identify and evaluate prognostic somatically mutated genes using different types of somatic mutation profiles, we used a pipeline (Fig. 1) adapted from Shukla et al.'s RNA-seq pipeline [3]. Clinical and controlled somatic mutation data for LUAD and COAD was gathered from TCGA [12,13]. The data ( Table 1) was partitioned randomly into training (n = 328 for LUAD and n = 286 for COAD) and validation (n = 167 for LUAD and n = 141 COAD) datasets and somatic mutation profiles generated.
Four different types of somatic mutation profiles were considered: MFR, TMF, BM, and NetNorM profiles. The somatic mutation profile of a single patient is a vector with an element for every gene. The BM profile of a patient consists of a sparse binary vector where an element denotes if a gene is somatically mutated or not. The NetNorM profile was generated from the BM profile by normalizing the number of mutated genes via the removal or addition of somatically mutated genes [8]. While the NetNorM profile remains binary in nature, its process mitigates the sparsity problem of somatic mutation data by incorporating gene-gene network prior knowledge.
Additionally, we propose the usage of MFR and TMF profiles, which to the best of our knowledge, has not be considered previously in the literature to confront the difficulties of working with sparse binary data. TMF profiles incorporate the tumor data on the number of reads supporting the mutation vs. the reference genome. The MFR takes it a step further and considers the mutation frequency ratio of the tumor sample vs. the paired normal tissue sample. Both TMF and MFR profiles use continuous rather than binary values for somatic mutation profile representation.
Individually for each type of somatic mutation profile and tumor type, somatic mutation based prognostic signatures are generated using the pipeline outlined in Fig. 1. Univariate Cox proportional hazards regression is first performed on the training dataset to short list prospective genes with a FDR cutoff of 0.05. The prospective genes are then subjected to bidirectional stepwise multivariate Cox proportional hazards regression model selection to the determine the final prognostic signature ( Table 2 and Table 3). We verified that all of the final prognostic signatures do not violate the proportional hazards assumption using the Schoenfeld Residual Test.

Comparison of risk stratification
Kaplan-Meier (KM) survival curves are used to assess and compare the different types of somatic mutation profiles in both the training and validation datasets. Using the final Cox model for risk scoring, the high-risk threshold for stratification in both the training and validation datasets was chosen to be the 75th percentile of the risk scores in the training dataset.
We observed that all somatic mutation profile types achieve significant risk stratification on the training dataset (log rank test p-value ≈ 0) for both LUAD and COAD (Fig. 2, Fig. 3). For both LUAD and COAD, however, only the stratification generated by MFR profiles is statistically significant in the validation datasets (log rank test p-value = 0.0101 for LUAD, 0.0314 for COAD) (Fig. 2a, Fig 3a), while all    AOC1  FALSE  TRUE  TRUE  FALSE   TLR9  FALSE  TRUE  FALSE  FALSE   MOSPD2  FALSE  TRUE  TRUE  TRUE   EPHA2  FALSE  TRUE  TRUE  TRUE   ZNF880  FALSE  TRUE  FALSE  FALSE   TAS2R39  FALSE  TRUE  FALSE  FALSE   DNTTIP1  FALSE  TRUE  FALSE  FALSE   HHAT  FALSE  TRUE  TRUE  FALSE   ALOXE3  FALSE  TRUE  TRUE  FALSE   PRMT5  FALSE  TRUE  TRUE  FALSE   FAM83B  FALSE  TRUE  FALSE  FALSE   BEST4  FALSE  TRUE  FALSE  FALSE   BCAS3  FALSE  TRUE  FALSE  FALSE   MAP3K1  FALSE  TRUE  FALSE  FALSE   GPR52  FALSE  TRUE  FALSE  FALSE   DNAJC10  FALSE  TRUE  FALSE  FALSE   ADGRG7  FALSE  TRUE  FALSE  FALSE   CDRT15  FALSE  TRUE  FALSE  FALSE   MOCS3  FALSE  TRUE  FALSE  FALSE   C5  FALSE  TRUE  FALSE  FALSE   CNTN1  FALSE  TRUE  FALSE  FALSE   CLCN2  FALSE  TRUE  FALSE  FALSE   CBLB  FALSE  TRUE  TRUE  TRUE   MSH3  FALSE  TRUE  FALSE  FALSE other profiles, including NetNorM, are not statistically significant (Fig. 2b, c and d, Fig. 3b, c and d). Furthermore, the final prognostic signatures generated by each type of somatic mutation profile only minimally overlap for both LUAD and COAD cases (Fig. 4).
The results suggest that the MFR profile's prognostic signature is more robust, while the other types of profiles are subject to harsh overfitting that is typical in contexts with a larger number of covariates than samples. This is consistent with the observation that NetNorM profiles typically do not perform statistically different from binary profiles [8]. De-sparsifying somatic mutation data using gene-gene network prior information does not necessarily lead to improved prognostic and risk stratification perfor mance.

Somatic mutation gene signatures
A PubMed search of the individual genes and a network analysis of the full signatures using Ingenuity Pathway Analysis (QIAGEN Inc., https://www.qiagenbioinformatics.com/   HIF1AN  FALSE  TRUE  FALSE  FALSE   IGHA1  TRUE  FALSE  FALSE  FALSE   IQCH  TRUE  FALSE  FALSE  FALSE   KANSL3  TRUE  TRUE  FALSE  FALSE   KRT73  FALSE  FALSE  FALSE  TRUE   MARCH11  TRUE  FALSE  TRUE  FALSE   MEOX1  TRUE  FALSE  FALSE  FALSE   METTL21C  TRUE  FALSE  TRUE  TRUE   MICA  TRUE  TRUE  TRUE  FALSE   NAV1  FALSE  FALSE  TRUE  FALSE   NKD1  TRUE  TRUE  TRUE  FALSE   NTSR1  FALSE  TRUE  FALSE  FALSE   OGFR  FALSE  FALSE  FALSE  TRUE   OR10A7  FALSE  TRUE  FALSE  FALSE   OR10H2  FALSE  FALSE  FALSE  TRUE   OR11H1  FALSE  FALSE  FALSE  TRUE   OR13C8  FALSE  TRUE  FALSE  FALSE   OR1D5  FALSE  FALSE  TRUE  TRUE   PDHB  TRUE  FALSE  FALSE  FALSE products/ingenuity-pathway-analysis/, accessed: Feb. 14, 2018) was performed to assess the biological relevancy of the final prognostic gene signatures generated by MFR profiles. A network containing 16 of the 20 genes in the LUAD prognostic signature (Table 4) was found (Fig. 5). The network is associated with cell death and survival, and cellular movement. All genes in the prognostic signature are positively associated with risk (denoted in red in Fig. 5). SDHA is the gene with the largest coefficient in the risk model (hazard ratio (HR) = 1.877). SDHA is a tumor suppressor and is implicated in paraganglioma and gastrointestinal stromal tumors [14]. While association of SDHA copy number variation to prognosis was found in lung squamous cell carcinoma [15], we have found no literature exploring the connection of SDHA to lung adenocarcinoma. Four additional genes in the LUAD signature also have known associations with lung cancer. PFKM has mutations associated with survival outcomes in lung squamous cell carcinoma [16]. MADD promotes survival of LUAD cells and is a potential therapeutic target [17]. SERPINI2 is tumor suppressor gene and is associated with squamous cell lung cancer [18]. Finally, it has been found that certain MMP8 mutations are correlated with risk of developing lung cancer [19].
Eight of the remaining genes in the LUAD signature are associated with other cancer types and their connection to LUAD is yet uncharacterized. ABCB6 [20,21], ZNF768 [22], and the TP53-mediated tumor suppressor gene NDN [23] are all associated with colorectal cancers. MSANTD3 is an oncogene in salivary gland acinic cell carcinoma [24]. FGD3 is implicated in breast cancer [25] and ARHGAP4 in ovarian tumors [26]. It has been observed that increased expression of HSD17B4 is correlated with poor prognosis in prostate cancer List of somatically mutated genes selected by the pipeline for COAD using each type of somatic mutation profiles [27]. Lastly, correlation of HEATR1 with shorter overall survival has been shown in pancreatic ductal adenocarcinoma [28]. For the COAD prognostic signature (Table 5), we found that 30 of the 32 genes were involved in two different networks. The first network contains 16 of the 32 genes in the COAD prognostic signature (Fig. 6) and is associated with embryonic, organismal, and tissue development. The second network contains 14 of the 32 genes in the COAD prognostic signature (Fig. 7) and is associated with cancer and organismal injury and abnormalities. Unlike the LUAD signature where all genes were positively associated with increased risk, mutations in seven of the genes are associated with reduced risk (USP50, UBTD1, ZNF83, FBX038, C11orf53, IQCH, and CHI3L1) and are denoted in green in Figs. 6 and 7.
Ten of the genes in the COAD signature are implicated in colorectal cancers (CRC). MICA has high cell-surface expression in cancers of the digestive system and have been found to be correlated with increased survival [29]. Copy number variation of RERG is correlated with CRC risk [30]. in CRC and other cancers [31]. Lower expression of UBTD1 correlates with worse prognosis [32]. SER-PINB3 has a driving role in more aggressive cellular phenotypes of CRC [33]. DMKN has been previously proposed as a biomarker of early-stage CRC [34]. PDHB diminishes the oncogenic effects of miR-146b-5p on the growth and invasion of CRC [35]. C11orf53 is a potential gene involved in CRC etiology [36]. CHI3L1 promotes macrophage recruitment and angiogenesis in CRC [37]. Lastly, alterations of CDH24 contribute to tumorigenesis, as CDH24 is important to the maintenance of cell adhesion [38].
Another nine genes of the COAD signature have known associations with other types of cancers, but not with CRC yet. DNALI1 [39] and MEOX1 [40] are associated with breast cancer. In particular, MEOX1 is correlated with poor survival of breast cancer patients. MARCH11 has been used as a biomarker in a methylation panel for early cancer detection and prognosis prediction in non-small cell lung cancer [41]. ARHGAP15 is correlated with survival in early-stage pancreatic ductal adenocarcinoma [42]. IGHA1 is associated with gastric tumorigenesis [43]. CER1 is associated with glioma [44]. SDR9C7 promotes lymph node metastasis in esophageal squamous cell carcinoma [45]. PRKG2 is associated with acute mast cell leukemia [46]. Finally, ZNF133 is potential biomarker for osteosarcoma [47].

Discussion
Cancer genomic data is increasingly becoming a hot topic in precision cancer medicine research, including the identification of therapeutic targets, biomarker-based clinical trials, and the study of genomic determinants of therapy response [48]. The signatures found in the present retrospective study are promising and their potential clinical integration should be further investigated with a prospective study.
While the results are promising, there are limitations to this initial work. Demographic and clinical data were not incorporated into the prognostic models. Gene expression data is also available for TCGA LUAD and COAD datasets. Integration of all data types could potentially improve prognostic and risk stratification performance and provide further biological insights. Furthermore, all types of cancer in TCGA should be analyzed for a future pan-cancer study.
The present study was also done at the gene level. There is potential that specific mutations to a gene may have different prognostic effects. However, with the sample size of TCGA data, it is not feasible to observe statistically significant results due to the increased sparsity of somatic mutation data at the specific mutation level. Further data or methods to mitigate the increased sparsity is required for further study.
The present work demonstrated the robustness of prognostic signatures using MFR profiles within TCGA LUAD and COAD VarScan-based somatic mutation data [49] by the partitioning of the data into training and validation datasets. As a result, the experimental and analysis protocols are consistent. The robustness with respect to different somatic mutation calling software within TCGA should be conducted, as calls from MuSE [50], MuTect [51], and SomaticSniper [52] are provided in addition to VarScan. Furthermore, the methods robustness to data generated from different experimental protocols, such as by investigating data generated by different institutions and projects, should be studied in the future.

Conclusions
To improve clinical tools and biological understanding of LUAD and COAD, we demonstrated a methodology to generating robust prognostic somatic mutation-based gene signatures. We demonstrated the robustness of MFR profiles and its potential to be a powerful prognostic tool in cancer, unlike other alternative types of somatic mutation profiles, TMF, BM, and NetNorM, that did not achieve statistically significant risk stratification in validation datasets. The genes selected by the methodology using MFR profiles was shown to be biologically relevant and has potential for use in effective management LUAD and COAD.

Somatic mutation data and profiles
Controlled TCGA somatic mutation data (VarScan MAF files [49]) were downloaded from NCI's Genomic Data Commons (https://gdc.cancer.gov/, accessed: Feb. 14, 2018) for LUAD and COAD (Project ID 17109, A Pan-Cancer Analysis of Somatic Mutation Profiles for Tumor Immunogenicity and Prognosis). The data were filtered, keeping only missense mutations. The missense mutations were then aggregated into gene level mutation profiles. For BM profiles, the gene is flagged as mutated if it contains any missense mutation. The NetNorM normalization method was used as a representative of somatic mutation profiles using gene-gene network information [8]. NetNorM uses networks from Pathway Commons (http:// www.pathwaycommons.org), which feature an integrated network data of public pathway and interaction databases. The user-specified parameter for NetNorM is the target number of mutated genes k. This parameter is set to the median number of mutated genes in the training dataset, which is 193 and 151 for LUAD and COAD, respectively. NetNorM ranks genes based on their mutation status and network connectedness. A patient's somatic mutation profile is normalized by setting only the top k genes as being mutated. Since mutated genes are always ranked higher than non-mutated genes, patients with more than k mutated genes will have lower ranked mutated genes set to unmutated, while patients with less than k mutated genes will obtain artificial proxy mutated genes.

Mutation frequency ratio and tumor frequency profiles
For patient i, the MAF files contain the number of reads supporting the reference allele for mutation j, TRC ij and NRC ij for tumor and normal samples, respectively. Analogously, denote the number of reads supporting the alternate allele, TAC ij and NAC ij for tumor and normal samples, respectively. The tumor and normal sample mutation frequencies, TMF ij and NMF ij , are computed using Eqs. (1) and (2), respectively. The mutation frequency ratio MFR ij is then simply the ratio of the tumor to normal sample mutation frequencies. To generate a patient's gene level MFR and TMF profiles, the mutations are aggregated by gene using the mean ratio or frequency within that gene.
Signature generation and statistical analysis TCGA clinical data were downloaded from NCI's Genomic Data Commons (https://gdc.cancer.gov/, accessed: Feb. 14, 2018) for LUAD and COAD. These data were partitioned randomly into training (n = 328 for LUAD and n = 286 for COAD) and validation (n = 167 for LUAD and n = 141 COAD) datasets. Rarely mutated genes in somatic mutation profiles were omitted when less than 1% of patients in a sample have the mutation. MFR and TMF profiles, which are continuous valued, were log2 transformed. Univariate Cox proportional hazards regression was used to assess association with overall survival using R survival package (R v3. 4