Identification of m6A-Related Biomarkers Associated with Prognosis of Colorectal Cancer.

BACKGROUND Colorectal cancer (CRC) is the second most deadly cancer in the world according to GLOBOCAN 2020 data. Accumulating evidence suggests that RNA methylation modification is also misregulated in human cancers and may be a potential ideal target for cancer treatment. MATERIAL AND METHODS m6A-related differentially expressed genes (DEGs) were identified from colon adenocarcinoma and rectum adenocarcinoma esophageal carcinoma patients with different pathological stages. The protein-protein interaction (PPI) network construction, Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of DEGs were conducted. Cox regression analysis was applied to the screening of m6A-related DEGs significantly associated with the overall survival (OS), and those selected genes were used for LASSO regression analysis to construct prognostic signature and calculate patients' risk scores. RESULTS We identified 673 m6A-related DEGs from CRC patients in different pathologic stages, and 146 of them were associated with OS. CTNNB1, TRIM37, RAB7A, CASC5/KNL1, CENPE, CCNB1, UBE2H, HSPA8, KIF1A, and FBXW4 were hub genes of the PPI network. Nine m6A-related genes were screened out to build the prognostic risk model. TNM stage, vascular invasion, and the risk score were independently related to the OS of CRC patients. CONCLUSIONS Nine candidate m6A-related mRNA biomarkers (LRRC17, NFKB1, NOS2, PCDHB2, RAB7A, RPS6KA1, RRNAD1, TLE6, and UBE2H) were found to be closely related to the clinicopathology and prognosis of colorectal cancer, indicating that they could be potential prognostic biomarkers for patients with colorectal cancer.

Background

Colorectal cancer, also known as bowel cancer, colon cancer, or rectal cancer, is a cancer formed by uncontrolled cell growth in the colon or the rectum (part of the large intestine), or in the appendix. According to GLOBOCAN 2020 data, colorectal cancer (CRC) is the second most deadly and third most prevalent cancer in the world [1]. The incidence of CRC worldwide has been steadily increasing, accounting for 10% of all cancer diagnoses [2]. By 2020, CRC was the second most deadly cancer in the world, with an estimated 1.9 million new cases and 935 000 deaths. Meanwhile, the global burden of CRC is expected to increase by 60%, with nearly 2.2 million new cases and 1.1 million deaths expected by 2030 [3]. About 20% of CRC patients present with synchronous metastases, most commonly in the liver, and up to 60% of the patients develop distant metastases within 5 years [4]. Therefore, it is urgent to determine genetic and environmental risk factors that may affect the prognosis of CRC patients in clinical practice. Recently, with the help of high-throughput sequencing technology, a breakthrough has been made in the identification of CRC biomarkers at the cellular and molecular levels, which would potentially improve the prognostic prediction accuracy and introduce new therapeutic targets for CRC patients.

Since the discovery of first structurally modified nucleoside, pseudouridine, in the 1950s, more than 150 different chemical modifications have been identified on cellular RNA so far [5]. m6A methylation is the most characteristic mRNA modification and has been the most extensively studied since its discovery [6,7]. The most common methylation modification in eukaryotic mRNA is N6-methyladenosine (m6A), which accounts for more than 80% of all RNA base methylation and exists in various species [8]. The abundance and effects of m6A on RNA depend on the dynamic interplay between its methyltransferase (“writers”, such as METTL3, METTL14, WTAP, KIAA1429, ZC3H13, and METTL16), binding protein (“readers”, such as YTH domain-containing proteins and MRB1) and demethylase (“erasers”, such as FTO and ALKBH5) [9]. Accumulating evidence suggests that RNA methylation modification is also misregulated in human cancers and may be a potential ideal target for cancer treatment [10]. m6A methylation modification affects multiple aspects of RNA metabolism, ranging from RNA processing, nuclear export, RNA translation, to decay [11]. In addition, research proved that m6A methylation modification of mRNA and non-coding RNA plays an important role in a variety of common cancers, including solid tumors and non-solid tumors, and regulates cell proliferation and migration in cancer by affecting the biological functions of cells, tumor cell differentiation, and homeostasis [12,13]. Therefore, by finding m6A-related genes and m6A RNA methylation modification sites in cancer, new therapeutic targets could be provided for cancer treatment. Recently, studies focused on the role of m6A-related genes and their methylation regulators have revealed that METTL3 interacts with the microprocessor protein DGCR8 and actively regulates the pri-miR221/222 process in an m6A-dependent manner, which possibly has a carcinogenic effect in bladder cancer [14]. The reduction of RNA m6A methylation can activate oncogenic Wnt/PI3K-Akt signaling and promote malignant phenotypes of gastric cancer cells [15].

In the study, signature analysis with clinical information was performed for determining the m6A-related genes expression profile of patients with COAD (colon adenocarcinoma) and READ (rectum adenocarcinoma esophageal carcinoma) from TCGA database and Gene Expression Omnibus (GEO) databases. LASSO and Cox regression aided identification of potential m6A-related genes to predict the survival of patients with colorectal cancer. We determined 9 prognostic m6A-related genes from TCGA dataset and further validated the model in GEO dataset. The results obtained in this study would help to predict the prognosis of CRC patients and improve personalized treatment and management.

Material and Methods

Datasets Acquisition

RNA expression data and clinical information of primary CRC tissues were downloaded from TCGA database (https://portal.gdc.cancer.gov/). The dataset included survival data and FPKM expression value from 644 CRC tumor samples (478 COAD+166 READ). Normalized array data (GSE39582, GPL570 Affymetrix Human Genome U133 Plus 2.0 Array, France) and sample annotation files for 579 primary CRC tissues were obtained from Gene Expression Omnibus (GEO).

Identification of m6A-Related Gene Set

Currently known m6A RNA methylation regulators include methyltransferase (METTL3, METTL14, METTL16, WTAP, KIAA1429, RBM15, and ZC3H13), binding protein (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, and HNRNPC), and demethylase (ALKBH5 and FTO). There wer 3412 m6A-related genes correlated with CRC identified from the m6Avar database (http://rmvar.renlab.org/) [16]. The candidate m6A-related gene set was obtained by removing duplicate genes and the genes with no expression value or expression value less than 80% of total expression value in the samples. After refinement, 3161 candidate genes were kept for further analysis.

m6A-Related DEGs from Patients with Different Pathological Stages

To investigate the expression difference of m6A-related genes in different pathological stages, one-way ANOVA was used for identifying differentially expressed genes (DEGs) among 4 pathological stages. After removing samples without specific pathological stages, 620 tumor samples of CRC patients were left for further evaluation. Subsequently, all the differential expression analysis was executed with the threshold of P value <0.05. The heatmap which showed the expression differences of m6A-related genes in 4 different pathological stages was plotted using R package pheatmap.

Survival Analysis

Univariate Cox regression analysis was employed to identify m6A-related genes associated with the prognosis of CRC patients. According to the expression of genes, we classified the patients into high or low expression groups. With the survival time and survival status (live or dead) as input files, the overall survival (OS) probability between high or low expression groups was calculated and compared. OS of every CRC tumor sample was estimated by Kaplan-Meier method. The log-rank test was used to determine the significance of OS probability between different subgroups.

After that, the OS related genes were used for LASSO-Cox regression analysis to construct the prognostic model of these m6A-related genes with R package glmnet [23-25]. Ten-fold cross-validation minimum criteria were used to select the least-squares minimum value (min) with the minimum mean across validation error. Finally, we used to prognostic model to screen the gene set and calculate each patient’s risk score by a standard formula, which combines the expression levels of m6A-related genes with LASSO-Cox regression coefficients. To confirm whether the risk score independently affected the patients’ OS, the multivariate Cox regression model was used to assessed the association of pathoclinical features with the OS.

Construction of PPI Network

The search tool for the retrieval of interacting genes (STRING database, V11; http://string-db.org/) was employed to predict the protein-protein interactions network of prognostic m6A-related DEGs [17]. On the STRING website, after put prognostic m6A-related DEGs, the website will show PPI network files according to its internal database. Subsequently, Cytoscape software (http://cytoscape.org/) was applied to visualize and analyze biological networks and node degrees of the 146 candidate genes based on a confidence score >0.4 [18].

Gene Ontology and KEGG Pathway Enrichment Analysis

Gene ontology (GO) is a tool for gene annotation using a dynamic, controlled vocabulary that classifies genes into 3 categories: biological process, molecular function, and cellular component [19,20]. GO analysis uses different genes to annotate gene functions based on the GO database. After obtaining all the functions involved in the genes, Fisher’s exact test and multiple comparison test are used to calculate the significance level (P value) and false positive rate (FDR) of each function, to screen out the significant functions embodied by the differentially expressed genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) database is used to assign gene sets to specific pathway maps of molecular interactions, reactions, and relation networks [21]. At present, KEGG Pathway is divided into 8 categories: overall network, metabolic processes, genetic information transmission, environmental information transmission, intracellular biological processes, biological systems, human diseases, and drug development. Pathway analysis was based on the KEGG database, using Fisher’s exact test and chi-square test for differentially expressed genes to analyze the significance of the pathway participated by the target gene.

GO functional annotation and KEGG pathway enrichment analyses of DEGs were performed by R package clusterProfiler [22] with P value <0.05 as statistically significant to further explore the functions and involved pathways of differentially expressed m6A-related genes.

Statistical Analysis

All statistical analyses were conducted in R program (version 4.0.3; https://www.r-project.org/). Survival analysis was performed using Kaplan-Meier method and log-rank test with R package “survival” [26], while survival curves were plotted by “SurvMiner” package [27]. Cox proportional risk regression model was used for multivariate analysis. For all statistical tests, a P value of 0.05 was considered significant.

Results

Identification of Differentially Expressed m6A-Related Genes from Patients with Different Pathological Stages

We selected 15 m6A RNA methylation regulators and 3412 m6A RNA methylation-related genes. After removing duplicate genes and the genes with no expression value or expression values less than 80% of total expression value in the samples, 3161 candidate genes were kept for further analysis. The workflow was shown in Figure 1. Subsequently, a total of 673 differentially expressed genes (DEGs) were identified by one-way ANOVA analysis among different stages (Figure 2). Notably, 4 m6A RNA methylation regulatory factors, including methyltransferase like 14 (METTL14), YTH domain containing 2 (YTHDC2), YTH N6-methyladenosine RNA binding protein 2 (YTHDF2), and zinc finger CCCH-type containing 13 (ZC3H13), were significantly differentially expressed among the 4 pathological subgroups. The expression profile of METTL14, YTHDC2, and YTHDF2 exhibited sustained decreasing with the progression of CRC, while the expression level of ZC3H13 was steadily elevated during CRC development (Supplementary Figure 1).

Figure 1

The workflow of this study.

Figure 2

Heatmaps of 673 m6A-related DEGs from patients at different pathological stages. It was plotted using R package pheatmap (V1.0.12).

Screening of Prognostic m6A-Related Genes

To obtain prognostic m6A-related genes, univariate Cox regression analysis was performed with the threshold of P value <0.05. As a result, a total of 146 m6A-related genes were found to be associated with the overall survival of CRC patients (Supplementary Table 1).

Protein–Protein Interaction Network of the m6A-Related DEGs

The PPI networks for m6A-related genes were constructed according to the STRING database. The NetworkAnalyzer module was used to analyze and compare the network, hub genes were identified by analyzing the network topology. The top 10 hub genes were catenin beta 1 (CTNNB1), tripartite motif containing 37 (TRIM37), member RAS oncogene family (RAB7A), kinetochore scaffold 1 (CASC5/KNL1), centromere protein E (CENPE), G2/mitotic-specific cyclin-B1 (CCNB1), ubiquitin conjugating enzyme E2 H (UBE2H), heat shock protein family A (Hsp70) member 8 (HSPA8), kinesin family member 1A (KIF1A), and F-box and WD repeat domain containing 4 (FBXW4) (Figure 3).

Figure 3

PPI network of prognostic m6A-related genes constructed using STRING database (V11). In the diagram, genes are represented by nodes and their interactions are linked by lines. Genes with red color and large circle had higher degree values in the network, while genes with yellow color and small circle had lower degrees in the network. Cytoscape software (V3.5) was applied to visualize and analyze biological networks and node degrees of the 146 candidate genes.

Functional Annotation of m6A-Related Genes

To further explore the biological functions and involved pathways of prognostic m6A-related genes, we used R package clusterProfiler to perform GO function annotation and KEGG pathway enrichment analysis. The most significant (P value <0.05) GO terms and KEGG pathways are shown in Figure 4. Results indicated that cerebral cortex neuron differentiation, forebrain development, digestive tract development, digestive system development, and odontogenesis of dentin-containing teeth were the top 5 significantly enriched GO terms in biological process category, while pathways of neurodegeneration-multiple diseases, small cell lung cancer, Amoebiasis, Tuberculosis and HIF-1 signaling pathway were the top 5 significantly enriched KEGG pathways.

Figure 4

(A, B) Gene Ontology functional annotation and KEGG enrichment of m6A-related DEGs. GO functional annotation and KEGG pathway enrichment analyses of DEGs were performed by R package clusterProfiler (V3.14.3).

Construction of the Prognostic Risk Model

After 1000 resamples by Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression analysis through “glmnet” and “survival” R package, 9 m6A-related genes were selected to construct the prognostic model (Figure 5). Based on the gene signature with LASSO coefficients, the risk score=(0.02327*LRRC17)+ (0.8170*NFKB1)+ (0.07971*NOS2)+ (0.1193* PCDHB2)+ (0.006866*RAB7A)+ (0.07122*RPS6KA1)+ (−0.0897*RRNAD1)+ (0.3128*TLE6)+ (−0.06165*UBE2H). The TCGA cohort was then split into the high- and low-risk groups according to the median value of prognostic risk score. Kaplan-Meier survival analysis showed that patients with lower risk scores had significantly better overall survival than those with higher risk scores (P<0.0001) (Figure 5A). As shown in Figure 5B, the expression of RRNAD1 and UBE2H was negatively associated with the risk score of colorectal cancer patients, while the expression of other genes was positively associated with the risk score of patients with colorectal cancer. We also visualized the risk score distribution according to the length of follow-up months (Figure 5C). In the TCGA datasets, a significant positive correlation was found between the expression of NOS2 and TLE6 (Figure 5D). Figure 6A and 6B showed screening process of 9 m6A-related prognostic genes by LASSO-Cox regression analysis and random permutation. In Figure 6A, each line represents one gene, and the gene with non-zero coefficients was kept. In Figure 6B, when the line trends to flat, the lamda value would be chosen to conduct regression model. To verify whether these candidate prognostic gene models were influenced by clinical factors such as age, TNM stage and sex, we performed multivariate Cox regression analysis. The results showed that TNM stage, the status of vascular invasion, and the risk score were all independently related to the OS of TCGA CRC patients (Figure 6C).

Figure 5

The construction of m6A-related prognostic risk model. (A) The Kaplan-Meier survival curve describes a significant survival difference between the high-risk and low-risk groups in the prognostic model. (B) The coefficient of each selected maker. (C) The risk score curve of CRC patients and the survival status and survival time distribution according to the risk score. (D) Correlations among the 9 maker genes.

Figure 6

Lasso-Cox regression and multivariate Cox regression analysis of prognostic m6A-related genes. (A, B) 9 m6A-related prognostic genes were screened by LASSO-Cox regression analysis and random permutation. (C, D) The hazard ration (HR) and P values of the training set (TCGA) and the validation set (GSE39582) were calculated by multivariate Cox regression. LASSO-Cox regression analysis to construct the prognostic model of these m6A-related genes with R package glmnet (V4.1).

Validation of the Prognostic Model

Concurrently, the prognostic risk score of the m6A-related signature was also validated in the GSE39582 dataset. In the exploration, the median risk score of 573 patients in GSE39582 was calculated as 2.188 according to the risk score formula, among which 286 patients with a score <2.188 were assigned as the low-risk group, while other 287 patients with a score ≥2.188 specified as the high-risk group. Patients with lower risk score had significantly longer overall survival than those with higher risk score in the validation set (P value=0.0082, Figure 7). Multivariate Cox regression analysis also showed that risk score was a prognostic factor for CRC patients’ OS (Figure 6D).

Figure 7

The validation of constructed prognostic model using GEO dataset. (A) The Kaplan-Meier survival curve describes a significant survival difference between the high-risk and low-risk groups in the prognostic model. (B) The risk score curve of CRC patients in GSE39582 and the survival status and survival time distribution according to the risk score.

Discussion

Colorectal cancer is one of the most lethal solid tumors, with complex molecular and cellular heterogeneity. In the past few decades, there has been a great deal of research focusing on the molecular mechanisms of colorectal cancer, but most of them have concentrated on the aberration of protein-coding genes, leaving post-transcriptional processes mysterious. However, post-transcriptional alterations play a significant role in the preservation of tumor cells by modulating every hallmark in cancer [28,29]. RNA methylation modifications compose over 60% of all RNA modifications, and the most common type of RNA methylation modification is N6-methyladenosine (m6A) RNA methylation [30]. Accumulating evidence suggests that m6A modification plays a critical role not only in hypertension and cardiovascular disease, but also in tumor genesis and metastasis [31]. Therefore, the identification of m6A-related genes and m6A RNA methylation regulators abnormal expression may improve our understanding of colorectal cancer and provide us with valuable therapeutic targets.

Previous study about CRC and m6A methylation focused on m6A regulators [32], one subtype of CRC, such as COAD [33], or the comparison of tumor and normal tissues [34]. In the study, we systematically analyzed RNA sequencing data of CRC (including COAD and READ) patients from TCGA and GEO database. The m6A RNA methylation-related genes and m6A regulators were both brought into this investigation. We found that 673 candidate m6A RNA methylation-related genes were abnormally expressed among different pathological stages. Moreover, based on TCGA dataset, 146 m6A-related genes were related to patients’ overall survival based on the result of univariate Cox regression analysis. By constructing the PPI network, we identified 10 hub genes, namely CTNNB1, TRIM37, RAB7A, CASC5/KNL1, CENPE, CCNB1, UBE2H, HSPA8, KIF1A and FBXW4. CTNNB1 is a key component of the Wnt/β-catenin signaling pathway and was identified as the regulator of m6A modification in hepatoblastoma induced by METTL3 [35]. Liu et al found that Sec62 is upregulated by the METTL3-mediated m6A modification and promotes the stemness and chemoresistance of CRC by binding to β-catenin and enhancing Wnt signaling [36]. That might indicate that the WNT/β-catenin pathway should receive more attention in further study. A case-control study of breast cancer demonstrated that KIF1A promoter methylation can distinguish breast cancer (BC) cases from controls in plasma and was inversely associated with DNA repair ability (DRC) levels [37]. Studies confirmed that CCNB1 silencing can activate the p53 signaling pathway, further inhibit cell proliferation, and promote cell senescence in pancreatic cancer [38]. Meanwhile, in a systematic analysis of melanoma, it was proved to be positively correlated with either YTHDF1 or HNRNPA2B1, suggesting that both genes may affect m6A modification by CCNB1 gene [39].

To further investigate the influence of m6A RNA methylation regulatory factors on the prognosis of colorectal cancer, we used LASSO-Cox regression to establish a prognostic risk model based on 9 m6A RNA methylation-related genes. Survival analysis showed that the high-risk and low-risk subgroups classified by the model did have different prognostic destination in both the TCGA training group and the validation set. Although there has been no earlier research to prove that the biomarkers, we found in the prognostic model a close relationship with RNA modification process, but some of them (GSR and S1PR3) were correlated with the apoptosis and the adaption of acidic microenvironment of CRC cells [40,41]. Further experimental validation of the relationship between those biomarker genes and m6A RNA methylation regulators is needed to test the feasibility of our prognostic model.

Conclusions

In summary, our study systematically analyzed the expression profile of m6A RNA methylation-related genes, their prognostic significance, potential functions and pathways, and protein-protein interactions from CRC patients with the help of TCGA and GEO databases. Nine candidate m6A-related mRNA biomarkers (LRRC17, NFKB1, NOS2, PCDHB2, RAB7A, RPS6KA1, RRNAD1, TLE6, and UBE2H) were found to be closely related to the clinicopathology and prognosis of colorectal cancer. This study not only suggests the potential value of m6A-related genes as novel prognostic biomarkers for colorectal cancer, but also provides important clues for the diagnosis and treatment of colorectal cancer patients.

The expression of m6A RNA methylation regulatory factor genes from patients at different pathological stages.

Supplementary Table 1

Detailed information (including univariates Cox regression P value and HR score) of 146 m6A-related genes that were associated with the overall survival of CRC patients.

Ensembl ID	Cox p value	HR	chr	Source	Type	Start	End	Strand	Gene type	Gene name
ENSG00000058404	0	3	chr7	HAVANA	Gene	44217150	44334577	−	Protein coding	CAMK2B
ENSG00000069696	0	1.4	chr11	HAVANA	Gene	637293	640706	+	Protein coding	DRD4
ENSG00000072832	0	1.2	chr4	HAVANA	Gene	5748084	5893058	−	Protein coding	CRMP1
ENSG00000077454	0	1.17	chr7	HAVANA	Gene	1.01E+08	1.01E+08	−	Protein coding	LRCH4
ENSG00000078018	0	1.69	chr2	HAVANA	Gene	2.09E+08	2.1E+08	+	Protein coding	MAP2
ENSG00000130294	0	1.21	chr2	HAVANA	Gene	2.41E+08	2.41E+08	−	Protein coding	KIF1A
ENSG00000133069	0	2.31	chr1	HAVANA	Gene	2.05E+08	2.05E+08	+	Protein coding	TMCC2
ENSG00000142549	0	1.45	chr19	HAVANA	Gene	51311848	51330354	+	Protein coding	IGLON5
ENSG00000160117	0	1.64	chr19	HAVANA	Gene	17281645	17287646	+	Protein coding	ANKLE1
ENSG00000167311	0	1.86	chr11	HAVANA	Gene	3638503	3642316	−	Protein coding	ART5
ENSG00000169169	0	1.65	chr19	HAVANA	Gene	49690898	49713731	+	Protein coding	CPT1C
ENSG00000187730	0	1.27	chr1	HAVANA	Gene	2019298	2030758	+	Protein coding	GABRD
ENSG00000104687	0.001	0.98	chr8	HAVANA	Gene	30678061	30727999	−	Protein coding	GSR
ENSG00000115694	0.001	1.08	chr2	HAVANA	Gene	2.41E+08	2.42E+08	−	Protein coding	STK25
ENSG00000148948	0.001	4.76	chr11	HAVANA	Gene	40114203	41459773	−	Protein coding	LRRC4C
ENSG00000151640	0.001	1.59	chr10	HAVANA	Gene	1.32E+08	1.32E+08	+	Protein coding	DPYSL4
ENSG00000168994	0.001	1.09	chr6	HAVANA	Gene	3722614	3752026	−	Protein coding	PXDC1
ENSG00000250510	0.001	1.97	chr12	HAVANA	Gene	6821545	6829972	+	Protein coding	GPR162
ENSG00000076344	0.002	1.23	chr16	HAVANA	Gene	268301	275980	−	Protein coding	RGS11
ENSG00000109606	0.002	0.95	chr4	HAVANA	Gene	24517441	24584550	−	Protein coding	DHX15
ENSG00000112852	0.002	1.4	chr5	HAVANA	Gene	1.41E+08	1.41E+08	+	Protein coding	PCDHB2
ENSG00000130227	0.002	0.93	chr8	HAVANA	Gene	21919671	22006585	+	Protein coding	XPO7
ENSG00000143184	0.002	1.34	chr1	HAVANA	Gene	1.69E+08	1.69E+08	+	Protein coding	XCL1
ENSG00000143303	0.002	1.08	chr1	HAVANA	Gene	1.57E+08	1.57E+08	+	Protein coding	RRNAD1
ENSG00000156006	0.002	0.9	chr8	HAVANA	Gene	18391245	18401218	+	Protein coding	NAT2
ENSG00000172340	0.002	0.98	chr3	HAVANA	Gene	67360460	67654614	−	Protein coding	SUCLG2
ENSG00000178209	0.002	1.01	chr8	HAVANA	Gene	1.44E+08	1.44E+08	−	Protein coding	PLEC
ENSG00000103168	0.003	1.09	chr16	HAVANA	Gene	84177847	84187070	−	Protein coding	TAF1C
ENSG00000121940	0.003	0.8	chr1	HAVANA	Gene	1.09E+08	1.09E+08	−	Protein coding	CLCC1
ENSG00000204947	0.003	2.03	chr7	HAVANA	Gene	1.49E+08	1.49E+08	−	Protein coding	ZNF425
ENSG00000047188	0.004	0.78	chr5	HAVANA	Gene	1.14E+08	1.14E+08	+	Protein coding	YTHDC2
ENSG00000112333	0.004	100.3	chr6	HAVANA	Gene	1.08E+08	1.08E+08	+	Protein coding	NR2E1
ENSG00000149289	0.004	0.71	chr11	HAVANA	Gene	1.1E+08	1.1E+08	+	Protein coding	ZC3H12C
ENSG00000151611	0.004	0.61	chr4	HAVANA	Gene	1.46E+08	1.46E+08	+	Protein coding	MMAA
ENSG00000164307	0.004	0.95	chr5	HAVANA	Gene	96760810	96808100	−	Protein coding	ERAP1
ENSG00000254585	0.004	3.44	chr15	HAVANA	Gene	23643544	23647841	−	Protein coding	MAGEL2
ENSG00000257923	0.004	1.06	chr7	HAVANA	Gene	1.02E+08	1.02E+08	+	Protein coding	CUX1
ENSG00000105737	0.005	1.53	chr19	HAVANA	Gene	41998321	42069498	−	Protein coding	GRIK5
ENSG00000146830	0.005	1.04	chr7	HAVANA	Gene	1.01E+08	1.01E+08	−	Protein coding	GIGYF1
ENSG00000172137	0.005	1.04	chr16	HAVANA	Gene	71358713	71390438	+	Protein coding	CALB2
ENSG00000186591	0.005	1.04	chr7	HAVANA	Gene	1.3E+08	1.3E+08	−	Protein coding	UBE2H
ENSG00000188549	0.005	1.15	chr15	HAVANA	Gene	40331452	40340967	−	Protein coding	C15orf52
ENSG00000197217	0.005	0.91	chr8	HAVANA	Gene	23385783	23457695	−	Protein coding	ENTPD4
ENSG00000104953	0.006	1.73	chr19	HAVANA	Gene	2977446	2995184	+	Protein coding	TLE6
ENSG00000143434	0.006	1.48	chr1	HAVANA	Gene	1.51E+08	1.51E+08	−	Protein coding	SEMA6C
ENSG00000116254	0.007	3.64	chr1	HAVANA	Gene	6101793	6180123	−	Protein coding	CHD5
ENSG00000124074	0.007	1.09	chr16	HAVANA	Gene	67662945	67667265	−	Protein coding	ENKD1
ENSG00000129244	0.007	1.6	chr17	HAVANA	Gene	7646627	7657768	+	Protein coding	ATP1B2
ENSG00000141756	0.007	1.01	chr17	HAVANA	Gene	41812680	41823217	+	Protein coding	FKBP10
ENSG00000168036	0.007	0.99	chr3	HAVANA	Gene	41194741	41260096	+	Protein coding	CTNNB1
ENSG00000233251	0.007	1.79	chr2	HAVANA	Gene	56173534	56185770	−	Antisense	AC007743.1
ENSG00000155975	0.008	0.85	chr8	HAVANA	Gene	17246571	17302427	+	Protein coding	VPS37A
ENSG00000169436	0.008	1.21	chr8	HAVANA	Gene	1.39E+08	1.39E+08	−	Protein coding	COL22A1
ENSG00000173137	0.008	1.04	chr8	HAVANA	Gene	1.44E+08	1.44E+08	+	Protein coding	ADCK5
ENSG00000188559	0.008	0.95	chr20	HAVANA	Gene	20389552	20712488	−	Protein coding	RALGAPA2
ENSG00000205090	0.008	1.83	chr1	HAVANA	Gene	1535174	1540453	−	Protein coding	TMEM240
ENSG00000061273	0.009	1.1	chr12	HAVANA	Gene	47782722	47833132	−	Protein coding	HDAC7
ENSG00000068784	0.009	0.83	chr2	HAVANA	Gene	45388680	45612165	−	Protein coding	SRBD1
ENSG00000128268	0.009	1.08	chr22	HAVANA	Gene	39457344	39492194	+	Protein coding	MGAT3
ENSG00000154359	0.009	0.82	chr8	HAVANA	Gene	12721894	12756073	−	Protein coding	LONRF1
ENSG00000163064	0.009	3.33	chr2	HAVANA	Gene	1.19E+08	1.19E+08	−	Protein coding	EN1
ENSG00000171246	0.009	1.29	chr17	HAVANA	Gene	80467148	80477843	−	Protein coding	NPTX1
ENSG00000213694	0.009	1.15	chr9	HAVANA	Gene	88991447	89005010	+	Protein coding	S1PR3
ENSG00000038210	0.01	0.94	chr4	HAVANA	Gene	25233975	25279092	+	Protein coding	PI4K2B
ENSG00000066557	0.01	0.88	chr1	HAVANA	Gene	70144805	70205620	−	Protein coding	LRRC40
ENSG00000075785	0.01	1.01	chr3	HAVANA	Gene	1.29E+08	1.29E+08	+	Protein coding	RAB7A
ENSG00000115275	0.01	1.03	chr2	HAVANA	Gene	74461057	74465410	−	Protein coding	MOGS
ENSG00000117155	0.01	0.86	chr1	HAVANA	Gene	84643707	84690803	−	Protein coding	SSX2IP
ENSG00000120868	0.01	0.86	chr12	HAVANA	Gene	98645141	98735433	+	Protein coding	APAF1
ENSG00000138604	0.01	0.95	chr15	HAVANA	Gene	69160584	69272217	+	Protein coding	GLCE
ENSG00000175505	0.01	1.08	chr11	HAVANA	Gene	67364168	67374177	−	Protein coding	CLCF1
ENSG00000033100	0.011	1.04	chr7	HAVANA	Gene	1.51E+08	1.51E+08	+	Protein coding	CHPF2
ENSG00000126903	0.011	1.02	chrX	HAVANA	Gene	1.54E+08	1.54E+08	−	Protein coding	SLC10A3
ENSG00000177692	0.011	0.19	chr21	HAVANA	Gene	33485530	33491720	−	Protein coding	DNAJC28
ENSG00000018510	0.012	0.93	chr2	HAVANA	Gene	1.77E+08	1.78E+08	+	Protein coding	AGPS
ENSG00000151623	0.012	0.89	chr4	HAVANA	Gene	1.48E+08	1.48E+08	−	Protein coding	NR3C2
ENSG00000196700	0.012	1.28	chr20	HAVANA	Gene	63956702	63969865	−	Protein coding	ZNF512B
ENSG00000082684	0.013	2.77	chr3	HAVANA	Gene	1.23E+08	1.23E+08	−	Protein coding	SEMA5B
ENSG00000106290	0.013	1.08	chr7	HAVANA	Gene	1E+08	1E+08	−	Protein coding	TAF6
ENSG00000164142	0.013	0.67	chr4	HAVANA	Gene	1.51E+08	1.52E+08	+	Protein coding	FAM160A1
ENSG00000167384	0.013	0.69	chr19	HAVANA	Gene	44474428	44500524	−	Protein coding	ZNF180
ENSG00000109971	0.014	1	chr11	HAVANA	Gene	1.23E+08	1.23E+08	−	Protein coding	HSPA8
ENSG00000148291	0.014	1.02	chr9	HAVANA	Gene	1.33E+08	1.33E+08	+	Protein coding	SURF2
ENSG00000064651	0.015	0.99	chr5	HAVANA	Gene	1.28E+08	1.28E+08	+	Protein coding	SLC12A2
ENSG00000149925	0.015	1	chr16	HAVANA	Gene	30064164	30070457	+	Protein coding	ALDOA
ENSG00000007171	0.016	0.98	chr17	HAVANA	Gene	27756766	27800499	−	Protein coding	NOS2
ENSG00000014164	0.016	1.04	chr8	HAVANA	Gene	1.43E+08	1.44E+08	−	Protein coding	ZC3H3
ENSG00000109320	0.016	0.94	chr4	HAVANA	Gene	1.03E+08	1.03E+08	+	Protein coding	NFKB1
ENSG00000077254	0.018	0.92	chr1	HAVANA	Gene	77695987	77759852	−	Protein coding	USP33
ENSG00000121361	0.018	1.12	chr12	HAVANA	Gene	21764955	21775581	−	Protein coding	KCNJ8
ENSG00000128606	0.018	1.23	chr7	HAVANA	Gene	1.03E+08	1.03E+08	+	Protein coding	LRRC17
ENSG00000160293	0.018	1.03	chr9	HAVANA	Gene	1.34E+08	1.34E+08	−	Protein coding	VAV2
ENSG00000099260	0.019	1.31	chr1	HAVANA	Gene	99645943	99694541	+	Protein coding	PALMD
ENSG00000138131	0.019	1.32	chr10	HAVANA	Gene	98247690	98268250	−	Protein coding	LOXL4
ENSG00000167994	0.02	1.09	chr11	HAVANA	Gene	61897301	61920269	−	Protein coding	RAB3IL1
ENSG00000103966	0.021	0.95	chr15	HAVANA	Gene	41895939	41972578	−	Protein coding	EHD4
ENSG00000160445	0.021	1.06	chr9	HAVANA	Gene	1.29E+08	1.29E+08	−	Protein coding	ZER1
ENSG00000162148	0.022	1.48	chr11	HAVANA	Gene	61481120	61490931	+	Protein coding	PPP1R32
ENSG00000117676	0.023	0.97	chr1	HAVANA	Gene	26529761	26575030	+	Protein coding	RPS6KA1
ENSG00000187244	0.023	1.02	chr19	HAVANA	Gene	44809059	44821421	+	Protein coding	BCAM
ENSG00000109805	0.024	0.91	chr4	HAVANA	Gene	17810902	17844862	+	Protein coding	NCAPG
ENSG00000158106	0.024	1.03	chr8	HAVANA	Gene	1.43E+08	1.43E+08	+	Protein coding	RHPN1
ENSG00000178409	0.024	0.76	chr6	HAVANA	Gene	1.07E+08	1.07E+08	−	Protein coding	BEND3
ENSG00000134057	0.025	0.99	chr5	HAVANA	Gene	69167010	69178245	+	Protein coding	CCNB1
ENSG00000187775	0.025	2.49	chr17	HAVANA	Gene	78423697	78577394	−	Protein coding	DNAH17
ENSG00000204540	0.025	1.2	chr6	HAVANA	Gene	31114750	31140092	+	Protein coding	PSORS1C1
ENSG00000162600	0.026	0.9	chr1	HAVANA	Gene	58415384	58546802	−	Protein coding	OMA1
ENSG00000163132	0.027	1.03	chr4	HAVANA	Gene	4859666	4863936	+	Protein coding	MSX1
ENSG00000137812	0.028	0.85	chr15	HAVANA	Gene	40594020	40664342	+	Protein coding	CASC5
ENSG00000105321	0.029	1.08	chr19	HAVANA	Gene	47255980	47273701	+	Protein coding	CCDC9
ENSG00000130702	0.029	1.02	chr20	HAVANA	Gene	62307955	62367312	−	Protein coding	LAMA5
ENSG00000198026	0.029	1.1	chr20	HAVANA	Gene	45948653	45972172	−	Protein coding	ZNF335
ENSG00000039523	0.031	1.04	chr16	HAVANA	Gene	67518418	67546788	+	Protein coding	FAM65A
ENSG00000090686	0.031	0.89	chr1	HAVANA	Gene	21678298	21783606	−	Protein coding	USP48
ENSG00000189184	0.031	0.88	chr4	HAVANA	Gene	1.38E+08	1.38E+08	−	Protein coding	PCDH18
ENSG00000112367	0.032	0.88	chr6	HAVANA	Gene	1.1E+08	1.1E+08	+	Protein coding	FIG4
ENSG00000143845	0.032	1.22	chr1	HAVANA	Gene	2.04E+08	2.04E+08	−	Protein coding	ETNK2
ENSG00000130158	0.033	1.06	chr19	HAVANA	Gene	11199295	11262481	−	Protein coding	DOCK6
ENSG00000159314	0.033	1.06	chr17	HAVANA	Gene	45393902	45434421	−	Protein coding	ARHGAP27
ENSG00000107829	0.035	1.06	chr10	HAVANA	Gene	1.02E+08	1.02E+08	−	Protein coding	FBXW4
ENSG00000128482	0.035	1.62	chr17	HAVANA	Gene	19411125	19417276	+	Protein coding	RNF112
ENSG00000164181	0.035	0.94	chr5	HAVANA	Gene	60751791	60844389	−	Protein coding	ELOVL7
ENSG00000127838	0.036	1.02	chr2	HAVANA	Gene	2.18E+08	2.18E+08	+	Protein coding	PNKD
ENSG00000180287	0.036	702.2	chr1	HAVANA	Gene	2.42E+08	2.43E+08	−	Protein coding	PLD5
ENSG00000132563	0.037	1.22	chr5	HAVANA	Gene	1.38E+08	1.38E+08	+	Protein coding	REEP2
ENSG00000155760	0.037	1.03	chr2	HAVANA	Gene	2.02E+08	2.02E+08	+	Protein coding	FZD7
ENSG00000183323	0.039	0.83	chr5	HAVANA	Gene	69280175	69332809	−	Protein coding	CCDC125
ENSG00000070031	0.04	1.01	chr11	HAVANA	Gene	626431	627143	−	Protein coding	SCT
ENSG00000115616	0.041	0.94	chr2	HAVANA	Gene	1.03E+08	1.03E+08	+	Protein coding	SLC9A2
ENSG00000134443	0.041	1.07	chr18	HAVANA	Gene	59220168	59230774	+	Protein coding	GRP
ENSG00000143412	0.041	1.06	chr1	HAVANA	Gene	1.51E+08	1.51E+08	+	Protein coding	ANXA9
ENSG00000108395	0.042	0.9	chr17	HAVANA	Gene	58982638	59106921	−	Protein coding	TRIM37
ENSG00000135406	0.042	1.35	chr12	HAVANA	Gene	49293252	49298686	+	Protein coding	PRPH
ENSG00000155846	0.043	0.75	chr5	HAVANA	Gene	1.5E+08	1.5E+08	+	Protein coding	PPARGC1B
ENSG00000008710	0.044	1.1	chr16	HAVANA	Gene	2088710	2135898	−	Protein coding	PKD1
ENSG00000119487	0.044	1.06	chr9	HAVANA	Gene	1.25E+08	1.26E+08	−	Protein coding	MAPKAP1
ENSG00000166033	0.044	1.01	chr10	HAVANA	Gene	1.22E+08	1.23E+08	+	Protein coding	HTRA1
ENSG00000091164	0.045	0.9	chr18	HAVANA	Gene	56597208	56651600	−	Protein coding	TXNL1
ENSG00000067445	0.046	1.51	chrX	HAVANA	Gene	54920462	54931431	+	Protein coding	TRO
ENSG00000081923	0.046	0.98	chr18	HAVANA	Gene	57646426	57803101	−	Protein coding	ATP8B1
ENSG00000138686	0.046	0.78	chr4	HAVANA	Gene	1.22E+08	1.22E+08	−	Protein coding	BBS7
ENSG00000138778	0.046	0.85	chr4	HAVANA	Gene	1.03E+08	1.03E+08	−	Protein coding	CENPE
ENSG00000142235	0.046	1.18	chr19	HAVANA	Gene	48485271	48513189	−	Protein coding	LMTK3
ENSG00000136052	0.048	0.91	chr12	HAVANA	Gene	1.05E+08	1.05E+08	−	Protein coding	SLC41A2
ENSG00000160714	0.048	1.03	chr1	HAVANA	Gene	1.55E+08	1.55E+08	−	Protein coding	UBE2Q1
ENSG00000179041	0.049	1.01	chr8	HAVANA	Gene	66429028	66430733	+	Protein coding	RRS1