Extensive increase of microarray signals in cancers calls for novel normalization assumptions

A unique circulating microRNA pairs signature serves as a superior tool for early diagnosis of pan-cancer

Cancer remains a major burden globally and the critical role of early diagnosis is self-evident. Although various miRNA-based signatures have been developed in past decades, clinical utilization is limited due to a lack of precise cutoff value. Here, we innovatively developed a signature based on pairwise expression of miRNAs (miRPs) for pan-cancer diagnosis using machine learning approach. We analyzed miRNA spectrum of 15832 patients, who were divided into training, validation, test, and external test sets, with 13 different cancers from 10 cohorts. Five different machine-learning (ML) algorithms (XGBoost, SVM, RandomForest, LASSO, and Logistic) were adopted for signature construction. The best ML algorithm and the optimal number of miRPs included were identified using area under the curve (AUC) and youden index in validation set. The AUC of the best model was compared to previously published 25 signatures. Overall, Random Forest approach including 31 miRPs (31-miRP) was developed, proving highly efficient in cancer diagnosis across different datasets and cancer types (AUC range: 0.980-1.000). Regarding diagnosis of cancers at early stage, 31-miRP also exhibited high capacities, with AUC ranging from 0.961 to 0.998. Moreover, 31-miRP exhibited advantages in differentiating cancers from normal tissues (AUC range: 0.976-0.998) as well as differentiating cancers from corresponding benign lesions. Encouragingly, comparing to previously published 25 different signatures, 31-miRP also demonstrated clear advantages. In conclusion, 31-miRP acts as a powerful model for cancer diagnosis, characterized by high specificity and sensitivity as well as a clear cutoff value, thereby holding potential as a reliable tool for cancer diagnosis at early stage.

Voluntary running exercise modifies astrocytic population and features in the peri-infarct cortex

Rehabilitative exercise following a brain stroke has beneficial effects on the morphological plasticity of neurons. Particularly, voluntary running exercise after focal cerebral ischemia promotes functional recovery and ameliorates ischemia-induced dendritic spine loss in the peri-infarct motor cortex layer 5. Moreover, neuronal morphology is affected by changes in the perineuronal environment. Glial cells, whose phenotypes may be altered by exercise, are known to play a pivotal role in the formation of this perineuronal environment. Herein, we investigated the effects of voluntary running exercise on glial cells after middle cerebral artery occlusion. Voluntary running exercise increased the population of glial fibrillary acidic protein-positive astrocytes born between post-operative days (POD) 0 and 3 on POD15 in the peri-infarct cortex. After exercise, transcriptomic analysis of post-ischemic astrocytes revealed 10 upregulated and 70 downregulated genes. Furthermore, gene ontology analysis showed that the 70 downregulated genes were significantly associated with neuronal morphology. In addition, exercise reduced the number of astrocytes expressing lipocalin 2, a regulator of dendritic spine density, on POD15. Our results suggest that exercise modifies the composition of astrocytic population and their phenotype.

Application of individualized differential expression analysis in human cancer proteome

Liquid chromatography–mass spectrometry-based quantitative proteomics can measure the expression of thousands of proteins from biological samples and has been increasingly applied in cancer research. Identifying differentially expressed proteins (DEPs) between tumors and normal controls is commonly used to investigate carcinogenesis mechanisms. While differential expression analysis (DEA) at an individual level is desired to identify patient-specific molecular defects for better patient stratification, most statistical DEP analysis methods only identify deregulated proteins at the population level. To date, robust individualized DEA algorithms have been proposed for ribonucleic acid data, but their performance on proteomics data is underexplored. Herein, we performed a systematic evaluation on five individualized DEA algorithms for proteins on cancer proteomic datasets from seven cancer types. Results show that the within-sample relative expression orderings (REOs) of protein pairs in normal tissues were highly stable, providing the basis for individualized DEA for proteins using REOs. Moreover, individualized DEA algorithms achieve higher precision in detecting sample-specific deregulated proteins than population-level methods. To facilitate the utilization of individualized DEA algorithms in proteomics for prognostic biomarker discovery and personalized medicine, we provide Individualized DEP Analysis IDEPAXMBD (XMBD: Xiamen Big Data, a biomedical open software initiative in the National Institute for Data Science in Health and Medicine, Xiamen University, China.) (https://github.com/xmuyulab/IDEPA-XMBD), which is a user-friendly and open-source Python toolkit that integrates individualized DEA algorithms for DEP-associated deregulation pattern recognition.

Comparison of Data Normalization Strategies for Array-Based MicroRNA Profiling Experiments and Identification and Validation of Circulating MicroRNAs as Endogenous Controls in Hypertension

Introduction: MicroRNAs are small noncoding RNAs with potential regulatory roles in hypertension and drug response. The presence of many of these RNAs in biofluids has spurred investigation into their role as possible biomarkers for use in precision approaches to healthcare. One of the major challenges in clinical translation of circulating miRNA biomarkers is the limited replication across studies due to lack of standards for data normalization techniques for array-based approaches and a lack of consensus on an endogenous control normalizer for qPCR-based candidate miRNA profiling studies.

Methods: We conducted genome-wide profiling of 754 miRNAs in baseline plasma of 36 European American individuals with uncomplicated hypertension selected from the PEAR clinical trial, who had been untreated for hypertension for at least one month prior to sample collection. After appropriate quality control with amplification score and missingness filters, we tested different normalization strategies such as normalization with global mean of imputed and unimputed data, mean of restricted set of miRNAs, quantile normalization, and endogenous control miRNA normalization to identify the method that best reduces the technical/experimental variability in the data. We identified best endogenous control candidates with expression pattern closest to the mean miRNA expression in the sample, as well as by assessing their stability using a combination of NormFinder, geNorm, Best Keeper and Delta Ct algorithms under the Reffinder software. The suitability of the four best endogenous controls was validated in 50 hypertensive African Americans from the same trial with reverse-transcription–qPCR and by evaluating their stability ranking in that cohort.

Results: Among the compared normalization strategies, quantile normalization and global mean normalization performed better than others in terms of reducing the standard deviation of miRNAs across samples in the array-based data. Among the four strongest candidate miRNAs from our selection process (miR-223-3p, 19b, 106a, and 126-5p), miR-223-3p and miR-126-5p were consistently expressed with the best stability ranking in the validation cohort. Furthermore, the combination of miR-223-3p and 126-5p showed better stability ranking when compared to single miRNAs.

Conclusion: We identified quantile normalization followed by global mean normalization to be the best methods in reducing the variance in the data. We identified the combination of miR-223-3p and 126-5p as potential endogenous control in studies of hypertension.

How to do quantile normalization correctly for gene expression data analyses

Quantile normalization is an important normalization technique commonly used in high-dimensional data analysis. However, it is susceptible to class-effect proportion effects (the proportion of class-correlated variables in a dataset) and batch effects (the presence of potentially confounding technical variation) when applied blindly on whole data sets, resulting in higher false-positive and false-negative rates. We evaluate five strategies for performing quantile normalization, and demonstrate that good performance in terms of batch-effect correction and statistical feature selection can be readily achieved by first splitting data by sample class-labels before performing quantile normalization independently on each split (“Class-specific”). Via simulations with both real and simulated batch effects, we demonstrate that the “Class-specific” strategy (and others relying on similar principles) readily outperform whole-data quantile normalization, and is robust-preserving useful signals even during the combined analysis of separately-normalized datasets. Quantile normalization is a commonly used procedure. But when carelessly applied on whole datasets without first considering class-effect proportion and batch effects, can result in poor performance. If quantile normalization must be used, then we recommend using the “Class-specific” strategy.

Avoid Oversimplifications in Machine Learning: Going beyond the Class-Prediction Accuracy

Class-prediction accuracy provides a quick but superficial way of determining classifier performance. It does not inform on the reproducibility of the findings or whether the selected or constructed features used are meaningful and specific. Furthermore, the class-prediction accuracy oversummarizes and does not inform on how training and learning have been accomplished: two classifiers providing the same performance in one validation can disagree on many future validations. It does not provide explainability in its decision-making process and is not objective, as its value is also affected by class proportions in the validation set. Despite these issues, this does not mean we should omit the class-prediction accuracy. Instead, it needs to be enriched with accompanying evidence and tests that supplement and contextualize the reported accuracy. This additional evidence serves as augmentations and can help us perform machine learning better while avoiding naive reliance on oversimplified metrics.

There is a huge potential for machine learning, but blind reliance on oversimplified metrics can mislead. Class-prediction accuracy is a common metric used for determining classifier performance. This article provides examples to show how the class-prediction accuracy is superficial and even misleading. We propose some augmentative measures to supplement the class-prediction accuracy. This in turn helps us to better understand the quality of learning of the classifier.

A qualitative transcriptional signature for the histological reclassification of lung squamous cell carcinomas and adenocarcinomas

Background

Targeted therapy for non-small cell lung cancer is histology dependent. However, histological classification by routine pathological assessment with hematoxylin-eosin staining and immunostaining for poorly differentiated tumors, particularly those from small biopsies, is still challenging. Additionally, the effectiveness of immunomarkers is limited by technical inconsistencies of immunostaining and lack of standardization for staining interpretation.

Results

Using gene expression profiles of pathologically-determined lung adenocarcinomas and squamous cell carcinomas, denoted as pADC and pSCC respectively, we developed a qualitative transcriptional signature, based on the within-sample relative gene expression orderings (REOs) of gene pairs, to distinguish ADC from SCC. The signature consists of two genes, KRT5 and AGR2, which has the stable REO pattern of KRT5 > AGR2 in pSCC and KRT5 < AGR2 in pADC. In the two test datasets with relative unambiguous NSCLC types, the apparent accuracy of the signature were 94.44 and 98.41%, respectively. In the other integrated dataset for frozen tissues, the signature reclassified 4.22% of the 805 pADC patients as SCC and 12% of the 125 pSCC patients as ADC. Similar results were observed in the clinical challenging cases, including FFPE specimens, mixed tumors, small biopsy specimens and poorly differentiated specimens. The survival analyses showed that the pADC patients reclassified as SCC had significantly shorter overall survival than the signature-confirmed pADC patients (log-rank p = 0.0123, HR = 1.89), consisting with the knowledge that SCC patients suffer poor prognoses than ADC patients. The proliferative activity, subtype-specific marker genes and consensus clustering analyses also supported the correctness of our signature.

Conclusions

The non-subjective qualitative REOs signature could effectively distinguish ADC from SCC, which would be an auxiliary test for the pathological assessment of the ambiguous cases.

Circumvent the uncertainty in the applications of transcriptional signatures to tumor tissues sampled from different tumor sites

The expression measurements of thousands of genes are correlated with the proportions of tumor epithelial cell (PTEC) in clinical samples. Thus, for a tumor diagnostic or prognostic signature based on a summarization of expression levels of the signature genes, the risk score for a patient may dependent on the tumor tissues sampled from different tumor sites with diverse PTEC for the same patient. Here, we proposed that the within-samples relative expression orderings (REOs) based gene pairs signatures should be insensitive to PTEC variations. Firstly, by analysis of paired tumor epithelial cell and stromal cell microdissected samples from 27 cancer patients, we showed that above 80% of gene pairs had consistent REOs between the two cells, indicating these REOs would be independent of PTEC in cancer tissues. Then, by simulating tumor tissues with different PTEC using each of the 27 paired samples, we showed that about 90% REOs of gene pairs in tumor epithelial cells were maintained in tumor samples even when PTEC decreased to 30%. Especially, the REOs of gene pairs with larger expression differences in tumor epithelial cells tend to be more robust against PTEC variations. Finally, as a case study, we developed a gene pair signature which could robustly distinguish colorectal cancer tissues with various PTEC from normal tissues. We concluded that the REOs-based signatures were robust against PTEC variations.

Robust transcriptional tumor signatures applicable to both formalin-fixed paraffin-embedded and fresh-frozen samples

Formalin-fixed paraffin-embedded (FFPE) samples represent a valuable resource for clinical researches. However, FFPE samples are usually considered an unreliable source for gene expression analysis due to the partial RNA degradation. In this study, through comparing gene expression profiles between FFPE samples and paired fresh-frozen (FF) samples for three cancer types, we firstly showed that expression measurements of thousands of genes had at least two-fold change in FFPE samples compared with paired FF samples. Therefore, for a transcriptional signature based on risk scores summarized from the expression levels of the signature genes, the risk score thresholds trained from FFPE (or FF) samples could not be applied to FF (or FFPE) samples. On the other hand, we found that more than 90% of the relative expression orderings (REOs) of gene pairs in the FF samples were maintained in their paired FFPE samples and largely unaffected by the storage time. The result suggested that the REOs of gene pairs were highly robust against partial RNA degradation in FFPE samples. Finally, as a case study, we developed a REOs-based signature to distinguish liver cirrhosis from hepatocellular carcinoma (HCC) using FFPE samples. The signature was validated in four datasets of FFPE samples and eight datasets of FF samples. In conclusion, the valuable FFPE samples can be fully exploited to identify REOs-based diagnostic and prognostic signatures which could be robustly applicable to both FF samples and FFPE samples with degraded RNA.

Discriminating cancer-related and cancer-unrelated chemoradiation-response genes for locally advanced rectal cancers

For patients with locally advanced rectal cancer (LARC) treated with preoperation chemoradiation (pCRT), identifying differentially expressed (DE) genes between non-responders and responders is a common approach for investigating mechanisms of chemoradiation resistance. However, some of such DE genes might be irrelevant to cancer itself but simply reflect the pharmacokinetic differences of the normal tissues. In this study, we adopted the RankComp algorithm to identify DE genes for each of LARC sample compared with its own normal state. Then, we identified genes with significantly different deregulation frequencies between the non-responders and responders, defined as cancer-related pCRT-response genes. Pathway enrichment and protein-protein interaction analyses showed that these genes specifically and intensively interacted with currently known effective genes of pCRT, involving in DNA replication, cell cycle and DNA repair. In contrast, after excluding the cancer-related pCRT-response genes, the other DE genes between non-responders and responders were enriched in many pathways of drug and protein metabolisms and transports, and interacted with both the known effective genes and pharmacokinetic genes. Hence, these two types of DE genes should be distinguished for investigating mechanisms of pCRT response in LARCs.

Differential expression analysis for individual cancer samples based on robust within-sample relative gene expression orderings across multiple profiling platforms

The highly stable within-sample relative expression orderings (REOs) of gene pairs in a particular type of human normal tissue are widely reversed in the cancer condition. Based on this finding, we have recently proposed an algorithm named RankComp to detect differentially expressed genes (DEGs) for individual disease samples measured by a particular platform. In this paper, with 461 normal lung tissue samples separately measured by four commonly used platforms, we demonstrated that tens of millions of gene pairs with significantly stable REOs in normal lung tissue can be consistently detected in samples measured by different platforms. However, about 20% of stable REOs commonly detected by two different platforms (e.g., Affymetrix and Illumina platforms) showed inconsistent REO patterns due to the differences in probe design principles. Based on the significantly stable REOs (FDR<0.01) for normal lung tissue consistently detected by the four platforms, which tended to have large rank differences, RankComp detected averagely 1184, 1335 and 1116 DEGs per sample with averagely 96.51%, 95.95% and 94.78% precisions in three evaluation datasets with 25, 57 and 58 paired lung cancer and normal samples, respectively. Individualized pathway analysis revealed some common and subtype-specific functional mechanisms of lung cancer. Similar results were observed for colorectal cancer. In conclusion, based on the cross-platform significantly stable REOs for a particular normal tissue, differentially expressed genes and pathways in any disease sample measured by any of the platforms can be readily and accurately detected, which could be further exploited for dissecting the heterogeneity of cancer.

Revealing potential molecular targets bridging colitis and colorectal cancer based on multidimensional integration strategy

Chronic inflammation may play a vital role in the pathogenesis of inflammation-associated tumors. However, the underlying mechanisms bridging ulcerative colitis (UC) and colorectal cancer (CRC) remain unclear. Here, we integrated multidimensional interaction resources, including gene expression profiling, protein-protein interactions (PPIs), transcriptional and post-transcriptional regulation data, and virus-host interactions, to tentatively explore potential molecular targets that functionally link UC and CRC at a systematic level. In this work, by deciphering the overlapping genes, crosstalking genes and pivotal regulators of both UC- and CRC-associated functional module pairs, we revealed a variety of genes (including FOS and DUSP1, etc.), transcription factors (including SMAD3 and ETS1, etc.) and miRNAs (including miR-155 and miR-196b, etc.) that may have the potential to complete the connections between UC and CRC. Interestingly, further analyses of the virus-host interaction network demonstrated that several virus proteins (including EBNA-LP of EBV and protein E7 of HPV) frequently inter-connected to UC- and CRC-associated module pairs with their validated targets significantly enriched in both modules of the host. Together, our results suggested that multidimensional integration strategy provides a novel approach to discover potential molecular targets that bridge the connections between UC and CRC, which could also be extensively applied to studies on other inflammation-related cancers.

A functional module-based exploration between inflammation and cancer in esophagus

Inflammation contributing to the underlying progression of diverse human cancers has been generally appreciated, however, explorations into the molecular links between inflammation and cancer in esophagus are still at its early stage. In our study, we presented a functional module-based approach, in combination with multiple data resource (gene expression, protein-protein interactions (PPI), transcriptional and post-transcriptional regulations) to decipher the underlying links. Via mapping differentially expressed disease genes, functional disease modules were identified. As indicated, those common genes and interactions tended to play important roles in linking inflammation and cancer. Based on crosstalk analysis, we demonstrated that, although most disease genes were not shared by both kinds of modules, they might act through participating in the same or similar functions to complete the molecular links. Additionally, we applied pivot analysis to extract significant regulators for per significant crosstalk module pair. As shown, pivot regulators might manipulate vital parts of the module subnetworks, and then work together to bridge inflammation and cancer in esophagus. Collectively, based on our functional module analysis, we demonstrated that shared genes or interactions, significant crosstalk modules, and those significant pivot regulators were served as different functional parts underlying the molecular links between inflammation and cancer in esophagus.

Individualized identification of disease-associated pathways with disrupted coordination of gene expression

Current pathway analysis approaches are primarily dedicated to capturing deregulated pathways at the population level and cannot provide patient-specific pathway deregulation information. In this article, the authors present a simple approach, called individPath, to detect pathways with significantly disrupted intra-pathway relative expression orderings for each disease sample compared with the stable, normal intra-pathway relative expression orderings pre-determined in previously accumulated normal samples. Through the analysis of multiple microarray data sets for lung and breast cancer, the authors demonstrate individPath's effectiveness for detecting cancer-associated pathways with disrupted relative expression orderings at the individual level and dissecting the heterogeneity of pathway deregulation among different patients. The portable use of this simple approach in clinical contexts is exemplified by the identification of prognostic intra-pathway gene pair signatures to predict overall survival of resected early-stage lung adenocarcinoma patients and signatures to predict relapse-free survival of estrogen receptor-positive breast cancer patients after tamoxifen treatment.

Individual-level analysis of differential expression of genes and pathways for personalized medicine

MOTIVATION

The differential expression analysis focusing on inter-group comparison can capture only differentially expressed genes (DE genes) at the population level, which may mask the heterogeneity of differential expression in individuals. Thus, to provide patient-specific information for personalized medicine, it is necessary to conduct differential expression analysis at the individual level.

RESULTS

We proposed a method to detect DE genes in individual disease samples by using the disrupted ordering in individual disease samples. In both simulated data and real paired cancer-normal sample data, this method showed excellent performance. It was found to be insensitive to experimental batch effects and data normalization. The landscape of stable gene pairs in a particular type of normal tissue could be predetermined using previously accumulated data, based on which dysregulated genes and pathways for any disease sample can be readily detected. The usefulness of the RankComp method in clinical settings was exemplified by the identification and application of prognostic markers for lung cancer.

AVAILABILITY AND IMPLEMENTATION

RankComp is implemented in R script that is freely available from Supplementary Materials.

Deciphering global signal features of high-throughput array data from cancers

Normalization of array data relies on the assumption that most genes are not altered, which means that the signals for different samples should be scaled to have similar median or average values. However, accumulating evidence suggests that gene expression could be widely up-regulated in cancers. Our previous results and subsequent findings have shown that violation of the assumption led to erroneous interpretation of microarray data. To decipher the global signal features of microarray data from cancer samples, we empirically evaluated a large collection of gene and miRNA expression profiles and copy-number variation arrays. Our results showed that, at the transcriptomic level, genes and miRNAs are widely over-expressed in a large proportion of cancers. In contrast, at the genomic level, global raw signal intensities for methylation and copy number variation show negligible differences between cancer and normal samples. These results force us to re-evaluate the proper use of normalization procedures under different experimental conditions and for different array platforms.

Functional comparison between genes dysregulated in ulcerative colitis and colorectal carcinoma

Background

Patients with ulcerative colitis (UC) are predisposed to colitis-associated colorectal cancer (CAC). However, the transcriptional mechanism of the transformation from UC to CAC is not fully understood.

Methodology

Firstly, we showed that CAC and non-UC-associated CRC were very similar in gene expression. Secondly, based on multiple datasets for UC and CRC, we extracted differentially expressed (DE) genes in UC and CRC versus normal controls, respectively. Thirdly, we compared the dysregulation directions (upregulation or downregulation) between DE genes of UC and CRC in CRC-related functions overrepresented with the DE genes of CRC, and proposed a regulatory model to explain the CRC-like dysregulation of genes in UC. A case study for “positive regulation of immune system process” was done to reveal the functional implication of DE genes with reversal dysregulations in these two diseases.

Principal Findings

In all the 44 detected CRC-related functions except for “viral transcription”, the dysregulation directions of DE genes in UC were significantly similar with their counterparts in CRC, and such CRC-like dysregulation in UC could be regulated by transcription factors affected by pro-inflammatory stimuli for colitis. A small portion of genes in each CRC-related function were dysregulated in opposite directions in the two diseases. The case study showed that genes related to humoral immunity specifically expressed in B cells tended to be upregulated in UC but downregulated in CRC.

Conclusions

The CRC-like dysregulation of genes in CRC-related functions in UC patients provides hints for understanding the transcriptional basis for UC to CRC transition. A small portion of genes with distinct dysregulation directions in each of the CRC-related functions in the two diseases implicate that their reversal dysregulations might be critical for UC to CRC transition. The cases study indicates that the humoral immune response might be inhibited during the transformation from UC to CRC.

Identification of human HK genes and gene expression regulation study in cancer from transcriptomics data analysis

The regulation of gene expression is essential for eukaryotes, as it drives the processes of cellular differentiation and morphogenesis, leading to the creation of different cell types in multicellular organisms. RNA-Sequencing (RNA-Seq) provides researchers with a powerful toolbox for characterization and quantification of transcriptome. Many different human tissue/cell transcriptome datasets coming from RNA-Seq technology are available on public data resource. The fundamental issue here is how to develop an effective analysis method to estimate expression pattern similarities between different tumor tissues and their corresponding normal tissues. We define the gene expression pattern from three directions: 1) expression breadth, which reflects gene expression on/off status, and mainly concerns ubiquitously expressed genes; 2) low/high or constant/variable expression genes, based on gene expression level and variation; and 3) the regulation of gene expression at the gene structure level. The cluster analysis indicates that gene expression pattern is higher related to physiological condition rather than tissue spatial distance. Two sets of human housekeeping (HK) genes are defined according to cell/tissue types, respectively. To characterize the gene expression pattern in gene expression level and variation, we firstly apply improved K-means algorithm and a gene expression variance model. We find that cancer-associated HK genes (a HK gene is specific in cancer group, while not in normal group) are expressed higher and more variable in cancer condition than in normal condition. Cancer-associated HK genes prefer to AT-rich genes, and they are enriched in cell cycle regulation related functions and constitute some cancer signatures. The expression of large genes is also avoided in cancer group. These studies will help us understand which cell type-specific patterns of gene expression differ among different cell types, and particularly for cancer.

Comparison of different normalization assumptions for analyses of DNA methylation data from the cancer genome

Nowadays, some researchers normalized DNA methylation arrays data in order to remove the technical artifacts introduced by experimental differences in sample preparation, array processing and other factors. However, other researchers analyzed DNA methylation arrays without performing data normalization considering that current normalizations for methylation data may distort real differences between normal and cancer samples because cancer genomes may be extensively subject to hypomethylation and the total amount of CpG methylation might differ substantially among samples. In this study, using eight datasets by Infinium HumanMethylation27 assay, we systemically analyzed the global distribution of DNA methylation changes in cancer compared to normal control and its effect on data normalization for selecting differentially methylated (DM) genes. We showed more differentially methylated (DM) genes could be found in the Quantile/Lowess-normalized data than in the non-normalized data. We found the DM genes additionally selected in the Quantile/Lowess-normalized data showed significantly consistent methylation states in another independent dataset for the same cancer, indicating these extra DM genes were effective biological signals related to the disease. These results suggested normalization can increase the power of detecting DM genes in the context of diagnostic markers which were usually characterized by relatively large effect sizes. Besides, we evaluated the reproducibility of DM discoveries for a particular cancer type, and we found most of the DM genes additionally detected in one dataset showed the same methylation directions in the other dataset for the same cancer type, indicating that these DM genes were effective biological signals in the other dataset. Furthermore, we showed that some DM genes detected from different studies for a particular cancer type were significantly reproducible at the functional level.