A New R Package for Categorizing Coding and Non-Coding Genes

2. Implementation

Statistical framework

Let $M=\{m_1,m_2,\dots ,m_X\}$, be a set of mutational catalogs for $X$ cancers, and $R=\left\{r_1,r_2,\dots ,r_Y\right\}$ be the set of $Y$ input genomic regions. Using all mutational catalogs, a p-value for the region $r_y$ in cancer $x$ is calculated using a one-sided Fisher’s exact test with the following contingency table:

Table 1. Contingency table for Fisher’s exact test.

#samples in $m_x$ that have a mutation in the region $r_y$.	#samples in $m_x$ that do not have a mutation in the region $r_y$.
#samples in all catalogs except $m_x$ that have a mutation in the region $r_y$.	#samples in all catalogs except $m_x$ that do not have a mutation in the region $r_y$.

Table 1. Contingency table for Fisher’s exact test.

#samples in $m_x$ that have a mutation in the region $r_y$.	#samples in $m_x$ that do not have a mutation in the region $r_y$.
#samples in all catalogs except $m_x$ that have a mutation in the region $r_y$.	#samples in all catalogs except $m_x$ that do not have a mutation in the region $r_y$.

This Fisher’s exact test (which is applied for all regions of interest in $R$) is followed by a permutation test (typically for 1,000,000 random permutations - can be changed by the user) to estimate the probability that any of the associations emerges by chance. In the permutation process, the order of sample IDs in the original mutational catalogs is shuffled to generate a list of simulated datasets of mutations. Then the Fisher’s exact test is applied to each simulated dataset to obtain a list of simulated P-values for each $r_y$ (the number of simulated P-values equals the number of permutations). Finally, each $r_y$ is determined to be significant if its original P-value (obtained from Fisher’s exact test on the original data) is less than the permutations P-value at confidence interval 99%. Users can define the significance level for the Fisher’s exact test and the confidence interval for the permutation test. Using this statistical framework, the user can extract a significant list of genes (from an initial list) that compared to other cancers, are recurrently and specifically mutated in cancer $x$.

Annotations scheme

Let $R=\{r_1,r_2,\dots ,r_Y\}$ be the set of $Y$ input genomic regions and $A=\{a_1,a_2,\dots ,a_Z\}$ be a set of $Z$ genomic annotation entries. $A$ can be indicated as one of five general types of annotation (e.g., tissue-specific histone marks, encode predicated promoters and enhancers, GWAS, eQTL) that attribute specific genomic features to a set of genomic regions. This package aides the user to interpret an arbitrary set of regions $R$ with an optional annotation $A$ (one of five general types).

In the case of annotation of type 1, which is a simple catalog of genomic regions such as the list of active enhancers by FANTOM5 (24) or genomic positions of disease-associated genomic variants by GWAS Catalog (25), the annotation output is a list designating the occurrence of overlap/the number of overlaps, of annotation entries with each $r_y$ regions.

The 2nd type of annotation assigns a category to each genomic annotation entries $a_z$, (e.g., the annotation of chromatin state segmentation (26)) and the result for each region $r_Y$ indicates the categories overlapping with $r_Y$ along with the percentages of overlaps.

In the case of 3rd annotation, a catalog assigns a single score to each annotation entry (such as the list of genomic regions enriched for histone modifications (27)) and the annotation result for each regions $r_Y$ will be an aggregated score and the percentage of its length that overlaps with annotation entries.

The 4th type of annotation assigns a group of IDs (which we call sub_ids) and their corresponding scores (which we call sub_scores) to each annotation entry (e.g., DNase clusters (28)). In this case, the output for each region $r_Y$ will represent the overlapping sub_ids and the corresponding aggregated sub_scores.

The 5th type annotation is defined for chromatin interaction datasets (such as Hi-C). In this case, the user provides two commentaries, one for chromatin interactions and one for target genomic regions, such as $T=\{t_1,t_2,\dots ,t_G\}$. As a result, the package will identify the interactions between entries of original genomic regions $R$ and target genomic regions $T$ (see supplementary materials for more details). For instance, the user can investigate the interactions between a set of lncRNAs (entries of $R$) that interact with a group of protein-coding genes (entries of $T$) through a dataset of Hi-C interactions in breast cancer. More details about the calculation of overlaps, scores, and sub_scores are explained in the supplementary materials.

3. Application of SomaGene

An application of SomaGene can be found in a recent comprehensive study aiming to prioritize non-coding RNA genes in the context of cancer Catalog (29). In this study, a list of 65,000 non-coding genes and the mutational catalogs of 19 cancer types downloaded from the ICGC consortia, were used as input to SomaGene for identifying those non-coding genes that significantly mutated in breast cancer samples. SomaGene identified 929 non-coding genes as significantly mutated genes in breast cancer samples (confidence interval 99% on 1,000,000 permutations provided by SomaGene). Interestingly, the candidate non-coding RNAs have significantly greater fraction of breast tissue related ENCODE enhancer and promoter marks (Figure 2a), FANTOM5 breast tissue differentially expressed enhancers (Figure 2b), ENCODE chromatin active histone marks (Figure 2c), and breast cancer associated genome-wide association SNPs (GWAS) (Figure 2d).

Also, we sorted non-coding RNAs based on their mutational p-value and repeated the enrichment analyses for the second, third and last set of non-coding RNAs (i.e., those non-coding RNAs that did not encompass any mutations). Interestingly, the enrichment for regulatory features have been also seen in the second and third lists of non-coding RNAs, but in much smaller fractions compared to the significant list of non-coding RNAs (Figure 2). As expected, no enrichment was seen for the last list of non-coding RNAs (Figure 2).

Finally, we used genomic and epigenetic annotations described in the above section and provided a prioritized list of non-coding RNAs with potential enhancer roles in breast cancer as well as five other cancers that can be accessed through our online resource http://ncrna.ictic.sharif.edu.

The medical and biological datasets are increasing rapidly. To analyse such big and complex data, artificial intelligence and integrative pipelines become most popular (5, 30-43). Here, we have developed SomaGene, a novel, user-friendly, interactive, open-access pipeline for identifying and annotating noncoding genes (or genomic coordinates) that significantly mutate in a cancer genome.