Characterizing efficient feature selection for single-cell expression analysis.

Unsupervised feature selection is a critical step for efficient and accurate analysis of single-cell RNA-seq data. Previous benchmarks used two different criteria to compare feature selection methods: (i) proportion of ground-truth marker genes included in the selected features and (ii) accuracy of cell clustering using ground-truth cell types. Here, we systematically compare the performance of 11 feature selection methods for both criteria. We first demonstrate the discordance between these criteria and suggest using the latter. We then compare the distribution of selected genes in their means between feature selection methods. We show that lowly expressed genes exhibit seriously high coefficients of variation and are mostly excluded by high-performance methods. In particular, high-deviation- and high-expression-based methods outperform the widely used in Seurat package in clustering cells and data visualization. We further show they also enable a clear separation of the same cell type from different tissues as well as accurate estimation of cell trajectories.

netGO: R-Shiny package for network-integrated pathway enrichment analysis.

SUMMARY: We present an R-Shiny package, netGO, for novel network-integrated pathway enrichment analysis. The conventional Fisher's exact test (FET) considers the extent of overlap between target genes and pathway gene-sets, while recent network-based analysis tools consider only network interactions between the two. netGO implements an intuitive framework to integrate both the overlap and networks into a single score, and adaptively resamples genes based on network degrees to assess the pathway enrichment. In benchmark tests for gene expression and genome-wide association study (GWAS) data, netGO captured the relevant gene-sets better than existing tools, especially when analyzing a small number of genes. Specifically, netGO provides user-interactive visualization of the target genes, enriched gene-set and their network interactions for both netGO and FET results for further analysis. For this visualization, we also developed a standalone R-Shiny package shinyCyJS to connect R-shiny and the JavaScript version of cytoscape. AVAILABILITY AND IMPLEMENTATION: netGO R-Shiny package is freely available from github, https://github.com/unistbig/netGO. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Glucosylceramide synthase regulates adipo-osteogenic differentiation through synergistic activation of PPARγ with GlcCer.

Dysregulation of the adipo-osteogenic differentiation balance of mesenchymal stem cells (MSCs), which are common progenitor cells of adipocytes and osteoblasts, has been associated with many pathophysiologic diseases, such as obesity, osteopenia, and osteoporosis. Growing evidence suggests that lipid metabolism is crucial for maintaining stem cell homeostasis and cell differentiation; however, the detailed underlying mechanisms are largely unknown. Here, we demonstrate that glucosylceramide (GlcCer) and its synthase, glucosylceramide synthase (GCS), are key determinants of MSC differentiation into adipocytes or osteoblasts. GCS expression was increased during adipogenesis and decreased during osteogenesis. Targeting GCS using RNA interference or a chemical inhibitor enhanced osteogenesis and inhibited adipogenesis by controlling the transcriptional activity of peroxisome proliferator-activated receptor γ (PPARγ). Treatment with GlcCer sufficiently rescued adipogenesis and inhibited osteogenesis in GCS knockdown MSCs. Mechanistically, GlcCer interacted directly with PPARγ through A/B domain and synergistically enhanced rosiglitazone-induced PPARγ activation without changing PPARγ expression, thereby treatment with exogenous GlcCer increased adipogenesis and inhibited osteogenesis. Animal studies demonstrated that inhibiting GCS reduced adipocyte formation in white adipose tissues under normal chow diet and high-fat diet feeding and accelerated bone repair in a calvarial defect model. Taken together, our findings identify a novel lipid metabolic regulator for the control of MSC differentiation and may have important therapeutic implications.

Diffusible Signaling Factor, a Quorum-Sensing Molecule, Interferes with and Is Toxic Towards Bdellovibrio bacteriovorus 109J.

Bdellovibrio bacteriovorus 109J is a predatory bacterium which lives by predating on other Gram-negative bacteria to obtain the nutrients it needs for replication and survival. Here, we evaluated the effects two classes of bacterial signaling molecules (acyl homoserine lactones (AHLs) and diffusible signaling factor (DSF)) have on B. bacteriovorus 109J behavior and viability. While AHLs had a non-significant impact on predation rates, DSF considerably delayed predation and bdelloplast lysis. Subsequent experiments showed that 50 µM DSF also reduced the motility of attack-phase B. bacteriovorus 109J cells by 50% (38.2 ± 14.9 vs. 17 ± 8.9 µm/s). Transcriptomic analyses found that DSF caused genome-wide changes in B. bacteriovorus 109J gene expression patterns during both the attack and intraperiplasmic phases, including the significant downregulation of the flagellum assembly genes and numerous serine protease genes. While the former accounts for the reduced speeds observed, the latter was confirmed experimentally with 50 µM DSF completely blocking protease secretion from attack-phase cells. Additional experiments found that 30% of the total cellular ATP was released into the supernatant when B. bacteriovorus 109J was exposed to 200 µM DSF, implying that this QS molecule negatively impacts membrane integrity.

Biclustering analysis of transcriptome big data identifies condition-specific microRNA targets.

We present a novel approach to identify human microRNA (miRNA) regulatory modules (mRNA targets and relevant cell conditions) by biclustering a large collection of mRNA fold-change data for sequence-specific targets. Bicluster targets were assessed using validated messenger RNA (mRNA) targets and exhibited on an average 17.0% (median 19.4%) improved gain in certainty (sensitivity + specificity). The net gain was further increased up to 32.0% (median 33.4%) by incorporating functional networks of targets. We analyzed cancer-specific biclusters and found that the PI3K/Akt signaling pathway is strongly enriched with targets of a few miRNAs in breast cancer and diffuse large B-cell lymphoma. Indeed, five independent prognostic miRNAs were identified, and repression of bicluster targets and pathway activity by miR-29 was experimentally validated. In total, 29 898 biclusters for 459 human miRNAs were collected in the BiMIR database where biclusters are searchable for miRNAs, tissues, diseases, keywords and target genes.

Efficient pathway enrichment and network analysis of GWAS summary data using GSA-SNP2.

Pathway-based analysis in genome-wide association study (GWAS) is being widely used to uncover novel multi-genic functional associations. Many of these pathway-based methods have been used to test the enrichment of the associated genes in the pathways, but exhibited low powers and were highly affected by free parameters. We present the novel method and software GSA-SNP2 for pathway enrichment analysis of GWAS P-value data. GSA-SNP2 provides high power, decent type I error control and fast computation by incorporating the random set model and SNP-count adjusted gene score. In a comparative study using simulated and real GWAS data, GSA-SNP2 exhibited high power and best prioritized gold standard positive pathways compared with six existing enrichment-based methods and two self-contained methods (alternative pathway analysis approach). Based on these results, the difference between pathway analysis approaches was investigated and the effects of the gene correlation structures on the pathway enrichment analysis were also discussed. In addition, GSA-SNP2 is able to visualize protein interaction networks within and across the significant pathways so that the user can prioritize the core subnetworks for further studies. GSA-SNP2 is freely available at https://sourceforge.net/projects/gsasnp2.

GScluster: network-weighted gene-set clustering analysis.

BACKGROUND: Gene-set analysis (GSA) has been commonly used to identify significantly altered pathways or functions from omics data. However, GSA often yields a long list of gene-sets, necessitating efficient post-processing for improved interpretation. Existing methods cluster the gene-sets based on the extent of their overlap to summarize the GSA results without considering interactions between gene-sets. RESULTS: Here, we presented a novel network-weighted gene-set clustering that incorporates both the gene-set overlap and protein-protein interaction (PPI) networks. Three examples were demonstrated for microarray gene expression, GWAS summary, and RNA-sequencing data to which different GSA methods were applied. These examples as well as a global analysis show that the proposed method increases PPI densities and functional relevance of the resulting clusters. Additionally, distinct properties of gene-set distance measures were compared. The methods are implemented as an R/Shiny package GScluster that provides gene-set clustering and diverse functions for visualization of gene-sets and PPI networks. CONCLUSIONS: Network-weighted gene-set clustering provides functionally more relevant gene-set clusters and related network analysis.

Gene dispersion is the key determinant of the read count bias in differential expression analysis of RNA-seq data.

BACKGROUND: In differential expression analysis of RNA-sequencing (RNA-seq) read count data for two sample groups, it is known that highly expressed genes (or longer genes) are more likely to be differentially expressed which is called read count bias (or gene length bias). This bias had great effect on the downstream Gene Ontology over-representation analysis. However, such a bias has not been systematically analyzed for different replicate types of RNA-seq data. RESULTS: We show that the dispersion coefficient of a gene in the negative binomial modeling of read counts is the critical determinant of the read count bias (and gene length bias) by mathematical inference and tests for a number of simulated and real RNA-seq datasets. We demonstrate that the read count bias is mostly confined to data with small gene dispersions (e.g., technical replicates and some of genetically identical replicates such as cell lines or inbred animals), and many biological replicate data from unrelated samples do not suffer from such a bias except for genes with some small counts. It is also shown that the sample-permuting GSEA method yields a considerable number of false positives caused by the read count bias, while the preranked method does not. CONCLUSION: We showed the small gene variance (similarly, dispersion) is the main cause of read count bias (and gene length bias) for the first time and analyzed the read count bias for different replicate types of RNA-seq data and its effect on gene-set enrichment analysis.

Indole negatively impacts predation by Bdellovibrio bacteriovorus and its release from the bdelloplast.

Bdellovibrio bacteriovorus is a predatory bacterium that attacks a wide range of Gram-negative bacterial pathogens and is proposed to be a potential living antibiotic. In this study, we evaluated the effects of indole, a bacterial signalling molecule commonly produced within the gut, on the predatory ability of B. bacteriovorus HD100. Indole significantly delayed predation on Escherichia coli MG1655 and Salmonella enterica KACC 11595 at physiological concentrations (0.25 to 1 mM) and completely inhibited predation when present at 2 mM. Microscopic analysis revealed that indole blocked the predator from attacking the prey. Furthermore, indole was not toxic to the predator but slowed down its motility. Microarray and reverse transcription quantitative polymerase chain reaction (RT-qPCR) analyses confirmed that as the gene group showing the greatest downregulation in the presence of indole was flagellar assembly genes. Indole also caused a wide spectrum changes in gene expression including general downregulation of genes involved in ribosome assembly. Furthermore, indole addition to the predatory culture after the entrance of B. bacteriovorus into the prey periplasm slowed down bdelloplast lysis. In conclusion, indole can have significant impacts on the predation efficiency, which should be taken into consideration especially if B. bacteriovorus is to be applied as a probiotic or living antibiotic.

REGNET: mining context-specific human transcription networks using composite genomic information.

BACKGROUND: Genome-wide expression profiles reflect the transcriptional networks specific to the given cell context. However, most statistical models try to estimate the average connectivity of the networks from a collection of gene expression data, and are unable to characterize the context-specific transcriptional regulations. We propose an approach for mining context-specific transcription networks from a large collection of gene expression fold-change profiles and composite gene-set information. RESULTS: Using a composite gene-set analysis method, we combine the information of transcription factor binding sites, Gene Ontology or pathway gene sets and gene expression fold-change profiles for a variety of cell conditions. We then collected all the significant patterns and constructed a database of context-specific transcription networks for human (REGNET). As a result, context-specific roles of transcription factors as well as their functional targets are readily explored. To validate the approach, nine predicted targets of E2F1 in HeLa cells were tested using chromatin immunoprecipitation assay. Among them, five (Gadd45b, Dusp6, Mll5, Bmp2 and E2f3) were successfully bound by E2F1. c-JUN and the EMT transcription networks were also validated from literature. CONCLUSIONS: REGNET is a useful tool for exploring the ternary relationships among the transcription factors, their functional targets and the corresponding cell conditions. It is able to provide useful clues for novel cell-specific transcriptional regulations. The REGNET database is available at http://mgrc.kribb.re.kr/regnet.

WegoLoc: accurate prediction of protein subcellular localization using weighted Gene Ontology terms.

SUMMARY: We present an accurate and fast web server, WegoLoc for predicting subcellular localization of proteins based on sequence similarity and weighted Gene Ontology (GO) information. A term weighting method in the text categorization process is applied to GO terms for a support vector machine classifier. As a result, WegoLoc surpasses the state-of-the-art methods for previously used test datasets. WegoLoc supports three eukaryotic kingdoms (animals, fungi and plants) and provides human-specific analysis, and covers several sets of cellular locations. In addition, WegoLoc provides (i) multiple possible localizations of input protein(s) as well as their corresponding probability scores, (ii) weights of GO terms representing the contribution of each GO term in the prediction, and (iii) a BLAST E-value for the best hit with GO terms. If the similarity score does not meet a given threshold, an amino acid composition-based prediction is applied as a backup method. AVAILABILITY: WegoLoc and User's guide are freely available at the website http://www.btool.org/WegoLoc CONTACT: smchiks@ks.ac.kr; dougnam@unist.ac.kr SUPPLEMENTARY INFORMATION: Supplementary data is available at http://www.btool.org/WegoLoc.

ADGO 2.0: interpreting microarray data and list of genes using composite annotations.

ADGO 2.0 is a web-based tool that provides composite interpretations for microarray data comparing two sample groups as well as lists of genes from diverse sources of biological information. Some other tools also incorporate composite annotations solely for interpreting lists of genes but usually provide highly redundant information. This new version has the following additional features: first, it provides multiple gene set analysis methods for microarray inputs as well as enrichment analyses for lists of genes. Second, it screens redundant composite annotations when generating and prioritizing them. Third, it incorporates union and subtracted sets as well as intersection sets. Lastly, users can upload their own gene sets (e.g. predicted miRNA targets) to generate and analyze new composite sets. The first two features are unique to ADGO 2.0. Using our tool, we demonstrate analyses of a microarray dataset and a list of genes for T-cell differentiation. The new ADGO is available at http://www.btool.org/ADGO2.

Benchmarking integration of single-cell differential expression.

Integration of single-cell RNA sequencing data between different samples has been a major challenge for analyzing cell populations. However, strategies to integrate differential expression analysis of single-cell data remain underinvestigated. Here, we benchmark 46 workflows for differential expression analysis of single-cell data with multiple batches. We show that batch effects, sequencing depth and data sparsity substantially impact their performances. Notably, we find that the use of batch-corrected data rarely improves the analysis for sparse data, whereas batch covariate modeling improves the analysis for substantial batch effects. We show that for low depth data, single-cell techniques based on zero-inflation model deteriorate the performance, whereas the analysis of uncorrected data using limmatrend, Wilcoxon test and fixed effects model performs well. We suggest several high-performance methods under different conditions based on various simulation and real data analyses. Additionally, we demonstrate that differential expression analysis for a specific cell type outperforms that of large-scale bulk sample data in prioritizing disease-related genes.

GSA-SNP: a general approach for gene set analysis of polymorphisms.

Genome-wide association (GWA) study aims to identify the genetic factors associated with the traits of interest. However, the power of GWA analysis has been seriously limited by the enormous number of markers tested. Recently, the gene set analysis (GSA) methods were introduced to GWA studies to address the association of gene sets that share common biological functions. GSA considerably increased the power of association analysis and successfully identified coordinated association patterns of gene sets. There have been several approaches in this direction with some limitations. Here, we present a general approach for GSA in GWA analysis and a stand-alone software GSA-SNP that implements three widely used GSA methods. GSA-SNP provides a fast computation and an easy-to-use interface. The software and test datasets are freely available at http://gsa.muldas.org. We provide an exemplary analysis on adult heights in a Korean population.

MicroRNA-29 Ameliorates Fibro-Inflammation and Insulin Resistance in HIF1α-Deficient Obese Adipose Tissue by Inhibiting Endotrophin Generation.

Dysregulation of extracellular matrix proteins in obese adipose tissue (AT) induces systemic insulin resistance. The metabolic roles of type VI collagen and its cleavage peptide endotrophin in obese AT are well established. However, the mechanisms regulating endotrophin generation remain elusive. Herein, we identified that several endotrophin-containing peptides (pre-endotrophins) were generated from the COL6A3 chain in a stepwise manner for the efficient production of mature endotrophin, partly through the action of hypoxia-induced matrix metalloproteinases (MMPs), including MMP2, MMP9, and MMP16. Hypoxia is an upstream regulator of COL6A3 expression and the proteolytic processing that regulates endotrophin generation. Hypoxia-inducible factor 1α (HIF1α) and the hypoxia-associated suppression of microRNA-29 (miR-29) cooperatively control the levels of COL6A3 and MMPs, which are responsible for endotrophin generation in hypoxic ATs. Adipocyte-specific Hif1α knock-out (APN-HIF1αKO) mice fed a chronic high-fat diet exhibited the significant amelioration of both local fibro-inflammation in AT and systemic insulin resistance compared with their control littermates, partly through the inhibition of endotrophin generation. Strikingly, adenovirus-mediated miR-29 overexpression in the ATs of APN-HIF1αKO mice in obesity significantly decreased endotrophin levels, suggesting that miR-29, combined with HIF1α inhibition in AT, could be a promising therapeutic strategy for treating obesity and related metabolic diseases.

De-correlating expression in gene-set analysis.

MOTIVATION: Group-wise pattern analysis of genes, known as gene-set analysis (GSA), addresses the differential expression pattern of biologically pre-defined gene sets. GSA exhibits high statistical power and has revealed many novel biological processes associated with specific phenotypes. In most cases, however, GSA relies on the invalid assumption that the members of each gene set are sampled independently, which increases false predictions. RESULTS: We propose an algorithm, termed DECO, to remove (or alleviate) the bias caused by the correlation of the expression data in GSAs. This is accomplished through the eigenvalue-decomposition of covariance matrixes and a series of linear transformations of data. In particular, moderate de-correlation methods that truncate or re-scale eigenvalues were proposed for a more reliable analysis. Tests of simulated and real experimental data show that DECO effectively corrects the correlation structure of gene expression and improves the prediction accuracy (specificity and sensitivity) for both gene- and sample-randomizing GSA methods. AVAILABILITY: The MATLAB codes and the tested data sets are available at ftp://deco.nims.re.kr/pub or from the author.

Powerful p-value combination methods to detect incomplete association.

Meta-analyses increase statistical power by combining statistics from multiple studies. Meta-analysis methods have mostly been evaluated under the condition that all the data in each study have an association with the given phenotype. However, specific experimental conditions in each study or genetic heterogeneity can result in "unassociated statistics" that are derived from the null distribution. Here, we show that power of conventional meta-analysis methods rapidly decreases as an increasing number of unassociated statistics are included, whereas the classical Fisher's method and its weighted variant (wFisher) exhibit relatively high power that is robust to addition of unassociated statistics. We also propose another robust method based on joint distribution of ordered p-values (ordmeta). Simulation analyses for t-test, RNA-seq, and microarray data demonstrated that wFisher and ordmeta, when only a small number of studies have an association, outperformed existing meta-analysis methods. We performed meta-analyses of nine microarray datasets (prostate cancer) and four association summary datasets (body mass index), where our methods exhibited high biological relevance and were able to detect genes that the-state-of-the-art methods missed. The metapro R package that implements the proposed methods is available from both CRAN and GitHub ( http://github.com/unistbig/metapro ).

Hepatic MIR20B promotes nonalcoholic fatty liver disease by suppressing PPARA.

Background: Non-alcoholic fatty liver disease (NAFLD) is characterized by excessive lipid accumulation and imbalances in lipid metabolism in the liver. Although nuclear receptors (NRs) play a crucial role in hepatic lipid metabolism, the underlying mechanisms of NR regulation in NAFLD remain largely unclear. Methods: Using network analysis and RNA-seq to determine the correlation between NRs and microRNA in human NAFLD patients, we revealed that MIR20B specifically targets PPARA. MIR20B mimic and anti-MIR20B were administered to human HepG2 and Huh-7 cells and mouse primary hepatocytes as well as high-fat diet (HFD)- or methionine-deficient diet (MCD)-fed mice to verify the specific function of MIR20B in NAFLD. We tested the inhibition of the therapeutic effect of a PPARα agonist, fenofibrate, by Mir20b and the synergic effect of combination of fenofibrate with anti-Mir20b in NAFLD mouse model. Results: We revealed that MIR20B specifically targets PPARA through miRNA regulatory network analysis of nuclear receptor genes in NAFLD. The expression of MIR20B was upregulated in free fatty acid (FA)-treated hepatocytes and the livers of both obesity-induced mice and NAFLD patients. Overexpression of MIR20B significantly increased hepatic lipid accumulation and triglyceride levels. Furthermore, MIR20B significantly reduced FA oxidation and mitochondrial biogenesis by targeting PPARA. In Mir20b-introduced mice, the effect of fenofibrate to ameliorate hepatic steatosis was significantly suppressed. Finally, inhibition of Mir20b significantly increased FA oxidation and uptake, resulting in improved insulin sensitivity and a decrease in NAFLD progression. Moreover, combination of fenofibrate and anti-Mir20b exhibited the synergic effect on improvement of NAFLD in MCD-fed mice. Conclusions: Taken together, our results demonstrate that the novel MIR20B targets PPARA, plays a significant role in hepatic lipid metabolism, and present an opportunity for the development of novel therapeutics for NAFLD. Funding: This research was funded by Korea Mouse Phenotyping Project (2016M3A9D5A01952411), the National Research Foundation of Korea (NRF) grant funded by the Korea government (2020R1F1A1061267, 2018R1A5A1024340, NRF-2021R1I1A2041463, 2020R1I1A1A01074940, 2016M3C9A394589324), and the Future-leading Project Research Fund (1.210034.01) of UNIST.

Gene-set approach for expression pattern analysis.

Recently developed gene set analysis methods evaluate differential expression patterns of gene groups instead of those of individual genes. This approach especially targets gene groups whose constituents show subtle but coordinated expression changes, which might not be detected by the usual individual gene analysis. The approach has been quite successful in deriving new information from expression data, and a number of methods and tools have been developed intensively in recent years. We review those methods and currently available tools, classify them according to the statistical methods employed, and discuss their pros and cons. We also discuss several interesting extensions to the methods.

Benchmarking RNA-seq differential expression analysis methods using spike-in and simulation data.

Benchmarking RNA-seq differential expression analysis methods using spike-in and simulated RNA-seq data has often yielded inconsistent results. The spike-in data, which were generated from the same bulk RNA sample, only represent technical variability, making the test results less reliable. We compared the performance of 12 differential expression analysis methods for RNA-seq data, including recent variants in widely used software packages, using both RNA spike-in and simulation data for negative binomial (NB) model. Performance of edgeR, DESeq2, and ROTS was particularly different between the two benchmark tests. Then, each method was tested under most extensive simulation conditions especially demonstrating the large impacts of proportion, dispersion, and balance of differentially expressed (DE) genes. DESeq2, a robust version of edgeR (edgeR.rb), voom with TMM normalization (voom.tmm) and sample weights (voom.sw) showed an overall good performance regardless of presence of outliers and proportion of DE genes. The performance of RNA-seq DE gene analysis methods substantially depended on the benchmark used. Based on the simulation results, suitable methods were suggested under various test conditions.