GDCRNATools: integrative analysis of protein coding genes, long non-coding genes, and microRNAs in GDC

Ruidong Li and Han Qu

Last update: 04 January, 2019

1 Introduction

GDCRNATools is an R package which provides a standard, easy-to-use and comprehensive pipeline for downloading, organizing, and integrative analyzing RNA expression data in the GDC portal with an emphasis on deciphering the lncRNA-mRNA related ceRNAs regulatory network in cancer.

Competing endogenous RNAs (ceRNAs) are RNAs that indirectly regulate other transcripts by competing for shared miRNAs. Although only a fraction of long non-coding RNAs has been functionally characterized, increasing evidences show that lncRNAs harboring multiple miRNA response elements (MREs) can act as ceRNAs to sequester miRNA activity and thus reduce the inhibition of miRNA on its targets. Deregulation of ceRNAs network may lead to human diseases.

The Genomic Data Commons (GDC) maintains standardized genomic, clinical, and biospecimen data from National Cancer Institute (NCI) programs including The Cancer Genome Atlas (TCGA) and Therapeutically Applicable Research To Generate Effective Treatments (TARGET), It also accepts high quality datasets from non-NCI supported cancer research programs, such as genomic data from the Foundation Medicine.

Many analyses can be perfomed using GDCRNATools, including differential gene expression analysis (limma(Ritchie et al. 2015), edgeR(Robinson, McCarthy, and Smyth 2010), and DESeq2(Love, Huber, and Anders 2014)), univariate survival analysis (CoxPH and KM), competing endogenous RNA network analysis (hypergeometric test, Pearson correlation analysis, regulation similarity analysis, sensitivity Pearson partial correlation(Paci, Colombo, and Farina 2014)), and functional enrichment analysis(GO, KEGG, DO). Besides some routine visualization methods such as volcano plot, scatter plot, and bubble plot, etc., three simple shiny apps are developed in GDCRNATools allowing users visualize the results on a local webpage. All the figures are plotted based on ggplot2 package unless otherwise specified.

This user-friendly package allows researchers perform the analysis by simply running a few functions and integrate their own pipelines such as molecular subtype classification, weighted correlation network analysis (WGCNA)(Langfelder and Horvath 2008), and TF-miRNA co-regulatory network analysis, etc. into the workflow easily. This could open a door to accelerate the study of crosstalk among different classes of RNAs and their regulatory relationships in cancer.

2 Citation

If you use GDCRNATools in your published research, please cite:

Li, R., Qu, H., Wang, S., Wei, J., Zhang, L., Ma, R., Lu, J., Zhu, J., Zhong, W., and Jia, Z. (2018). GDCRNATools: an R/Bioconductor package for integrative analysis of lncRNA, miRNA and mRNA data in GDC. Bioinformatics 34, 2515-2517. https://doi.org/10.1093/bioinformatics/bty124.

4 Quick start

In GDCRNATools, some functions are built for users to download and process GDC data efficiently. Users can also use their own data that is processed by other tools such as the UCSC Xena GDC hub, TCGAbiolinks(Colaprico et al. 2016), or TCGA-Assembler(Zhu, Qiu, and Ji 2014), etc.

Here we use a small dataset to show the most basic steps for ceRNAs network analysis. More detailed instruction of each step is in the Case Study section.

4.1 Data preparation

4.1.2 Parse metadata

4.2 ceRNAs network analysis

4.2.1 Identification of differentially expressed genes (DEGs)

4.2.2 ceRNAs network analysis of DEGs

## Step 1/3: Hypergenometric test done !
## Step 2/3: Correlation analysis done !
## Step 3/3: Regulation pattern analysis done !

4.2.3 Export ceRNAs network to Cytoscape

5 Case study: TCGA-CHOL

In this section, we use the whole datasets of TCGA-CHOL project as an example to illustrate how GDCRNATools works in detail.

5.1 Data download

Two methods are provided for downloading Gene Expression Quantification (HTSeq-Counts), Isoform Expression Quantification (BCGSC miRNA Profiling), and Clinical (Clinical Supplement) data:

5.1.2 Manual download

Users can also download data manually by providing the manifest file that is downloaded from the GDC cart

Step1: Download GDC Data Transfer Tool on the GDC website
Step2: Add data to the GDC cart, then download manifest file and metadata of the cart
Step3: Download data using gdcRNADownload() or gdcClinicalDownload() functions by providing the manifest file

5.2 Data organization and DE analysis

5.2.1 Parse metadata

Metadata can be parsed by either providing the metadata file (.json) that is downloaded in the data download step, or specifying the project.id and data.type in gdcParseMetadata() function to obtain information of data in the manifest file to facilitate data organization and basic clinical information of patients such as age, stage and gender, etc. for data analysis.

Only one sample would be kept if the sample had been sequenced more than once by gdcFilterDuplicate(). Samples that are neither Primary Tumor (code: 01) nor Solid Tissue Normal (code: 11) would be filtered out by gdcFilterSampleType()

5.2.2 Merge raw counts data

gdcRNAMerge() merges raw counts data of RNAseq to a single expression matrix with rows are Ensembl id and columns are samples. Total read counts for 5p and 3p strands of miRNAs can be processed from isoform quantification files and then merged to a single expression matrix with rows are miRBase v21 identifiers and columns are samples. If data for different samples are in separate folders, specify organized=FALSE, otherwise, specify organized=TRUE.

5.2.3 Merge clinical data

By specifying key.info=TRUE, only common clinical information will be organized and reported. Otherwise, all the clinical information from the XML files will be extracted.

5.2.4 TMM normalization and voom transformation

By running gdcVoomNormalization() function, raw counts data would be normalized by TMM method implemented in edgeR(Robinson, McCarthy, and Smyth 2010) and further transformed by the voom method provided in limma(Ritchie et al. 2015). Low expression genes (logcpm < 1 in more than half of the samples) will be filtered out by default. All the genes can be kept by setting filter=TRUE in the gdcVoomNormalization().

5.2.5 Differential gene expression analysis

Usually, people are interested in genes that are differentially expressed between different groups (eg. Primary Tumor vs. Solid Tissue Normal). gdcDEAnalysis(), a convenience wrapper, provides three widely used methods limma(Ritchie et al. 2015), edgeR(Robinson, McCarthy, and Smyth 2010), and DESeq2(Love, Huber, and Anders 2014) to identify differentially expressed genes (DEGs) or miRNAs between any two groups defined by users. Note that DESeq2(Love, Huber, and Anders 2014) maybe slow with a single core. Multiple cores can be specified with the nCore argument if DESeq2(Love, Huber, and Anders 2014) is in use. Users are encouraged to consult the vignette of each method for more detailed information.

All DEGs, DE long non-coding genes, DE protein coding genes and DE miRNAs could be reported separately by setting geneType argument in gdcDEReport(). Gene symbols and biotypes based on the Ensembl 90 annotation are reported in the output.

5.2.6 DEG visualization

Volcano plot and bar plot are used to visualize DE analysis results in different manners by gdcVolcanoPlot() and gdcBarPlot() functions, respectively . Hierarchical clustering on the expression matrix of DEGs can be analyzed and plotted by the gdcHeatmap() function.

5.2.6.1 Volcano plot

5.2.6.3 Heatmap

Heatmap is generated based on the heatmap.2() function in gplots package.

5.3 Competing endogenous RNAs network analysis

Three criteria are used to determine the competing endogenous interactions between lncRNA-mRNA pairs:

  • The lncRNA and mRNA must share significant number of miRNAs
  • Expression of lncRNA and mRNA must be positively correlated
  • Those common miRNAs should play similar roles in regulating the expression of lncRNA and mRNA

5.3.1 Hypergeometric test

Hypergenometric test is performed to test whether a lncRNA and mRNA share many miRNAs significantly.

A newly developed algorithm spongeScan(Furi’o-Tar’i et al. 2016) is used to predict MREs in lncRNAs acting as ceRNAs. Databases such as starBase v2.0(Li et al. 2014), miRcode(Jeggari, Marks, and Larsson 2012) and mirTarBase release 7.0(Chou et al. 2017) are used to collect predicted and experimentally validated miRNA-mRNA and/or miRNA-lncRNA interactions. Gene IDs in these databases are updated to the latest Ensembl 90 annotation of human genome and miRNAs names are updated to the new release miRBase 21 identifiers. Users can also provide their own datasets of miRNA-lncRNA and miRNA-mRNA interactions.

The figure and equation below illustrate how the hypergeometric test works

\[p=1-\sum_{k=0}^m \frac{\binom{K}{k}\binom{N-K}{n-k}}{\binom{N}{n}} \] here \(m\) is the number of shared miRNAs, \(N\) is the total number of miRNAs in the database, \(n\) is the number of miRNAs targeting the lncRNA, \(K\) is the number of miRNAs targeting the protein coding gene.

5.3.2 Pearson correlation analysis

Pearson correlation coefficient is a measure of the strength of a linear association between two variables. As we all know, miRNAs are negative regulators of gene expression. If more common miRNAs are occupied by a lncRNA, less of them will bind to the target mRNA, thus increasing the expression level of mRNA. So expression of the lncRNA and mRNA in a ceRNA pair should be positively correlated.

5.3.3 Regulation pattern analysis

Two methods are used to measure the regulatory role of miRNAs on the lncRNA and mRNA:

  • Regulation similarity

We defined a measurement regulation similarity score to check the similarity between miRNAs-lncRNA expression correlation and miRNAs-mRNA expression correlation.

\[Regulation\ similarity\ score = 1-\frac{1}{M} \sum_{k=1}^M [{\frac{|corr(m_k,l)-corr(m_k,g)|}{|corr(m_k,l)|+|corr(m_k,g)|}}]^M\]

where \(M\) is the total number of shared miRNAs, \(k\) is the \(k\)th shared miRNAs, \(corr(m_k, l)\) and \(corr(m_k, g)\) represents the Pearson correlation between the \(k\)th miRNA and lncRNA, the \(k\)th miRNA and mRNA, respectively

  • Sensitivity correlation

Sensitivity correlation is defined by Paci et al.(Paci, Colombo, and Farina 2014) to measure if the correlation between a lncRNA and mRNA is mediated by a miRNA in the lncRNA-miRNA-mRNA triplet. We take average of all triplets of a lncRNA-mRNA pair and their shared miRNAs as the sensitivity correlation between a selected lncRNA and mRNA.

\[Sensitivity\ correlation = corr(l,g)-\frac{1}{M}\sum_{k=1}^M {\frac{corr(l,g)-corr(m_k,l)corr(m_k,g)}{\sqrt{1-corr(m_k,l)^2}\sqrt{1-corr(m_k,g)^2}}}\] where \(M\) is the total number of shared miRNAs, \(k\) is the \(k\)th shared miRNAs, \(corr(l,g)\), \(corr(m_k,l)\) and \(corr(m_k, g)\) represents the Pearson correlation between the long non-coding RNA and the protein coding gene, the kth miRNA and lncRNA, the kth miRNA and mRNA, respectively

5.3.4 ceRNAs network analysis

The hypergeometric test of shared miRNAs, expression correlation analysis of lncRNA-mRNA pair, and regulation pattern analysis of shared miRNAs are all implemented in the gdcCEAnalysis() function.

5.3.4.1 ceRNAs network analysis using internal databases

Users can use the internally incoporated databases of miRNA-mRNA (starBase v2.0, miRcode, and mirTarBase v7.0) and miRNA-lncRNA (starBase v2.0, miRcode, spongeScan) interactions to perform the ceRNAs network analysis.

5.3.4.2 ceRNAs network analysis using user-provided datasets

gdcCEAnalysis() can also take user-provided miRNA-mRNA and miRNA-lncRNA interaction datasets, such as miRNA-target interactions predicted by TargetScan, miRanda, and Diana Tools, etc. for the ceRNAs network analysis.

## $ENSG00000138829
##  [1] "hsa-miR-200b-3p" "hsa-miR-429"     "hsa-miR-101-3p" 
##  [4] "hsa-miR-137"     "hsa-miR-9-5p"    "hsa-miR-139-5p" 
##  [7] "hsa-miR-200c-3p" "hsa-miR-136-5p"  "hsa-miR-494-3p" 
## [10] "hsa-miR-495-3p"  "hsa-miR-154-5p"  "hsa-miR-410-3p" 
## [13] "hsa-miR-211-5p"  "hsa-miR-140-5p"  "hsa-miR-22-3p"  
## [16] "hsa-miR-33b-5p"  "hsa-miR-144-3p"  "hsa-miR-133a-3p"
## [19] "hsa-miR-23a-3p"  "hsa-miR-217"     "hsa-miR-33a-5p" 
## [22] "hsa-miR-218-5p"  "hsa-miR-133b"    "hsa-miR-876-5p" 
## [25] "hsa-miR-204-5p"  "hsa-miR-23b-3p"  "hsa-miR-23c"    
## 
## $ENSG00000113615
##  [1] "hsa-miR-200b-3p" "hsa-miR-200a-3p" "hsa-miR-429"    
##  [4] "hsa-miR-30e-5p"  "hsa-miR-30c-5p"  "hsa-miR-92b-3p" 
##  [7] "hsa-miR-199a-5p" "hsa-miR-181b-5p" "hsa-miR-181a-5p"
## [10] "hsa-miR-107"     "hsa-miR-200c-3p" "hsa-miR-141-3p" 
## [13] "hsa-miR-26a-5p"  "hsa-miR-16-5p"   "hsa-miR-15a-5p" 
## [16] "hsa-miR-1297"    "hsa-miR-92a-3p"  "hsa-miR-136-5p" 
## [19] "hsa-miR-300"     "hsa-miR-381-3p"  "hsa-miR-539-5p" 
## [22] "hsa-miR-7-5p"    "hsa-miR-132-3p"  "hsa-miR-212-3p" 
## [25] "hsa-miR-195-5p"  "hsa-miR-497-5p"  "hsa-miR-144-3p" 
## [28] "hsa-miR-27a-3p"  "hsa-miR-23a-3p"  "hsa-miR-181c-5p"
## [31] "hsa-miR-181d-5p" "hsa-miR-371a-5p" "hsa-miR-128-3p" 
## [34] "hsa-miR-26b-5p"  "hsa-miR-103a-3p" "hsa-miR-15b-5p" 
## [37] "hsa-miR-367-3p"  "hsa-miR-30a-5p"  "hsa-miR-653-5p" 
## [40] "hsa-miR-25-3p"   "hsa-miR-182-5p"  "hsa-miR-183-5p" 
## [43] "hsa-miR-490-3p"  "hsa-miR-30b-5p"  "hsa-miR-30d-5p" 
## [46] "hsa-miR-31-5p"   "hsa-miR-23b-3p"  "hsa-miR-27b-3p" 
## [49] "hsa-miR-32-5p"   "hsa-miR-199b-5p" "hsa-miR-23c"    
## [52] "hsa-miR-374b-5p" "hsa-miR-374a-5p" "hsa-miR-363-3p" 
## [55] "hsa-miR-424-5p" 
## 
## $ENSG00000112144
##  [1] "hsa-miR-200b-3p" "hsa-miR-429"     "hsa-miR-30e-5p" 
##  [4] "hsa-miR-30c-5p"  "hsa-miR-101-3p"  "hsa-miR-202-3p" 
##  [7] "hsa-miR-139-5p"  "hsa-miR-200c-3p" "hsa-miR-26a-5p" 
## [10] "hsa-miR-1297"    "hsa-miR-543"     "hsa-miR-300"    
## [13] "hsa-miR-382-5p"  "hsa-miR-410-3p"  "hsa-miR-144-3p" 
## [16] "hsa-miR-23a-3p"  "hsa-miR-217"     "hsa-miR-26b-5p" 
## [19] "hsa-miR-218-5p"  "hsa-miR-367-3p"  "hsa-miR-30a-5p" 
## [22] "hsa-miR-383-5p"  "hsa-miR-30b-5p"  "hsa-miR-30d-5p" 
## [25] "hsa-miR-23b-3p"  "hsa-miR-374b-5p" "hsa-miR-374a-5p"
## [28] "hsa-miR-448"

5.3.7 Correlation plot on a local webpage

shinyCorPlot() , a interactive plot function based on shiny package, can be easily operated by just clicking the genes in each drop down box (in the GUI window). By running shinyCorPlot() function, a local webpage would pop up and correlation plot between a lncRNA and mRNA would be automatically shown.

5.4 Other downstream analyses

Downstream analyses such as univariate survival analysis and functional enrichment analysis are developed in the GDCRNATools package to facilitate the identification of genes in the ceRNAs network that play important roles in prognosis or involve in important pathways.

5.4.1 Univariate survival analysis

Two methods are provided to perform univariate survival analysis: Cox Proportional-Hazards (CoxPH) model and Kaplan Meier (KM) analysis based on the survival package. CoxPH model considers expression value as continous variable while KM analysis divides patients into high-expreesion and low-expression groups by a user-defined threshold such as median or mean. gdcSurvivalAnalysis() take a list of genes as input and report the hazard ratio, 95% confidence intervals, and test significance of each gene on overall survival.

5.4.1.4 KM plot on a local webpage by shinyKMPlot

The shinyKMPlot() function is also a simple shiny app which allow users view KM plots (based on the R package survminer.) of all genes of interests on a local webpackage conveniently.

5.4.2 Functional enrichment analysis

gdcEnrichAnalysis() can perform Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and Disease Ontology (DO) functional enrichment analyses of a list of genes simultaneously. GO and KEGG analyses are based on the R/Bioconductor packages clusterProfilier(Yu et al. 2012) and DOSE(Yu et al. 2015). Redundant GO terms can be removed by specifying simplify=TRUE in the gdcEnrichAnalysis() function which uses the simplify() function in the clusterProfilier(Yu et al. 2012) package.

5.4.2.3 View pathway maps on a local webpage

shinyPathview() allows users view and download pathways of interests by simply selecting the pathway terms on a local webpage.

##  [1] "hsa05414~Dilated cardiomyopathy (DCM)"                          
##  [2] "hsa05410~Hypertrophic cardiomyopathy (HCM)"                     
##  [3] "hsa05412~Arrhythmogenic right ventricular cardiomyopathy (ARVC)"
##  [4] "hsa04512~ECM-receptor interaction"                              
##  [5] "hsa04510~Focal adhesion"                                        
##  [6] "hsa04360~Axon guidance"                                         
##  [7] "hsa04270~Vascular smooth muscle contraction"                    
##  [8] "hsa05205~Proteoglycans in cancer"                               
##  [9] "hsa04022~cGMP-PKG signaling pathway"                            
## [10] "hsa00480~Glutathione metabolism"

6 sessionInfo

## R version 3.5.2 (2018-12-20)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server 2012 R2 x64 (build 9600)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=C                          
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] DT_0.5            GDCRNATools_1.2.1
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.1.3             Hmisc_4.1-1                
##   [3] fastmatch_1.1-0             plyr_1.8.4                 
##   [5] igraph_1.2.2                lazyeval_0.2.1             
##   [7] splines_3.5.2               crosstalk_1.0.0            
##   [9] BiocParallel_1.16.5         pathview_1.22.1            
##  [11] GenomeInfoDb_1.18.1         ggplot2_3.1.0              
##  [13] urltools_1.7.1              digest_0.6.18              
##  [15] htmltools_0.3.6             GOSemSim_2.8.0             
##  [17] viridis_0.5.1               GO.db_3.7.0                
##  [19] gdata_2.18.0                magrittr_1.5               
##  [21] checkmate_1.8.5             memoise_1.1.0              
##  [23] cluster_2.0.7-1             limma_3.38.3               
##  [25] Biostrings_2.50.2           readr_1.3.1                
##  [27] annotate_1.60.0             matrixStats_0.54.0         
##  [29] enrichplot_1.2.0            prettyunits_1.0.2          
##  [31] colorspace_1.3-2            blob_1.1.1                 
##  [33] rappdirs_0.3.1              ggrepel_0.8.0              
##  [35] xfun_0.4                    dplyr_0.7.8                
##  [37] crayon_1.3.4                RCurl_1.95-4.11            
##  [39] jsonlite_1.6                graph_1.60.0               
##  [41] genefilter_1.64.0           bindr_0.1.1                
##  [43] zoo_1.8-4                   survival_2.43-3            
##  [45] glue_1.3.0                  survminer_0.4.3            
##  [47] GenomicDataCommons_1.6.0    gtable_0.2.0               
##  [49] zlibbioc_1.28.0             XVector_0.22.0             
##  [51] UpSetR_1.3.3                DelayedArray_0.8.0         
##  [53] Rgraphviz_2.26.0            BiocGenerics_0.28.0        
##  [55] scales_1.0.0                DOSE_3.8.0                 
##  [57] DBI_1.0.0                   edgeR_3.24.3               
##  [59] Rcpp_1.0.0                  cmprsk_2.2-7               
##  [61] viridisLite_0.3.0           xtable_1.8-3               
##  [63] progress_1.2.0              htmlTable_1.13             
##  [65] units_0.6-2                 gridGraphics_0.3-0         
##  [67] foreign_0.8-71              bit_1.1-14                 
##  [69] europepmc_0.3               km.ci_0.5-2                
##  [71] Formula_1.2-3               stats4_3.5.2               
##  [73] htmlwidgets_1.3             httr_1.4.0                 
##  [75] fgsea_1.8.0                 gplots_3.0.1               
##  [77] RColorBrewer_1.1-2          acepack_1.4.1              
##  [79] pkgconfig_2.0.2             XML_3.98-1.16              
##  [81] farver_1.1.0                nnet_7.3-12                
##  [83] locfit_1.5-9.1              labeling_0.3               
##  [85] ggplotify_0.0.3             tidyselect_0.2.5           
##  [87] rlang_0.3.0.1               reshape2_1.4.3             
##  [89] later_0.7.5                 AnnotationDbi_1.44.0       
##  [91] munsell_0.5.0               tools_3.5.2                
##  [93] generics_0.0.2              RSQLite_2.1.1              
##  [95] broom_0.5.1                 ggridges_0.5.1             
##  [97] evaluate_0.12               stringr_1.3.1              
##  [99] yaml_2.2.0                  org.Hs.eg.db_3.7.0         
## [101] knitr_1.21                  bit64_0.9-7                
## [103] survMisc_0.5.5              caTools_1.17.1.1           
## [105] purrr_0.2.5                 KEGGREST_1.22.0            
## [107] ggraph_1.0.2                bindrcpp_0.2.2             
## [109] nlme_3.1-137                mime_0.6                   
## [111] KEGGgraph_1.42.0            DO.db_2.9                  
## [113] xml2_1.2.0                  biomaRt_2.38.0             
## [115] compiler_3.5.2              rstudioapi_0.8             
## [117] png_0.1-7                   tibble_2.0.0               
## [119] tweenr_1.0.1                geneplotter_1.60.0         
## [121] stringi_1.2.4               lattice_0.20-38            
## [123] Matrix_1.2-15               KMsurv_0.1-5               
## [125] pillar_1.3.1                triebeard_0.3.0            
## [127] data.table_1.11.8           cowplot_0.9.3              
## [129] bitops_1.0-6                httpuv_1.4.5.1             
## [131] GenomicRanges_1.34.0        qvalue_2.14.0              
## [133] R6_2.3.0                    latticeExtra_0.6-28        
## [135] promises_1.0.1              KernSmooth_2.23-15         
## [137] gridExtra_2.3               IRanges_2.16.0             
## [139] gtools_3.8.1                MASS_7.3-51.1              
## [141] assertthat_0.2.0            SummarizedExperiment_1.12.0
## [143] rjson_0.2.20                DESeq2_1.22.2              
## [145] S4Vectors_0.20.1            GenomeInfoDbData_1.2.0     
## [147] parallel_3.5.2              hms_0.4.2                  
## [149] clusterProfiler_3.10.1      grid_3.5.2                 
## [151] prettydoc_0.2.1             rpart_4.1-13               
## [153] tidyr_0.8.2                 rmarkdown_1.11             
## [155] rvcheck_0.1.3               ggpubr_0.2                 
## [157] ggforce_0.1.3               Biobase_2.42.0             
## [159] shiny_1.2.0                 base64enc_0.1-3

References

Chou, Chih-Hung, Sirjana Shrestha, Chi-Dung Yang, Nai-Wen Chang, Yu-Ling Lin, Kuang-Wen Liao, Wei-Chi Huang, et al. 2017. “miRTarBase Update 2018: A Resource for Experimentally Validated microRNA-Target Interactions.” Nucleic Acids Research, November, gkx1067–gkx1067. https://doi.org/10.1093/nar/gkx1067.

Colaprico, Antonio, Tiago C. Silva, Catharina Olsen, Luciano Garofano, Claudia Cava, Davide Garolini, Thais S. Sabedot, et al. 2016. “TCGAbiolinks: An R/Bioconductor Package for Integrative Analysis of TCGA Data.” Nucleic Acids Research 44 (8):e71. https://doi.org/10.1093/nar/gkv1507.

Furi’o-Tar’i, Pedro, Sonia Tarazona, Toni Gabald’on, Anton J. Enright, and Ana Conesa. 2016. “spongeScan: A Web for Detecting microRNA Binding Elements in lncRNA Sequences.” Nucleic Acids Research 44 (Web Server issue):W176–W180. https://doi.org/10.1093/nar/gkw443.

Jeggari, Ashwini, Debora S Marks, and Erik Larsson. 2012. “miRcode: A Map of Putative microRNA Target Sites in the Long Non-Coding Transcriptome.” Bioinformatics 28 (15):2062–3. https://doi.org/10.1093/bioinformatics/bts344.

Langfelder, Peter, and Steve Horvath. 2008. “WGCNA: An R Package for Weighted Correlation Network Analysis.” BMC Bioinformatics 9 (December):559. https://doi.org/10.1186/1471-2105-9-559.

Li, Jun-Hao, Shun Liu, Hui Zhou, Liang-Hu Qu, and Jian-Hua Yang. 2014. “starBase V2.0: Decoding miRNA-ceRNA, miRNA-ncRNA and Protein–RNA Interaction Networks from Large-Scale CLIP-Seq Data.” Nucleic Acids Research 42 (Database issue):D92–D97. https://doi.org/10.1093/nar/gkt1248.

Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2.” Genome Biology 15 (December):550. https://doi.org/10.1186/s13059-014-0550-8.

Paci, Paola, Teresa Colombo, and Lorenzo Farina. 2014. “Computational Analysis Identifies a Sponge Interaction Network Between Long Non-Coding RNAs and Messenger RNAs in Human Breast Cancer.” BMC Systems Biology 8 (July):83. https://doi.org/10.1186/1752-0509-8-83.

Ritchie, Matthew E., Belinda Phipson, Di Wu, Yifang Hu, Charity W. Law, Wei Shi, and Gordon K. Smyth. 2015. “Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies.” Nucleic Acids Research 43 (7):e47. https://doi.org/10.1093/nar/gkv007.

Robinson, Mark D., Davis J. McCarthy, and Gordon K. Smyth. 2010. “edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1):139–40. https://doi.org/10.1093/bioinformatics/btp616.

Yu, Guangchuang, Li-Gen Wang, Yanyan Han, and Qing-Yu He. 2012. “clusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters.” OMICS : A Journal of Integrative Biology 16 (5):284–87. https://doi.org/10.1089/omi.2011.0118.

Yu, Guangchuang, Li-Gen Wang, Guang-Rong Yan, and Qing-Yu He. 2015. “DOSE: An R/Bioconductor Package for Disease Ontology Semantic and Enrichment Analysis.” Bioinformatics 31 (4):608–9. https://doi.org/10.1093/bioinformatics/btu684.

Zhu, Yitan, Peng Qiu, and Yuan Ji. 2014. “TCGA-Assembler: An Open-Source Pipeline for TCGA Data Downloading, Assembling, and Processing.” Nature Methods 11 (6):599–600. https://doi.org/10.1038/nmeth.2956.