Batch correction of microarray data substantially improves the identification of genes differentially expressed in rheumatoid arthritis and osteoarthritis

Evaluation of normalization methods for predicting quantitative phenotypes in metagenomic data analysis

Genotype-to-phenotype mapping is an essential problem in the current genomic era. While qualitative case-control predictions have received significant attention, less emphasis has been placed on predicting quantitative phenotypes. This emerging field holds great promise in revealing intricate connections between microbial communities and host health. However, the presence of heterogeneity in microbiome datasets poses a substantial challenge to the accuracy of predictions and undermines the reproducibility of models. To tackle this challenge, we investigated 22 normalization methods that aimed at removing heterogeneity across multiple datasets, conducted a comprehensive review of them, and evaluated their effectiveness in predicting quantitative phenotypes in three simulation scenarios and 31 real datasets. The results indicate that none of these methods demonstrate significant superiority in predicting quantitative phenotypes or attain a noteworthy reduction in Root Mean Squared Error (RMSE) of the predictions. Given the frequent occurrence of batch effects and the satisfactory performance of batch correction methods in predicting datasets affected by these effects, we strongly recommend utilizing batch correction methods as the initial step in predicting quantitative phenotypes. In summary, the performance of normalization methods in predicting metagenomic data remains a dynamic and ongoing research area. Our study contributes to this field by undertaking a comprehensive evaluation of diverse methods and offering valuable insights into their effectiveness in predicting quantitative phenotypes.

A primer on correlation-based dimension reduction methods for multi-omics analysis

The continuing advances of omic technologies mean that it is now more tangible to measure the numerous features collectively reflecting the molecular properties of a sample. When multiple omic methods are used, statistical and computational approaches can exploit these large, connected profiles. Multi-omics is the integration of different omic data sources from the same biological sample. In this review, we focus on correlation-based dimension reduction approaches for single omic datasets, followed by methods for pairs of omics datasets, before detailing further techniques for three or more omic datasets. We also briefly detail network methods when three or more omic datasets are available and which complement correlation-oriented tools. To aid readers new to this area, these are all linked to relevant R packages that can implement these procedures. Finally, we discuss scenarios of experimental design and present road maps that simplify the selection of appropriate analysis methods. This review will help researchers navigate emerging methods for multi-omics and integrating diverse omic datasets appropriately. This raises the opportunity of implementing population multi-omics with large sample sizes as omics technologies and our understanding improve.

Batch normalization followed by merging is powerful for phenotype prediction integrating multiple heterogeneous studies

Heterogeneity in different genomic studies compromises the performance of machine learning models in cross-study phenotype predictions. Overcoming heterogeneity when incorporating different studies in terms of phenotype prediction is a challenging and critical step for developing machine learning algorithms with reproducible prediction performance on independent datasets. We investigated the best approaches to integrate different studies of the same type of omics data under a variety of different heterogeneities. We developed a comprehensive workflow to simulate a variety of different types of heterogeneity and evaluate the performances of different integration methods together with batch normalization by using ComBat. We also demonstrated the results through realistic applications on six colorectal cancer (CRC) metagenomic studies and six tuberculosis (TB) gene expression studies, respectively. We showed that heterogeneity in different genomic studies can markedly negatively impact the machine learning classifier's reproducibility. ComBat normalization improved the prediction performance of machine learning classifier when heterogeneous populations are present, and could successfully remove batch effects within the same population. We also showed that the machine learning classifier's prediction accuracy can be markedly decreased as the underlying disease model became more different in training and test populations. Comparing different merging and integration methods, we found that merging and integration methods can outperform each other in different scenarios. In the realistic applications, we observed that the prediction accuracy improved when applying ComBat normalization with merging or integration methods in both CRC and TB studies. We illustrated that batch normalization is essential for mitigating both population differences of different studies and batch effects. We also showed that both merging strategy and integration methods can achieve good performances when combined with batch normalization. In addition, we explored the potential of boosting phenotype prediction performance by rank aggregation methods and showed that rank aggregation methods had similar performance as other ensemble learning approaches.

Transcriptomic data meta-analysis reveals common and injury model specific gene expression changes in the regenerating zebrafish heart

Zebrafish have the capacity to fully regenerate the heart after an injury, which lies in sharp contrast to the irreversible loss of cardiomyocytes after a myocardial infarction in humans. Transcriptomics analysis has contributed to dissect underlying signaling pathways and gene regulatory networks in the zebrafish heart regeneration process. This process has been studied in response to different types of injuries namely: ventricular resection, ventricular cryoinjury, and genetic ablation of cardiomyocytes. However, there exists no database to compare injury specific and core cardiac regeneration responses. Here, we present a meta-analysis of transcriptomic data of regenerating zebrafish hearts in response to these three injury models at 7 days post injury (7dpi). We reanalyzed 36 samples and analyzed the differentially expressed genes (DEG) followed by downstream Gene Ontology Biological Processes (GO:BP) analysis. We found that the three injury models share a common core of DEG encompassing genes involved in cell proliferation, the Wnt signaling pathway and genes that are enriched in fibroblasts. We also found injury-specific gene signatures for resection and genetic ablation, and to a lower extent the cryoinjury model. Finally, we present our data in a user-friendly web interface that displays gene expression signatures across different injury types and highlights the importance to consider injury-specific gene regulatory networks when interpreting the results related to cardiac regeneration in the zebrafish. The analysis is freely available at: https://mybinder.org/v2/gh/MercaderLabAnatomy/PUB_Botos_et_al_2022_shinyapp_binder/HEAD?urlpath=shiny/bus-dashboard/ .

The importance of batch sensitization in missing value imputation

Data analysis is complex due to a myriad of technical problems. Amongst these, missing values and batch effects are endemic. Although many methods have been developed for missing value imputation (MVI) and batch correction respectively, no study has directly considered the confounding impact of MVI on downstream batch correction. This is surprising as missing values are imputed during early pre-processing while batch effects are mitigated during late pre-processing, prior to functional analysis. Unless actively managed, MVI approaches generally ignore the batch covariate, with unknown consequences. We examine this problem by modelling three simple imputation strategies: global (M1), self-batch (M2) and cross-batch (M3) first via simulations, and then corroborated on real proteomics and genomics data. We report that explicit consideration of batch covariates (M2) is important for good outcomes, resulting in enhanced batch correction and lower statistical errors. However, M1 and M3 are error-generating: global and cross-batch averaging may result in batch-effect dilution, with concomitant and irreversible increase in intra-sample noise. This noise is unremovable via batch correction algorithms and produces false positives and negatives. Hence, careless imputation in the presence of non-negligible covariates such as batch effects should be avoided.

MultiBaC: An R package to remove batch effects in multi-omic experiments

Motivation

Batch effects in omics datasets are usually a source of technical noise that masks the biological signal and hampers data analysis. Batch effect removal has been widely addressed for individual omics technologies. However, multi-omic datasets may combine data obtained in different batches where omics type and batch are often confounded. Moreover, systematic biases may be introduced without notice during data acquisition, which creates a hidden batch effect. Current methods fail to address batch effect correction in these cases.

Results

In this article, we introduce the MultiBaC R package, a tool for batch effect removal in multi-omics and hidden batch effect scenarios. The package includes a diversity of graphical outputs for model validation and assessment of the batch effect correction.

Availability and implementation

MultiBaC package is available on Bioconductor (https://www.bioconductor.org/packages/release/bioc/html/MultiBaC.html) and GitHub (https://github.com/ConesaLab/MultiBaC.git). The data underlying this article are available in Gene Expression Omnibus repository (accession numbers GSE11521, GSE1002, GSE56622 and GSE43747).

Supplementary information

Supplementary data are available at Bioinformatics online.

AMDBNorm: an approach based on distribution adjustment to eliminate batch effects of gene expression data

Batch effects explain a large part of the noise when merging gene expression data. Removing irrelevant variations introduced by batch effects plays an important role in gene expression studies. To obtain reliable differential analysis results, it is necessary to remove the variation caused by technical conditions between different batches while preserving biological variation. Usually, merging data directly with batch effects leads to a sharp rise in false positives. Although some methods of batch correction have been developed, they have some drawbacks. In this study, we develop a new algorithm, adjustment mean distribution-based normalization (AMDBNorm), which is based on a probability distribution to correct batch effects while preserving biological variation. AMDBNorm solves the defects of the existing batch correction methods. We compared several popular methods of batch correction with AMDBNorm using two real gene expression datasets with batch effects and analyzed the results of batch correction from the visual and quantitative perspectives. To ensure the biological variation was well protected, the effects of the batch correction methods were verified by hierarchical cluster analysis. The results showed that the AMDBNorm algorithm could remove batch effects of gene expression data effectively and retain more biological variation than other methods. Our approach provides the researchers with reliable data support in the study of differential gene expression analysis and prognostic biomarker selection.

Robustifying genomic classifiers to batch effects via ensemble learning

MOTIVATION

Genomic data are often produced in batches due to practical restrictions, which may lead to unwanted variation in data caused by discrepancies across batches. Such 'batch effects' often have negative impact on downstream biological analysis and need careful consideration. In practice, batch effects are usually addressed by specifically designed software, which merge the data from different batches, then estimate batch effects and remove them from the data. Here, we focus on classification and prediction problems, and propose a different strategy based on ensemble learning. We first develop prediction models within each batch, then integrate them through ensemble weighting methods.

RESULTS

We provide a systematic comparison between these two strategies using studies targeting diverse populations infected with tuberculosis. In one study, we simulated increasing levels of heterogeneity across random subsets of the study, which we treat as simulated batches. We then use the two methods to develop a genomic classifier for the binary indicator of disease status. We evaluate the accuracy of prediction in another independent study targeting a different population cohort. We observed that in independent validation, while merging followed by batch adjustment provides better discrimination at low level of heterogeneity, our ensemble learning strategy achieves more robust performance, especially at high severity of batch effects. These observations provide practical guidelines for handling batch effects in the development and evaluation of genomic classifiers.

AVAILABILITY AND IMPLEMENTATION

The data underlying this article are available in the article and in its online supplementary material. Processed data is available in the Github repository with implementation code, at https://github.com/zhangyuqing/bea_ensemble.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Ten-gene signature reveals the significance of clinical prognosis and immuno-correlation of osteosarcoma and study on novel skeleton inhibitors regarding MMP9

OBJECTIVES

This study aimed to identify novel targets in the carcinogenesis, therapy and prognosis of osteosarcoma from genomic level, together with screening ideal lead compounds with potential inhibition regarding MMP-9.

METHODS

Gene expression profiles from GSE12865, GSE14359, GSE33382, GSE36001 and GSE99671 were obtained respectively from GEO database. Differentially expressed genes were identified, and functional enrichment analysis, such as GO, KEGG, GSEA, PPI were performed to make a comprehensive understanding of the hub genes. Next, a series of high-precision computational techniques were conducted to screen potential lead compounds targeting MMP9, including virtual screening, ADME, toxicity prediction, and accurate docking analysis.

RESULTS

10 genes, MMP9, CD74, SPP1, CXCL12, TYROBP, FCER1G, HCLS1, ARHGDIB, LAPTM5 and IGF1R were identified as hub genes in the initiation of osteosarcoma. Machine learning, multivariate Cox analysis, ssGSEA and survival analysis demonstrated that these genes had values in prognosis, immune-correlation and targeted treatment. Tow novel compounds, ZINC000072131515 and ZINC000004228235, were screened as potential inhibitor regarding MMP9, and they could bind to MMP9 with favorable interaction energy and high binding affinity. Meanwhile, they were precited to be efficient and safe drugs with low-ames mutagenicity, none weight evidence of carcinogenicity, as well as non-toxic with liver.

CONCLUSIONS

This study revealed the significance of 10-gene signature in the development of osteosarcoma. Besides, drug candidates identified in this study provided a solid basis on MMP9 inhibitors' development.

A transfer learning approach to facilitate ComBat-based harmonization of multicentre radiomic features in new datasets

PURPOSE

To facilitate the demonstration of the prognostic value of radiomics, multicenter radiomics studies are needed. Pooling radiomic features of such data in a statistical analysis is however challenging, as they are sensitive to the variability in scanner models, acquisition protocols and reconstruction settings, which is often unavoidable in a multicentre retrospective analysis. A statistical harmonization strategy called ComBat was utilized in radiomics studies to deal with the "center-effect". The goal of the present work was to integrate a transfer learning (TL) technique within ComBat-and recently developed alternate versions of ComBat with improved flexibility (M-ComBat) and robustness (B-ComBat)-to allow the use of a previously determined harmonization transform to the radiomic feature values of new patients from an already known center.

MATERIAL AND METHODS

The proposed TL approach were incorporated in the four versions of ComBat (standard, B, M, and B-M ComBat). The proposed approach was evaluated using a dataset of 189 locally advanced cervical cancer patients from 3 centers, with magnetic resonance imaging (MRI) and positron emission tomography (PET) images, with the clinical endpoint of predicting local failure. The impact performance of the TL approach was evaluated by comparing the harmonization achieved using only parts of the data to the reference (harmonization achieved using all the available data). It was performed through three different machine learning pipelines.

RESULTS

The proposed TL technique was successful in harmonizing features of new patients from a known center in all versions of ComBat, leading to predictive models reaching similar performance as the ones developed using the features harmonized with all the data available.

CONCLUSION

The proposed TL approach enables applying a previously determined ComBat transform to new, previously unseen data.

Multivariate meta-analysis reveals global transcriptomic signatures underlying distinct human naive-like pluripotent states

The ground or naive pluripotent state of human pluripotent stem cells (hPSCs), which was initially established in mouse embryonic stem cells (mESCs), is an emerging and tentative concept. To verify this vital concept in hPSCs, we performed a multivariate meta-analysis of major hPSC datasets via the combined analytic powers of percentile normalization, principal component analysis (PCA), t-distributed stochastic neighbor embedding (t-SNE), and SC3 consensus clustering. This robust bioinformatics approach has significantly improved the predictive values of our meta-analysis. Accordingly, we revealed various similarities or dissimilarities between some naive-like hPSCs (NLPs) generated from different laboratories. Our analysis confirms some previous studies and provides new evidence concerning the existence of three distinct naive-like pluripotent states. Moreover, our study offers global transcriptomic markers that define diverse pluripotent states under various hPSC growth protocols.

Identify potential miRNA-mRNA regulatory networks contributing to high-risk neuroblastoma

Neuroblastoma (NB) is a common tumor in children, usually in the retroperitoneum. After various treatments, low- and intermediate-risk patients have achieved good results, but the prognosis of high-risk patients is still very poor. Therefore, it is necessary to find new effective targets for the treatment of high-risk patients. In this study, comprehensive bioinformatics analysis was used to identify the differentially expressed genes (DEG and DEM) between high-risk patients and non-high-risk patients, and it was identified that ADRB2 may affect the survival status of high-risk patients due to miR -30a-5p regulation. The GSE49710, GSE73517, and GSE121513 datasets were downloaded from the Gene Expression Synthesis (GEO) database, and DEG and DEM were selected. Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were applied to the selected DEGs. The STRING database and Cytoscape software were used to construct protein-protein interaction (PPI) networks and perform modular analysis of the DEGs. The TARGET data set containing information on overall survival days were used for the prognostic analysis of central genes. We identified a total of 255 DEGs from GSE49710 and GSE73517, and 193 DEMs from GSE121513. We identified the 5 most important central genes from the PPI network, performed a prognostic analysis in the target data set, and verified their expression using RT-qPCR to select the most important ADRB2 gene to predict miRNA. Integrating the differential miRNA predicted by miRDB and miRSystem and GSE121513 between the targeted miRNA and the prognosis, miR-30a-5p was finally identified as the targeted miRNA of ADRB2.

Transcriptome profiling reveals target in primary myelofibrosis together with structural biology study on novel natural inhibitors regarding JAK2

This study aimed to identify effective targets for carcinogenesis of primary myelofibrosis (PMF), as well as to screen ideal lead compounds with potential inhibition effect on Janus kinase 2 to contribute to the medication design and development. Gene expression profiles of GSE26049, GSE53482, GSE61629 were obtained from the Gene Expression Omnibus database. The differentially expressed genes were identified, and functional enrichment analyses such as Gene Ontology, protein-protein interaction network etc., were performed step by step. Subsequently, highly-precise computational techniques were conducted to identify potential inhibitors of JAK2. A series of structural biology methods including virtual screening, ADMET (absorption, distribution, metabolism, excretion, and toxicity) prediction, molecule docking, molecular dynamics simulation etc., were implemented to discover novel natural compounds. Results elucidated that PMF patients had abnormal LCN2, JAK2, MMP8, CAMP, DEFA4, LTF, MPO, HBD, STAT4, EBF1 mRNA expression compared to normal patients. Functional enrichment analysis revealed that these genes were mainly enriched in erythrocyte differentiation, neutrophil degranulation and killing cells of other organisms. Two novel natural compounds, ZINC000013513540 and ZINC000004099068 were found binding to JAK2 with favorable interaction energy together with high binding affinity. They were predicted with non-Ames mutagenicity, low-rodent carcinogenicity, less developmental toxicity potential as well as non-toxicity with liver. Molecular dynamics simulation demonstrated that these two complexes: ZINC000013513540-JAK2 and ZINC000004099068-JAK2 could exist stably under natural circumstances. In conclusion, this study revealed hub genes in the carcinogenesis of PMF. ZINC000013513540 and ZINC000004099068 were promising drugs in dealing with PMF. This study may also accelerate exploration of new drugs.

DRIM: A Web-Based System for Investigating Drug Response at the Molecular Level by Condition-Specific Multi-Omics Data Integration

Pharmacogenomics is the study of how genes affect a person's response to drugs. Thus, understanding the effect of drug at the molecular level can be helpful in both drug discovery and personalized medicine. Over the years, transcriptome data upon drug treatment has been collected and several databases compiled before drug treatment cancer cell multi-omics data with drug sensitivity (IC 50, AUC) or time-series transcriptomic data after drug treatment. However, analyzing transcriptome data upon drug treatment is challenging since more than 20,000 genes interact in complex ways. In addition, due to the difficulty of both time-series analysis and multi-omics integration, current methods can hardly perform analysis of databases with different data characteristics. One effective way is to interpret transcriptome data in terms of well-characterized biological pathways. Another way is to leverage state-of-the-art methods for multi-omics data integration. In this paper, we developed Drug Response analysis Integrating Multi-omics and time-series data (DRIM), an integrative multi-omics and time-series data analysis framework that identifies perturbed sub-pathways and regulation mechanisms upon drug treatment. The system takes drug name and cell line identification numbers or user's drug control/treat time-series gene expression data as input. Then, analysis of multi-omics data upon drug treatment is performed in two perspectives. For the multi-omics perspective analysis, IC 50-related multi-omics potential mediator genes are determined by embedding multi-omics data to gene-centric vector space using a tensor decomposition method and an autoencoder deep learning model. Then, perturbed pathway analysis of potential mediator genes is performed. For the time-series perspective analysis, time-varying perturbed sub-pathways upon drug treatment are constructed. Additionally, a network involving transcription factors (TFs), multi-omics potential mediator genes, and perturbed sub-pathways is constructed, and paths to perturbed pathways from TFs are determined by an influence maximization method. To demonstrate the utility of our system, we provide analysis results of sub-pathway regulatory mechanisms in breast cancer cell lines of different drug sensitivity. DRIM is available at: http://biohealth.snu.ac.kr/software/DRIM/.

Silencing of SPP1 Suppresses Progression of Tongue Cancer by Mediating the PI3K/Akt Signaling Pathway

Background:

In the present study, we aimed to find an effective target for the treatment of tongue cancer using gene chip screening and signal pathway research.

Methods:

We used microarray screening and gene expression profile analyses to find important differentially expressed genes in tongue cancer. We constructed a protein-protein interaction network, and used enrichment analysis of the Kyoto Encyclopedia of Genes and Genomes to screen for important genes. We then silenced the genes of interest in SCC154 cells to study the relationship with the Phosphatidylinositol 3-kinase/Akt signal pathway. Western blot analyses, the 3-(4,5Dimethylthiazol-yl)-2,5Dimethylthiazol-2-yl)-2,5diphenyltetrazolium bromide (MTT) test, and immunofluorescence assays were used to compare the expression levels of Phosphatidylinositol 3-kinase/Akt signal pathway-related proteins, cell viability, and cell proliferation ability in normal SCC154 cells, Si-RNA SCC154 cells, and gene-silenced SCC154 cells. The scratch test, Transwell test, and western blotting were used to determine migration, invasion, and carcinogenesis.

Results:

Using GSE9844, GSE13601, and GSE31056 gene chips, we identified 93 upregulated genes and 76 downregulated genes in tongue cancer. Using the protein-protein interaction network and Kyoto Encyclopedia of Genes and Genomes enrichment analyses, we further identified 47 differentially expressed genes. Using Kaplan-Meier plotter online tools, we also identified 3 genes (SPP1, Recombinant Human Secreted Phosphoprotein 1; PLAU, plasminogen activator urinary; and APP, amyloid precursor protein). Compared with normal SCC154 cells and Si-RNA control SCC154 cells, the expressions of Phosphatidylinositol 3-kinase/Akt pathway proteins in si-SPP1 SCC154 cells were significantly decreased (*P < 0.05), and the protein activities and proliferation abilities were also significantly decreased (*P < 0.05), while the migration ability, invasion ability, and cancer forming ability were significantly increased (*P < 0.05).

Conclusion:

Inhibition of the SPP1 gene may have a therapeutic effect on tongue cancer, and could be an effective target for the treatment of this disorder.

Analysis and identification of novel biomarkers involved in neuroblastoma via integrated bioinformatics

Neuroblastoma (NB) is the most common extracranial solid tumor in children. Under various treatments, some patients still have a poor prognosis. Hence, it is necessary to find new valid targets for NB therapy. In this study, a comprehensive bioinformatic analysis was used to identify differentially expressed genes (DEGs) between NB and control cells, and to select hub genes associated with NB. GSE66586 and GSE78061 datasets were downloaded from the Gene Expression Omnibus (GEO) database and DEGs were selected. Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were applied to the selected DEGs. The STRING database and Cytoscape software were used to construct protein-protein interaction (PPI) networks and perform modular analysis of the DEGs. The R2 database was used for prognostic analysis. We identified a total of 238 DEGs from two microarray databases. GO enrichment analysis shows that these DEGs are mainly concentrated in the regulation of cell growth, cell migration, cell fate determination, and cell maturation. KEGG pathway analysis showed that these DEGs are mainly involved in focal adhesion, the TNF signaling pathway, cancer-related pathways, and signaling pathways regulating stem cell pluripotency. We identified the 15 most closely related DEGs from the PPI network, and performed R2 database prognostic analysis to select five hub genes - CTGF, EDN1, GATA2, LOX, and SERPINE1. This study distinguished hub genes and related signaling pathways that can potentially serve as diagnostic indicators and therapeutic biomarkers for NB, thereby improving understanding of the molecular mechanisms involved in NB.

Merging microarray studies to identify a common gene expression signature to several structural heart diseases

BACKGROUND

Heart disease is the leading cause of death worldwide. Knowing a gene expression signature in heart disease can lead to the development of more efficient diagnosis and treatments that may prevent premature deaths. A large amount of microarray data is available in public repositories and can be used to identify differentially expressed genes. However, most of the microarray datasets are composed of a reduced number of samples and to obtain more reliable results, several datasets have to be merged, which is a challenging task. The identification of differentially expressed genes is commonly done using statistical methods. Nonetheless, these methods are based on the definition of an arbitrary threshold to select the differentially expressed genes and there is no consensus on the values that should be used.

RESULTS

Nine publicly available microarray datasets from studies of different heart diseases were merged to form a dataset composed of 689 samples and 8354 features. Subsequently, the adjusted p-value and fold change were determined and by combining a set of adjusted p-values cutoffs with a list of different fold change thresholds, 12 sets of differentially expressed genes were obtained. To select the set of differentially expressed genes that has the best accuracy in classifying samples from patients with heart diseases and samples from patients with no heart condition, the random forest algorithm was used. A set of 62 differentially expressed genes having a classification accuracy of approximately 95% was identified.

CONCLUSIONS

We identified a gene expression signature common to different cardiac diseases and supported our findings by showing their involvement in the pathophysiology of the heart. The approach used in this study is suitable for the identification of gene expression signatures, and can be extended to different diseases.

Performance comparison of modified ComBat for harmonization of radiomic features for multicenter studies

Multicenter studies are needed to demonstrate the clinical potential value of radiomics as a prognostic tool. However, variability in scanner models, acquisition protocols and reconstruction settings are unavoidable and radiomic features are notoriously sensitive to these factors, which hinders pooling them in a statistical analysis. A statistical harmonization method called ComBat was developed to deal with the "batch effect" in gene expression microarray data and was used in radiomics studies to deal with the "center-effect". Our goal was to evaluate modifications in ComBat allowing for more flexibility in choosing a reference and improving robustness of the estimation. Two modified ComBat versions were evaluated: M-ComBat allows to transform all features distributions to a chosen reference, instead of the overall mean, providing more flexibility. B-ComBat adds bootstrap and Monte Carlo for improved robustness in the estimation. BM-ComBat combines both modifications. The four versions were compared regarding their ability to harmonize features in a multicenter context in two different clinical datasets. The first contains 119 locally advanced cervical cancer patients from 3 centers, with magnetic resonance imaging and positron emission tomography imaging. In that case ComBat was applied with 3 labels corresponding to each center. The second one contains 98 locally advanced laryngeal cancer patients from 5 centers with contrast-enhanced computed tomography. In that specific case, because imaging settings were highly heterogeneous even within each of the five centers, unsupervised clustering was used to determine two labels for applying ComBat. The impact of each harmonization was evaluated through three different machine learning pipelines for the modelling step in predicting the clinical outcomes, across two performance metrics (balanced accuracy and Matthews correlation coefficient). Before harmonization, almost all radiomic features had significantly different distributions between labels. These differences were successfully removed with all ComBat versions. The predictive ability of the radiomic models was always improved with harmonization and the improved ComBat provided the best results. This was observed consistently in both datasets, through all machine learning pipelines and performance metrics. The proposed modifications allow for more flexibility and robustness in the estimation. They also slightly but consistently improve the predictive power of resulting radiomic models.

Transcriptomics in Toxicogenomics, Part II: Preprocessing and Differential Expression Analysis for High Quality Data

Preprocessing of transcriptomics data plays a pivotal role in the development of toxicogenomics-driven tools for chemical toxicity assessment. The generation and exploitation of large volumes of molecular profiles, following an appropriate experimental design, allows the employment of toxicogenomics (TGx) approaches for a thorough characterisation of the mechanism of action (MOA) of different compounds. To date, a plethora of data preprocessing methodologies have been suggested. However, in most cases, building the optimal analytical workflow is not straightforward. A careful selection of the right tools must be carried out, since it will affect the downstream analyses and modelling approaches. Transcriptomics data preprocessing spans across multiple steps such as quality check, filtering, normalization, batch effect detection and correction. Currently, there is a lack of standard guidelines for data preprocessing in the TGx field. Defining the optimal tools and procedures to be employed in the transcriptomics data preprocessing will lead to the generation of homogeneous and unbiased data, allowing the development of more reliable, robust and accurate predictive models. In this review, we outline methods for the preprocessing of three main transcriptomic technologies including microarray, bulk RNA-Sequencing (RNA-Seq), and single cell RNA-Sequencing (scRNA-Seq). Moreover, we discuss the most common methods for the identification of differentially expressed genes and to perform a functional enrichment analysis. This review is the second part of a three-article series on Transcriptomics in Toxicogenomics.

MultiBaC: A strategy to remove batch effects between different omic data types

Diversity of omic technologies has expanded in the last years together with the number of omic data integration strategies. However, multiomic data generation is costly, and many research groups cannot afford research projects where many different omic techniques are generated, at least at the same time. As most researchers share their data in public repositories, different omic datasets of the same biological system obtained at different labs can be combined to construct a multiomic study. However, data obtained at different labs or moments in time are typically subjected to batch effects that need to be removed for successful data integration. While there are methods to correct batch effects on the same data types obtained in different studies, they cannot be applied to correct lab or batch effects across omics. This impairs multiomic meta-analysis. Fortunately, in many cases, at least one omics platform-i.e. gene expression- is repeatedly measured across labs, together with the additional omic modalities that are specific to each study. This creates an opportunity for batch analysis. We have developed MultiBaC (multiomic Multiomics Batch-effect Correction correction), a strategy to correct batch effects from multiomic datasets distributed across different labs or data acquisition events. Our strategy is based on the existence of at least one shared data type which allows data prediction across omics. We validate this approach both on simulated data and on a case where the multiomic design is fully shared by two labs, hence batch effect correction within the same omic modality using traditional methods can be compared with the MultiBaC correction across data types. Finally, we apply MultiBaC to a true multiomic data integration problem to show that we are able to improve the detection of meaningful biological effects.

Screening responsive or resistant biomarkers of immune checkpoint inhibitors based on online databases

Immune checkpoint inhibitors are a promising strategy in the treatment of cancer, especially advanced types. However, not all patients are responsive to immune checkpoint inhibitors. The response rate depends on the immune microenvironment, tumor mutational burden (TMB), expression level of immune checkpoint proteins, and molecular subtypes of cancers. Along with the Cancer Genome Project, various open access databases, including The Cancer Genome Atlas and Gene Expression Omnibus, provide large volumes of data, which allow researchers to explore responsive or resistant biomarkers of immune checkpoint inhibitors. In this review, we introduced some methodologies on database selection, biomarker screening, current progress of immune checkpoint blockade in solid tumor treatment, possible mechanisms of drug resistance, strategies of overcoming resistance, and indications for immune checkpoint inhibitor therapy.

Alternative empirical Bayes models for adjusting for batch effects in genomic studies

Background

Combining genomic data sets from multiple studies is advantageous to increase statistical power in studies where logistical considerations restrict sample size or require the sequential generation of data. However, significant technical heterogeneity is commonly observed across multiple batches of data that are generated from different processing or reagent batches, experimenters, protocols, or profiling platforms. These so-called batch effects often confound true biological relationships in the data, reducing the power benefits of combining multiple batches, and may even lead to spurious results in some combined studies. Therefore there is significant need for effective methods and software tools that account for batch effects in high-throughput genomic studies.

Results

Here we contribute multiple methods and software tools for improved combination and analysis of data from multiple batches. In particular, we provide batch effect solutions for cases where the severity of the batch effects is not extreme, and for cases where one high-quality batch can serve as a reference, such as the training set in a biomarker study. We illustrate our approaches and software in both simulated and real data scenarios.

Conclusions

We demonstrate the value of these new contributions compared to currently established approaches in the specified batch correction situations.

Electronic supplementary material

The online version of this article (10.1186/s12859-018-2263-6) contains supplementary material, which is available to authorized users.

Integrated bioinformatics analysis reveals key candidate genes and pathways in breast cancer

Breast cancer (BC) is the leading malignancy in women worldwide, yet relatively little is known about the genes and signaling pathways involved in BC tumorigenesis and progression. The present study aimed to elucidate potential key candidate genes and pathways in BC. Five gene expression profile data sets (GSE22035, GSE3744, GSE5764, GSE21422 and GSE26910) were downloaded from the Gene Expression Omnibus (GEO) database, which included data from 113 tumorous and 38 adjacent non-tumorous tissue samples. Differentially expressed genes (DEGs) were identified using t-tests in the limma R package. These DEGs were subsequently investigated by pathway enrichment analysis and a protein-protein interaction (PPI) network was constructed. The most significant module from the PPI network was selected for pathway enrichment analysis. In total, 227 DEGs were identified, of which 82 were upregulated and 145 were downregulated. Pathway enrichment analysis results revealed that the upregulated DEGs were mainly enriched in ‘cell division’, the ‘proteinaceous extracellular matrix (ECM)’, ‘ECM structural constituents’ and ‘ECM-receptor interaction’, whereas downregulated genes were mainly enriched in ‘response to drugs’, ‘extracellular space’, ‘transcriptional activator activity’ and the ‘peroxisome proliferator-activated receptor signaling pathway’. The PPI network contained 174 nodes and 1,257 edges. DNA topoisomerase 2-a, baculoviral inhibitor of apoptosis repeat-containing protein 5, cyclin-dependent kinase 1, G2/mitotic-specific cyclin-B1 and kinetochore protein NDC80 homolog were identified as the top 5 hub genes. Furthermore, the genes in the most significant module were predominantly involved in ‘mitotic nuclear division’, ‘mid-body’, ‘protein binding’ and ‘cell cycle’. In conclusion, the DEGs, relative pathways and hub genes identified in the present study may aid in understanding of the molecular mechanisms underlying BC progression and provide potential molecular targets and biomarkers for BC.

Why Batch Effects Matter in Omics Data, and How to Avoid Them

Effective integration and analysis of new high-throughput data, especially gene-expression and proteomic-profiling data, are expected to deliver novel clinical insights and therapeutic options. Unfortunately, technical heterogeneity or batch effects (different experiment times, handlers, reagent lots, etc.) have proven challenging. Although batch effect-correction algorithms (BECAs) exist, we know little about effective batch-effect mitigation: even now, new batch effect-associated problems are emerging. These include false effects due to misapplying BECAs and positive bias during model evaluations. Depending on the choice of algorithm and experimental set-up, biological heterogeneity can be mistaken for batch effects and wrongfully removed. Here, we examine these emerging batch effect-associated problems, propose a series of best practices, and discuss some of the challenges that lie ahead.

Unveiling gene trait relationship by cross-platform meta-analysis on Chinese hamster ovary cell transcriptome

In the past few years, transcriptome analysis has been increasingly employed to better understand the physiology of Chinese hamster ovary (CHO) cells at a global level. As more transcriptome data accumulated, meta-analysis on data sets collected from various sources can potentially provide better insights on common properties of those cells. Here, we performed meta-analysis on transcriptome data of different CHO cell lines obtained using NimbleGen or Affymetrix microarray platforms. Hierarchical clustering, non-negative matrix factorization (NMF) analysis, and principal component analysis (PCA) accordantly showed the samples were clustered into two groups: one consists of adherent cells in serum-containing medium, and the other suspension cells in serum-free medium. Genes that were differentially expressed between the two clusters were enriched in a few functional classes by Database for Annotation, Visualization, and Integrated Discovery (DAVID) of which many were common with the enriched gene sets identified by Gene Set Enrichment Analysis (GSEA), including extracellular matrix (ECM) receptor interaction, cell adhesion molecules (CAMs), and lipid related metabolism pathways. Despite the heterogeneous sources of the cell samples, the adherent and suspension growth characteristics and serum-supplementation appear to be a dominant feature in the transcriptome. The results demonstrated that meta-analysis of transcriptome could uncover features in combined data sets that individual data set might not reveal. As transcriptome data sets accumulate over time, meta-analysis will become even more revealing. Biotechnol. Bioeng. 2017;114: 1583-1592. © 2017 Wiley Periodicals, Inc.

Re-evaluating data quality of dog mitochondrial, Y chromosomal, and autosomal SNPs genotyped by SNP array

Quality deficiencies in single nucleotide polymorphism (SNP) analyses have important implications. We used missingness rates to investigate the quality of a recently published dataset containing 424 mitochondrial, 211 Y chromosomal, and 160 432 autosomal SNPs generated by a semicustom Illumina SNP array from 5 392 dogs and 14 grey wolves. Overall, the individual missingness rate for mitochondrial SNPs was ~43.8%, with 980 (18.1%) individuals completely missing mitochondrial SNP genotyping (missingness rate=1). In males, the genotype missingness rate was ~28.8% for Y chromosomal SNPs, with 374 males recording rates above 0.96. These 374 males also exhibited completely failed mitochondrial SNPs genotyping, indicative of a batch effect. Individual missingness rates for autosomal markers were greater than zero, but less than 0.5. Neither mitochondrial nor Y chromosomal SNPs achieved complete genotyping (locus missingness rate=0), whereas 5.9% of autosomal SNPs had a locus missingness rate=1. The high missingness rates and possible batch effect show that caution and rigorous measures are vital when genotyping and analyzing SNP array data for domestic animals. Further improvements of these arrays will be helpful to future studies.

A urinary microRNA signature can predict the presence of bladder urothelial carcinoma in patients undergoing surveillance

BACKGROUND

The objective of this study was to determine whether microRNA (miRNA) profiling of urine could identify the presence of urothelial carcinoma of the bladder (UCB) and to compare its performance characteristics to that of cystoscopy.

METHODS

In the discovery cohort we screened 81 patients, which included 21 benign controls, 30 non-recurrers and 30 patients with active cancer (recurrers), using a panel of 12 miRNAs. Data analysis was performed using a machine learning approach of a Support Vector Machine classifier with a Student's t-test feature selection procedure. This was trained using a three-fold cross validation approach and performance was measured using the area under the receiver operator characteristic curve (AUC). The miRNA signature was validated in an independent cohort of a further 50 patients.

RESULTS

The best predictor to distinguish patients with UCB from non-recurrers was achieved using a combination of six miRNAs (AUC=0.85). This validated in an independent cohort (AUC=0.74) and detected UCB with a high sensitivity (88%) and sufficient specificity (48%) with all significant cancers identified. The performance of the classifier was best in detecting clinically significant disease such as presence of T1 Stage disease (AUC=0.92) and high-volume disease (AUC=0.81). Cystoscopy rates in the validation cohort would have been reduced by 30%.

CONCLUSIONS

Urinary profiling using this panel of miRNAs shows promise for detection of tumour recurrence in the surveillance of UCB. Such a panel may be useful in reducing the morbidity and costs associated with cystoscopic surveillance, and now merits prospective evaluation.

Methods that remove batch effects while retaining group differences may lead to exaggerated confidence in downstream analyses

Removal of, or adjustment for, batch effects or center differences is generally required when such effects are present in data. In particular, when preparing microarray gene expression data from multiple cohorts, array platforms, or batches for later analyses, batch effects can have confounding effects, inducing spurious differences between study groups. Many methods and tools exist for removing batch effects from data. However, when study groups are not evenly distributed across batches, actual group differences may induce apparent batch differences, in which case batch adjustments may bias, usually deflate, group differences. Some tools therefore have the option of preserving the difference between study groups, e.g. using a two-way ANOVA model to simultaneously estimate both group and batch effects. Unfortunately, this approach may systematically induce incorrect group differences in downstream analyses when groups are distributed between the batches in an unbalanced manner. The scientific community seems to be largely unaware of how this approach may lead to false discoveries.

Association of Protein Translation and Extracellular Matrix Gene Sets with Breast Cancer Metastasis: Findings Uncovered on Analysis of Multiple Publicly Available Datasets Using Individual Patient Data Approach

Introduction

Microarray analysis has revolutionized the role of genomic prognostication in breast cancer. However, most studies are single series studies, and suffer from methodological problems. We sought to use a meta-analytic approach in combining multiple publicly available datasets, while correcting for batch effects, to reach a more robust oncogenomic analysis.

Aim

The aim of the present study was to find gene sets associated with distant metastasis free survival (DMFS) in systemically untreated, node-negative breast cancer patients, from publicly available genomic microarray datasets.

Methods

Four microarray series (having 742 patients) were selected after a systematic search and combined. Cox regression for each gene was done for the combined dataset (univariate, as well as multivariate – adjusted for expression of Cell cycle related genes) and for the 4 major molecular subtypes. The centre and microarray batch effects were adjusted by including them as random effects variables. The Cox regression coefficients for each analysis were then ranked and subjected to a Gene Set Enrichment Analysis (GSEA).

Results

Gene sets representing protein translation were independently negatively associated with metastasis in the Luminal A and Luminal B subtypes, but positively associated with metastasis in Basal tumors. Proteinaceous extracellular matrix (ECM) gene set expression was positively associated with metastasis, after adjustment for expression of cell cycle related genes on the combined dataset. Finally, the positive association of the proliferation-related genes with metastases was confirmed.

Conclusion

To the best of our knowledge, the results depicting mixed prognostic significance of protein translation in breast cancer subtypes are being reported for the first time. We attribute this to our study combining multiple series and performing a more robust meta-analytic Cox regression modeling on the combined dataset, thus discovering 'hidden' associations. This methodology seems to yield new and interesting results and may be used as a tool to guide new research.

Removing batch effects from purified plasma cell gene expression microarrays with modified ComBat

Background

Gene expression profiling (GEP) via microarray analysis is a widely used tool for assessing risk and other patient diagnostics in clinical settings. However, non-biological factors such as systematic changes in sample preparation, differences in scanners, and other potential batch effects are often unavoidable in long-term studies and meta-analysis. In order to reduce the impact of batch effects on microarray data, Johnson, Rabinovic, and Li developed ComBat for use when combining batches of gene expression microarray data.

We propose a modification to ComBat that centers data to the location and scale of a pre-determined, ‘gold-standard’ batch. This modified ComBat (M-Combat) is designed specifically in the context of meta-analysis and batch effect adjustment for use with predictive models that are validated and fixed on historical data from a ‘gold-standard’ batch.

Results

We combined data from MIRT across two batches (‘Old’ and ‘New’ Kit sample preparation) as well as external data sets from the HOVON-65/GMMG-HD4 and MRC-IX trials into a combined set, first without transformation and then with both ComBat and M-ComBat transformations. Fixed and validated gene risk signatures developed at MIRT on the Old Kit standard (GEP5, GEP70, and GEP80 risk scores) were compared across these combined data sets.

Both ComBat and M-ComBat eliminated all of the differences among probes caused by systematic batch effects (over 98% of all untransformed probes were significantly different by ANOVA with 0.01 q-value threshold reduced to zero significant probes with ComBat and M-ComBat). The agreement in mean and distribution of risk scores, as well as the proportion of high-risk subjects identified, coincided with the ‘gold-standard’ batch more with M-ComBat than with ComBat. The performance of risk scores improved overall using either ComBat or M-Combat; however, using M-ComBat and the original, optimal risk cutoffs allowed for greater ability in our study to identify smaller cohorts of high-risk subjects.

Conclusion

M-ComBat is a practical modification to an accepted method that offers greater power to control the location and scale of batch-effect adjusted data. M-ComBat allows for historical models to function as intended on future samples despite known, often unavoidable systematic changes to gene expression data.

Novel application of multi-stimuli network inference to synovial fibroblasts of rheumatoid arthritis patients

Background

Network inference of gene expression data is an important challenge in systems biology. Novel algorithms may provide more detailed gene regulatory networks (GRN) for complex, chronic inflammatory diseases such as rheumatoid arthritis (RA), in which activated synovial fibroblasts (SFBs) play a major role. Since the detailed mechanisms underlying this activation are still unclear, simultaneous investigation of multi-stimuli activation of SFBs offers the possibility to elucidate the regulatory effects of multiple mediators and to gain new insights into disease pathogenesis.

Methods

A GRN was therefore inferred from RA-SFBs treated with 4 different stimuli (IL-1 β, TNF- α, TGF- β, and PDGF-D). Data from time series microarray experiments (0, 1, 2, 4, 12 h; Affymetrix HG-U133 Plus 2.0) were batch-corrected applying ‘ComBat’, analyzed for differentially expressed genes over time with ‘Limma’, and used for the inference of a robust GRN with NetGenerator V2.0, a heuristic ordinary differential equation-based method with soft integration of prior knowledge.

Results

Using all genes differentially expressed over time in RA-SFBs for any stimulus, and selecting the genes belonging to the most significant gene ontology (GO) term, i.e., ‘cartilage development’, a dynamic, robust, moderately complex multi-stimuli GRN was generated with 24 genes and 57 edges in total, 31 of which were gene-to-gene edges. Prior literature-based knowledge derived from Pathway Studio or manual searches was reflected in the final network by 25/57 confirmed edges (44%). The model contained known network motifs crucial for dynamic cellular behavior, e.g., cross-talk among pathways, positive feed-back loops, and positive feed-forward motifs (including suppression of the transcriptional repressor OSR2 by all 4 stimuli.

Conclusion

A multi-stimuli GRN highly concordant with literature data was successfully generated by network inference from the gene expression of stimulated RA-SFBs. The GRN showed high reliability, since 10 predicted edges were independently validated by literature findings post network inference. The selected GO term ‘cartilage development’ contained a number of differentiation markers, growth factors, and transcription factors with potential relevance for RA. Finally, the model provided new insight into the response of RA-SFBs to multiple stimuli implicated in the pathogenesis of RA, in particular to the ‘novel’ potent growth factor PDGF-D.

Integrative genomic analysis identifies epigenetic marks that mediate genetic risk for epithelial ovarian cancer

BACKGROUND

Both genetic and epigenetic factors influence the development and progression of epithelial ovarian cancer (EOC). However, there is an incomplete understanding of the interrelationship between these factors and the extent to which they interact to impact disease risk. In the present study, we aimed to gain insight into this relationship by identifying DNA methylation marks that are candidate mediators of ovarian cancer genetic risk.

METHODS

We used 214 cases and 214 age-matched controls from the Mayo Clinic Ovarian Cancer Study. Pretreatment, blood-derived DNA was profiled for genome-wide methylation (Illumina Infinium HumanMethylation27 BeadArray) and single nucleotide polymorphisms (SNPs, Illumina Infinium HD Human610-Quad BeadArray). The Causal Inference Test (CIT) was implemented to distinguish CpG sites that mediate genetic risk, from those that are consequential or independently acted on by genotype.

RESULTS

Controlling for the estimated distribution of immune cells and other key covariates, our initial epigenome-wide association analysis revealed 1,993 significantly differentially methylated CpGs that between cases and controls (FDR, q < 0.05). The relationship between methylation and case-control status for these 1,993 CpGs was found to be highly consistent with the results of previously published, independent study that consisted of peripheral blood DNA methylation signatures in 131 pretreatment cases and 274 controls. Implementation of the CIT test revealed 17 CpG/SNP pairs, comprising 13 unique CpGs and 17 unique SNPs, which represent potential methylation-mediated relationships between genotype and EOC risk. Of these 13 CpGs, several are associated with immune related genes and genes that have been previously shown to exhibit altered expression in the context of cancer.

CONCLUSIONS

These findings provide additional insight into EOC etiology and may serve as novel biomarkers for EOC susceptibility.

Congestive heart failure effects on atrial fibroblast phenotype: differences between freshly-isolated and cultured cells

Introduction

Fibroblasts are important in the atrial fibrillation (AF) substrate resulting from congestive heart failure (CHF). We previously noted changes in in vivo indices of fibroblast function in a CHF dog model, but could not detect changes in isolated cells. This study assessed CHF-induced changes in the phenotype of fibroblasts freshly isolated from control versus CHF dogs, and examined effects of cell culture on these differences.

Methods/Results

Left-atrial fibroblasts were isolated from control and CHF dogs (ventricular tachypacing 240 bpm×2 weeks). Freshly-isolated fibroblasts were compared to fibroblasts in primary culture. Extracellular-matrix (ECM) gene-expression was assessed by qPCR, protein by Western blot, fibroblast morphology with immunocytochemistry, and K⁺-current with patch-clamp. Freshly-isolated CHF fibroblasts had increased expression-levels of collagen-1 (10-fold), collagen-3 (5-fold), and fibronectin-1 (3-fold) vs. control, along with increased cell diameter (13.4±0.4 µm vs control 8.4±0.3 µm) and cell spreading (shape factor 0.81±0.02 vs. control 0.87±0.02), consistent with an activated phenotype. Freshly-isolated control fibroblasts displayed robust tetraethylammonium (TEA)-sensitive K⁺-currents that were strongly downregulated in CHF. The TEA-sensitive K⁺-current differences between control and CHF fibroblasts were attenuated after 2-day culture and eliminated after 7 days. Similarly, cell-culture eliminated the ECM protein-expression and shape differences between control and CHF fibroblasts.

Conclusions

Freshly-isolated CHF and control atrial fibroblasts display distinct ECM-gene and morphological differences consistent with in vivo pathology. Culture for as little as 48 hours activates fibroblasts and obscures the effects of CHF. These results demonstrate potentially-important atrial-fibroblast phenotype changes in CHF and emphasize the need for caution in relating properties of cultured fibroblasts to in vivo systems.