Identification of MicroRNA-Related Tumorigenesis Variants and Genes in the Cancer Genome Atlas (TCGA) Data

Article type: Research Article

Keywords: TCGA, microRNA, 3′UTR, cancer, single nucleotide variant, dbNSFP

Authors: Chang Li, Brian Wu, Han Han, Jeff Zhao, Yongsheng Bai, Xiaoming Liu ⚬

Affiliations: USF Genomics & College of Public Health, University of South Florida, Tampa, FL 33612, USA; lic@usf.edu; Huron High School, Ann Arbor, MI 48105, USA; bwus127@gmail.com (B.W.); johnhanhan298@gmail.com (H.H.); Marine Military Academy, Harlingen, TX 78550, USA; jeffzhao100@gmail.com; Next-Gen Intelligent Science Training, Ann Arbor, MI 48105, USA; Department of Biology, Eastern Michigan University, Ypsilanti, MI 48197, USA

License: © 2020 by the authors. CC BY 4.0 Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Article links: DOI: 10.3390/genes11090953  |  PubMed: 32824926  |  PMC: PMC7565843

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (1.4 MB)

1. Introduction

MicroRNAs (miRNA) are small (21–22bp) noncoding RNAs that play a post-transcriptional regulatory role through targeting the 3′ untranslated regions (UTRs) of mRNA(s). MiRNAs regulate mRNA(s) via base pairing through complementary sequences within mRNA molecules. Genetic mutations, especially Single Nucleotide Variants (SNVs) located on the miRNA target site (MTS), can disrupt the binding of targeting miRNAs. This disruption has been reported to be associated with numerous diseases and cancers [ref. 1].

Given their functional importance, our understanding of the role of MTS SNVs towards cancer is very limited, and many such SNVs have not even been identified. The Cancer Genome Atlas (TCGA) project provides us with a great opportunity to tackle this problem. TCGA includes sequencing data for 33 cancer types with over 20,000 primary cancer and normal tissue samples. More importantly, it has expression information for mRNAs and miRNAs in both normal and tumor tissues for various cancer types. The anti-correlated miRNA–mRNA pairs in the same tissue between normal and tumor samples could provide us with unprecedented opportunities to investigate tumorigenesis through miRNA regulation and identify functional MTS SNVs associated with a specific cancer type.

Thus, in this study, we utilized such anti-correlated miRNA–mRNA pairs identified in a recent study, [ref. 2] to cross-check with publicly available databases—namely dbMTS (database of miRNA target site SNVs) [ref. 3] and dbNSFP (database for nonsynonymous SNPs’ functional predictions) [ref. 4]—for functionally annotated SNVs. This process validated existing cancer related SNVs and identified and prioritized potential SNVs for future functional and experimental validation. Specifically, our study has identified 23 genes (24 SNVs) from dbMTS and 33 genes (65 SNVs) from dbNSFP that passed our filtering criteria. Cancer association annotation has been conducted on these genes, and our preliminary results reported that several of the identified genes are associated with tumorigenesis.

2. Materials and Methods

Anti-correlated expressions of 15 cancer types’ miRNA–mRNA pairs from a previous study [ref. 2] were extracted from TCGA (http://cancergenome.nih.gov), and a total of 65,535 miRNA–mRNA pairs were taken as candidates to be used in subsequent analyses. Because miRNA can repress the expression of its targets, we hypothesize that these anti-correlated miRNA–mRNA pairs are likely to be dysregulated in the observed cancer types, which highlight their potential functional importance.

2.1. Retrieval and Screening of miRNA Target Site SNVs Using dbMTS

The dbMTS (database of MTS SNVs) is a database containing only single nucleotide variants (SNVs) that could potentially affect the binding of miRNAs. The dbMTS uses computer simulations with multiple popular miRNA target prediction tools to predict the effect of any potential SNV on miRNA. It includes more than 200 annotations for each of its SNVs, to help users assess their functional importance. One piece of important information in dbMTS that we considered was the name of the miRNA–mRNA pair that was disrupted or created by introducing the SNV.

Using our candidate miRNA–mRNA pairs and information of SNVs from dbMTS, we identified all the candidate SNVs in the 3′UTR that could help us select our candidate miRNA–mRNA pairs. A customized Python script was written to retrieve SNVs that could affect binding between the candidate miRNA–mRNA pairs in dbMTS. A brief summary of our method of identifying 3′ UTR SNVs is shown in Figure 1. Specifically, gene symbols in our candidate list were first standardized and converted to Ensembl transcript IDs using Ensembl Biomart (Ensembl Genes 100), as dbMTS includes transcript level predictions. Second, the suffixes ‘-3p’ and ‘-5p’ in dbMTS were ignored and removed, since the miRNAs in our candidate pairs are precursor miRNAs instead of mature miRNAs. Third, using GWAS Catalog [ref. 5], GRASP [ref. 6], and ClinVar [ref. 7], we annotated known cancer-related SNVs using dbMTS. We found 376 variants having at least one annotated cancer-related phenotype. Because a transcript’s 3′ UTR can overlap with coding regions of other transcripts, we then used VEP/Ensembl [ref. 8] annotation to remove those SNVs that could potentially reside in coding regions. The VEP/Ensembl annotations for each SNV were retrieved from dbMTS. This process resulted in 24 variants that were exclusively located in the 3′ UTRs of all possible transcripts. Additionally, to investigate if any SNVs identified from dbMTS (a germline SNV database) have been reported as somatic, we checked our identified SNVs that could potentially affect miRNA binding in SomamiR DB 2 [ref. 9]. SomamiR DB is a database exclusively designed to study somatic mutations altering miRNA–mRNA interactions. We retrieved somatic SNVs that were located in experimentally identified MTS using CLASH (Cross-linking, ligation, and sequencing of hybrids) or CLIP-seq (Cross-linking immunoprecipitation sequencing) methods.

2.2. Retrieval and Screening of Nonsynonymous SNVs

The dbNSFP database includes only nonsynonymous SNVs (SNVs that change the final protein product) in the human genome. It includes more than 300 unique functional annotations for each SNV.

We retrieved all SNVs in dbNSFP that meet the following criteria, as shown in Figure 2:

• SNVs located in genes that are included in our candidate gene list;
• SNVs that include one of the following strings reported for ClinVar phenotype in dbNSFP to ensure their functional relationship with cancer: (‘carcinoma’, ‘sarcoma’, ‘leukemia’, ‘tumor’, ‘cancer’, ‘lymphoma’, ‘myeloma’);
• FATHMM was used to filter for damaging variants by removing those SNVs that are predicted as ‘Neutral’;
• SNV allele frequency should not be high (minor allele frequency <0.05 in both 1000Gp3_AF and UK10K_AF) because pathogenic SNVs are not likely to be prevalent in the general population;
• SNVs that are pathogenic or likely pathogenic as reported by ClinVar clinical significance. These steps were achieved using a custom Python script.

2.3. Functional and Pathway Analyses

The Catalogue Of Somatic Mutations In Cancer (COSMIC) [ref. 10] database was used to identify the impact of somatic mutations in cancer. To identify somatic mutations in our findings, we downloaded COSMIC mutation data from its targeted and genome wide screens. Our candidate SNVs were checked for their somatic statuses using the variants’ HGVS-nomenclature as keys, and examined for gene–cancer association using COSMIC. The Kyoto Encyclopedia of Genes and Genomes (KEGG) [ref. 11] is a database that links gene catalogs to higher-level systemic functions of the cell, the organism, and the ecosystem. We used KEGG exclusively for gene level analysis. Diseases that were found to be associated with our 23 genes were then further analyzed through public websites to find any possible cancer associations. We also employed COSMIC and KEGG for function and pathway validation and conducted similar searches for cancer association of the 33 annotated genes obtained from dbNSFP. These steps of finding functional interpretations of our candidate SNV/gene lists were mainly achieved through manually searching the above-mentioned publicly available databases.

2.4. Survival Analyses

Survival analyses were performed using two web-based tools, GEPIA [ref. 12] and OncoLnc [ref. 13]. Both tools used multivariate Cox regression to campare survival outcomes between patient groups. Kaplan–Meier survival curves were used to visualize the survival differences between patient groups over time. One difference between these two tools is the way they pre-process the TCGA data. To check if this difference affects the outcome, we adopted both tools in the analysis. Our candidate gene list was submitted to these two web-based tools, and corresponding outputs were downloaded and compared. One thing to note regarding these candidate genes is that they would require further experimental studies to mechanically associate their functions to tumorigensis.

3. Results and Discussion

3.1. Twenty four Non-Coding SNVs were Identified at Potential miRNA Target Sites

As illustrated in Figure 1, we obtained 870,254 SNVs that could potentially affect our candidate miRNA–mRNA pairs after checking with dbMTS. Next, we identified 376 variants with at least one annotated phenotype (GWAS Catalog, GRASP, and/or ClinVar, Table S1). We identified 24 SNVs exclusively located in the 3′UTR of all available transcripts. These SNVs were predicted by TargetScan to cause either gain of miRNA regulation or substitution of existing miRNA regulation (Table S2). These candidate SNVs would be overlooked by most of the traditional in-silico functional prioritization tools, because these MTS SNVs are expected to show only moderate effect size to gene expression. This indicates the uniqueness of our approach in identifying cancer type-specific functional SNVs for further validation. These 24 candidate SNVs were in 23 different genes (Table 1), and all of these genes have 2 or more potential cancer associations according to our observed anti-correlated miRNA-gene pairs. We expect that genes involved in multiple cancer types (e.g., MAFB and SLC30A7) are especially useful in detecting common tumorigenesis pathways for various types of cancer under the influence of similar set of miRNAs. The candidate list of SVNs and genes could account for the cancer association as a result of the anti-correlation between the miRNA–mRNA pairs in the TCGA project, but further experimental and/or clinical validation is required to elucidate their contributions to tumorigenesis.

3.2. Two Hundred and Thirty Three Somatic SNVs were Identified in Experimentally Identified miRNA Target Sites

Among the 870,254 SNVs retrieved from dbMTS, we found 233 somatic SNVs in SomamiR DB. Eleven of these SNVs located in MTS were identified by CLASH based experiments (Table 2). These SNVs may be responsible for our observed anti-correlated miRNA–mRNA pairs. For example, gene FLAD1 has been reported to be overexpressed in breast cancer [ref. 14], and let-7b dysregulation has been reported to be associated with breast cancer [ref. 15]. While no previous study has linked let-7b and FLAD1 together with breast cancer, our analysis shows the possibility of a novel mechanism of tumorigenesis for breast cancer. Additionally, a recent study has linked CDH1 and the miR-30 family with pancreatic cancer [ref. 16], which is observed in our study (Table 2). This evidence highlights the ability of our approach to identify informative tumorigenesis SNVs and genes. Additional SNVs located in MTS identified by CLIP-seq-based experiments are provided in Table S3.

3.3. Sixty Five Nonsynonymous SNVs were Identified in the Candidate Gene List

We identified 28,764 variants from the dbNSFP database for all genes from the TCGA analysis results to conduct our filtering process. We utilized three filters in order to identify SNVs that could be associated with cancer. The first filter used was “damaging in the fathmm-XF_coding prediction”; this step kept 12,617 variants. Next, we filtered the variants using The 1000 Genomes Project Phase 3 data by requiring SNVs’ minor allele frequency to be less than 0.05, leaving us with 644 variants. Next, we filtered the variants by requiring minor allele frequency to be less than 0.05 in the UK10K data, leaving us with 139 variants. Lastly, we identified pathogenic or likely pathogenic variants using ClinVar significance, leaving us with 65 candidate variants that could be associated with different types of cancer (Table S4).

To check for the somatic status of these 65 identified SNVs, we examined the COSMIC database and identified 31 out of 65 SNVs as somatic, which matched 383 COSMIC records (some SNVs were observed in multiple samples). Interestingly, these identified somatic SNVs were predicted by FATHMM as potentially pathogenic, which is a promising starting point for future functional and clinical studies interested in cancer driver mutations (Table S5).

3.4. Biological Pathway Analysis

At variant level, among the 24 SNVs identified from dbMTS, we found out that one of the SNVs, rs3044, had been reported to be associated with colorectal cancer [ref. 17]. Another SNV, ENST00000467482.5:c.*134C>T, was observed in three patients with lung cancer in COSMIC. The functions of the other SNVs need further biological validation to determine if they can contribute to cancer susceptibility.

At gene level, we obtained 23 different genes that contain the 24 identified variants from dbMTS (gene XIAP has 2 variants). We used these 23 genes for cancer and disease association analysis using KEGG pathway annotations. We observed that two diseases are associated with gene BMPR1A (Juvenile polyposis syndrome [ref. 18], Hereditary mixed polyposis syndrome [ref. 19]), which causes an increased risk for colorectal carcinoma. In addition, X-linked lymphoproliferative syndrome is associated with gene XIAP [ref. 20], which causes lymphomas in approximately one-third of all patients.

Additionally, using KEGG Pathways, we checked 33 different genes containing 65 variants from dbNSFP for cancer and disease association. KEGG results indicated that several of these 33 genes had cancer-related disease associations. Diseases associated with genes SDHB and SDHD (Cowden Syndrome [ref. 21], Malignant paraganglioma [ref. 22]), TSC1 and TSC2 (Tuberous sclerosis complex [ref. 23]), ERCC6 (Disorders of nucleotide excision repair [ref. 24]), FLCN (Birt–Hogg–Dube Syndrome [ref. 25]), and CHEK2 (Li-Fraumeni Syndrome [ref. 26]) increase risk for various cancer types. Mismatch repair deficiency, associated with gene EPCAM [ref. 27], and familial adenomatous polyposis, associated with gene APC [ref. 28], can potentially increase risk for colorectal cancer. Nijmegen syndrome, associated with gene NBN [ref. 29], increases risk for lymphomas, and basal cell nevus syndrome, associated with gene PTCH1 [ref. 30], increases risk for basal cell carcinoma.

Among the 33 candidate genes identified from dbNSFP analysis, we found that there are nine genes that contain at least one cancer associated variant reported in COSMIC (Table 3). This result indicates that genes harboring multiple variants are likely to be associated with more cancer types, suggesting that some variants are cancer type specific.

3.5. Survival Analysis

We conducted a survival analysis of the two disease-associated genes from the pathway analysis, BMPR1A and XIAP. We first used an online tool, Gepia. For BMPR1A, among the three associated cancer types (Table 1), its high expression has significant survival benefits in KIRC (Kidney Renal Clear Cell Carcinoma) (p value: 0.00013, Figure 3). For XIAP, among the three associated cancer types (Table 1), its low expression has significant survival benefits in BRCA (Breast Invasive Carcinoma) (p value: 0.02, Figure 4). We also conducted survival analysis using another online tool, OncoLnc. The results were similar between these two different approaches, which indicate that the candidate genes we identified are potentially functional.

4. Conclusions

Given that two categories of variants (3′UTR and nonsynonymous) are important in changing gene expression, we analyzed the anti-correlated miRNA–mRNA target pairs and target genes obtained from TCGA dataset. We prioritized many candidate SNVs and genes harboring them using dbMTS, dbNSFP, and other functional annotation tools. We believe the disruption resulting from SNV in the miRNA target site is accountable for the observed anti-correlated miRNA-gene pairs, which contributes to tumorigenesis. Additionally, non-synonymous SNVs identified in dbNSFP can effectively change the amino acid sequence, which is also a possible cause for tumorigenesis. The candidate genes and candidate somatic and germline variants identified in this study could provide guidance for wet-lab scientists or clinical researchers in biomarker discovery. However, this will require extensive experimental work or validation from clinical data.

References

1. M.S. Nicoloso, H. Sun, R. Spizzo, H. Kim, P. Wickramasinghe, M. Shimizu, S.E. Wojcik, J. Ferdin, T. Kunej, L. Xiao. Single-Nucleotide Polymorphisms inside MicroRNA Target Sites Influence Tumor Susceptibility. Cancer Res., 2010. [DOI | PubMed]
2. X. Dai, L. Ding, H. Liu, Z. Xu, H. Jiang, S.K. Handelman, Y. Bai. Identifying Interaction Clusters for MiRNA and MRNA Pairs in TCGA Network. Genes, 2019. [DOI | PubMed]
3. C. Li, C. Mou, M.D. Swartz, B. Yu, Y. Bai, Y. Tu, X. Liu. DbMTS: A Comprehensive Database of Putative Human MicroRNA Target Site SNVs and Their Functional Predictions. Hum. Mutat., 2020. [DOI | PubMed]
4. X. Liu, X. Jian, E. Boerwinkle. DbNSFP: A Lightweight Database of Human Nonsynonymous SNPs and Their Functional Predictions. Hum. Mutat., 2011. [DOI | PubMed]
5. A. Buniello, J.A.L. Macarthur, M. Cerezo, L.W. Harris, J. Hayhurst, C. Malangone, A. McMahon, J. Morales, E. Mountjoy, E. Sollis. The NHGRI-EBI GWAS Catalog of Published Genome-Wide Association Studies, Targeted Arrays and Summary Statistics 2019. Nucleic Acids Res., 2019. [DOI | PubMed]
6. J.D. Eicher, C. Landowski, B. Stackhouse, A. Sloan, W. Chen, N. Jensen, J.P. Lien, R. Leslie, A.D. Johnson. GRASP v2.0: An Update on the Genome-Wide Repository of Associations between SNPs and Phenotypes. Nucleic Acids Res., 2015. [DOI | PubMed]
7. M.J. Landrum, J.M. Lee, G.R. Riley, W. Jang, W.S. Rubinstein, D.M. Church, D.R. Maglott. ClinVar: Public Archive of Relationships among Sequence Variation and Human Phenotype. Nucleic Acids Res., 2014. [DOI | PubMed]
8. W. McLaren, L. Gil, S.E. Hunt, H.S. Riat, G.R.S. Ritchie, A. Thormann, P. Flicek, F. Cunningham. The Ensembl Variant Effect Predictor. Genome Biol., 2016. [DOI | PubMed]
9. A. Bhattacharya, Y. Cui. SomamiR 2.0: A Database of Cancer Somatic Mutations Altering MicroRNA-CeRNA Interactions. Nucleic Acids Res., 2016. [DOI | PubMed]
10. J.G. Tate, S. Bamford, H.C. Jubb, Z. Sondka, D.M. Beare, N. Bindal, H. Boutselakis, C.G. Cole, C. Creatore, E. Dawson. COSMIC: The Catalogue of Somatic Mutations in Cancer. Nucleic Acids Res., 2019. [DOI | PubMed]
11. H. Ogata, S. Goto, K. Sato, W. Fujibuchi, H. Bono, M. Kanehisa. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res., 1999. [DOI | PubMed]
12. Z. Tang, C. Li, B. Kang, G. Gao, C. Li, Z. Zhang. GEPIA: A Web Server for Cancer and Normal Gene Expression Profiling and Interactive Analyses. Nucleic Acids Res., 2017. [DOI | PubMed]
13. J. Anaya. OncoLnc: Linking TCGA Survival Data to MRNAs, MiRNAs, and LncRNAs. PeerJ Comput. Sci., 2016. [DOI]
14. X. Jia, C. Wang, H. Huang, P. Zhang, Z. Yao, L. Xu. FLAD1 Is Overexpressed in Breast Cancer and Is a Potential Predictor of Prognosis and Treatment. Int. J. Clin. Exp. Med., 2019
15. C.K. Thammaiah, S. Jayaram. Role of Let-7 Family MicroRNA in Breast Cancer. Non-Coding RNA Res., 2016. [DOI | PubMed]
16. Y. Xiong, Y. Wang, L. Wang, Y. Huang, Y. Xu, L. Xu, Y. Guo, J. Lu, X. Li, M. Zhu. MicroRNA-30b Targets Snail to Impede Epithelial-Mesenchymal Transition in Pancreatic Cancer Stem Cells. J. Cancer, 2018. [DOI | PubMed]
17. R. Leslie, C.J. O’Donnell, A.D. Johnson. GRASP: Analysis of Genotype-Phenotype Results from 1390 Genome-Wide Association Studies and Corresponding Open Access Database. Bioinformatics, 2014. [DOI | PubMed]
18. L.A.A. Brosens, D. Langeveld, W.A. van Hattem, F.M. Giardiello, G.J.A. Offerhaus. Juvenile Polyposis Syndrome. World J. Gastroenterol., 2011. [DOI | PubMed]
19. X. Cao, K.W. Eu, M.P. Kumarasinghe, H.H. Li, C. Loi, P.Y. Cheah. Mapping of Hereditary Mixed Polyposis Syndrome (HMPS) to Chromosome 10q23 by Genomewide High-Density Single Nucleotide Polymorphism (SNP) Scan and Identification of BMPR1A Loss of Function. J. Med. Genet., 2005. [DOI]
20. S. Rigaud, M.C. Fondanèche, N. Lambert, B. Pasquier, V. Mateo, P. Soulas, L. Galicier, F. Le Deist, F. Rieux-Laucat, P. Revy. XIAP Deficiency in Humans Causes an X-Linked Lymphoproliferative Syndrome. Nature, 2006. [DOI | PubMed]
21. Y. Ni, K.M. Zbuk, T. Sadler, A. Patocs, G. Lobo, E. Edelman, P. Platzer, M.S. Orloff, K.A. Waite, C. Eng. Germline Mutations and Variants in the Succinate Dehydrogenase Genes in Cowden and Cowden-like Syndromes. Am. J. Hum. Genet., 2008. [DOI | PubMed]
22. G. Parenti, B. Zampetti, E. Rapizzi, T. Ercolino, V. Giach, M. Mannelli. Updated and New Perspectives on Diagnosis, Prognosis, and Therapy of Malignant Pheochromocytoma/Paraganglioma. J. Oncol., 2012. [DOI | PubMed]
23. J. Borkowska, R.A. Schwartz, K. Kotulska, S. Jozwiak. Tuberous Sclerosis Complex: Tumors and Tumorigenesis. Int. J. Dermatol., 2010. [DOI | PubMed]
24. J.E. Cleaver, E.T. Lam, I. Revet. Disorders of Nucleotide Excision Repair: The Genetic and Molecular Basis of Heterogeneity. Nat. Rev. Genet., 2009. [DOI | PubMed]
25. F.H. Menko, M. Van Steensel, S. Giraud, L.J. Friis-Hansen, S. Richard, S. Ungari, M. Nordenskjöld, T.V. Hansen, J. Solly, E. Maher. Birt-Hogg-Dubé Syndrome: Diagnosis and Management. Lancet Oncol., 2009. [DOI | PubMed]
26. L.L. Bachinski, S.E. Olufemi, X. Zhou, C.C. Wu, L. Yip, S. Shete, G. Lozano, C.I. Amos, L.C. Strong, R. Krahe. Genetic Mapping of a Third Li-Fraumeni Syndrome Predisposition Locus to Human Chromosome 1q23. Cancer Res., 2005. [PubMed]
27. H.T. Lynch, D.L. Riegert-Johnson, C. Snyder, J.F. Lynch, J. Hagenkord, C.R. Boland, J. Rhees, S.N. Thibodeau, L.A. Boardman, J. Davies. Lynch Syndrome-Associated Extracolonic Tumors Are Rare in Two Extended Families with the Same EPCAM Deletion. Am. J. Gastroenterol., 2011. [DOI | PubMed]
28. E. Half, D. Bercovich, P. Rozen. Familial Adenomatous Polyposis. Orphanet J. Rare Dis., 2009. [DOI | PubMed]
29. R. Waltes, R. Kalb, M. Gatei, A.W. Kijas, M. Stumm, A. Sobeck, B. Wieland, R. Varon, Y. Lerenthal, M.F. Lavin. Human RAD50 Deficiency in a Nijmegen Breakage Syndrome-like Disorder. Am. J. Hum. Genet., 2009. [DOI | PubMed]
30. C. Takahashi, N. Kanazawa, Y. Yoshikawa, R. Yoshikawa, Y. Saitoh, H. Chiyo, T. Tanizawa, T. Hashimoto-Tamaoki, Y. Nakano. Germline PTCH1 Mutations in Japanese Basal Cell Nevus Syndrome Patients. J. Hum. Genet., 2009. [DOI | PubMed]