LRcell: detecting the source of differential expression at the sub-cell-type level from bulk RNA-seq data.

Given most tissues are consist of abundant and diverse (sub-)cell types, an important yet unaddressed problem in bulk RNA-seq analysis is to identify at which (sub-)cell type(s) the differential expression occurs. Single-cell RNA-sequencing (scRNA-seq) technologies can answer the question, but they are often labor-intensive and cost-prohibitive. Here, we present LRcell, a computational method aiming to identify specific (sub-)cell type(s) that drives the changes observed in a bulk RNA-seq experiment. In addition, LRcell provides pre-embedded marker genes computed from putative scRNA-seq experiments as options to execute the analyses. We conduct a simulation study to demonstrate the effectiveness and reliability of LRcell. Using three different real datasets, we show that LRcell successfully identifies known cell types involved in psychiatric disorders. Applying LRcell to bulk RNA-seq results can produce a hypothesis on which (sub-)cell type(s) contributes to the differential expression. LRcell is complementary to cell type deconvolution methods.

Disease category-specific annotation of variants using an ensemble learning framework.

Understanding the impact of non-coding sequence variants on complex diseases is an essential problem. We present a novel ensemble learning framework-CASAVA, to predict genomic loci in terms of disease category-specific risk. Using disease-associated variants identified by GWAS as training data, and diverse sequencing-based genomics and epigenomics profiles as features, CASAVA provides risk prediction of 24 major categories of diseases throughout the human genome. Our studies showed that CASAVA scores at a genomic locus provide a reasonable prediction of the disease-specific and disease category-specific risk prediction for non-coding variants located within the locus. Taking MHC2TA and immune system diseases as an example, we demonstrate the potential of CASAVA in revealing variant-disease associations. A website (http://zhanglabtools.org/CASAVA) has been built to facilitate easily access to CASAVA scores.

Role of TET1-mediated epigenetic modulation in Alzheimer's disease.

Alzheimer's disease (AD) is a neurodegenerative disorder influenced by a complex interplay of environmental, epigenetic, and genetic factors. DNA methylation (5mC) and hydroxymethylation (5hmC) are DNA modifications that serve as tissue-specific and temporal regulators of gene expression. TET family enzymes dynamically regulate these epigenetic modifications in response to environmental conditions, connecting environmental factors with gene expression. Previous epigenetic studies have identified 5mC and 5hmC changes associated with AD. In this study, we performed targeted resequencing of TET1 on a cohort of early-onset AD (EOAD) and control samples. Through gene-wise burden analysis, we observed significant enrichment of rare TET1 variants associated with AD (p = 0.04). We also profiled 5hmC in human postmortem brain tissues from AD and control groups. Our analysis identified differentially hydroxymethylated regions (DhMRs) in key genes responsible for regulating the methylome: TET3, DNMT3L, DNMT3A, and MECP2. To further investigate the role of Tet1 in AD pathogenesis, we used the 5xFAD mouse model with a Tet1 KO allele to examine how Tet1 loss influences AD pathogenesis. We observed significant changes in neuropathology, 5hmC, and RNA expression associated with Tet1 loss, while the behavioral alterations were not significant. The loss of Tet1 significantly increased amyloid plaque burden in the 5xFAD mouse (p = 0.044) and lead to a non-significant trend towards exacerbated AD-associated stress response in 5xFAD mice. At the molecular level, we found significant DhMRs enriched in genes involved in pathways responsible for neuronal projection organization, dendritic spine development and organization, and myelin assembly. RNA-Seq analysis revealed a significant increase in the expression of AD-associated genes such as Mpeg1, Ctsd, and Trem2. In conclusion, our results suggest that TET enzymes, particularly TET1, which regulate the methylome, may contribute to AD pathogenesis, as the loss of TET function increases AD-associated pathology.

Identification of androgen response-related lncRNAs in prostate cancer.

BACKGROUND: Long noncoding RNAs (lncRNAs) are RNA molecules with over 200 nucleotides that do not code for proteins, but are known to be widely expressed and have key roles in gene regulation and cellular functions. They are also found to be involved in the onset and development of various cancers, including prostate cancer (PCa). Since PCa are commonly driven by androgen regulated signaling, mainly stimulated pathways, identification and determining the influence of lncRNAs in androgen response is useful and necessary. LncRNAs regulated by the androgen receptor (AR) can serve as potential biomarkers for PCa. In the present study, gene expression data analysis were performed to distinguish lncRNAs related to the androgen response pathway. METHODS AND RESULTS: We used publicly available RNA-sequencing and ChIP-seq data to identify lncRNAs that are associated with the androgen response pathway. Using Universal Correlation Coefficient (UCC) and Pearson Correlation Coefficient (PCC) analyses, we found 15 lncRNAs that have (a) highly correlated expression with androgen response genes in PCa and are (b) differentially expressed in the setting of treatment with an androgen agonist as well as antagonist compared to controls. Using publicly available ChIP-seq data, we investigated the role of androgen/AR axis in regulating expression of these lncRNAs. We observed AR binding in the promoter regions of 5 lncRNAs (MIR99AHG, DUBR, DRAIC, PVT1, and COLCA1), showing the direct influence of AR on their expression and highlighting their association with the androgen response pathway. CONCLUSION: By utilizing publicly available multiomics data and by employing in silico methods, we identified five candidate lncRNAs that are involved in the androgen response pathway. These lncRNAs should be investigated as potential biomarkers for PCa.

Widespread rearrangement of 3D chromatin organization underlies polycomb-mediated stress-induced silencing.

Chromosomes of metazoan organisms are partitioned in the interphase nucleus into discrete topologically associating domains (TADs). Borders between TADs are formed in regions containing active genes and clusters of architectural protein binding sites. The transcription of most genes is repressed after temperature stress in Drosophila. Here we show that temperature stress induces relocalization of architectural proteins from TAD borders to inside TADs, and this is accompanied by a dramatic rearrangement in the 3D organization of the nucleus. TAD border strength declines, allowing for an increase in long-distance inter-TAD interactions. Similar but quantitatively weaker effects are observed upon inhibition of transcription or depletion of individual architectural proteins. Heat shock-induced inter-TAD interactions result in increased contacts among enhancers and promoters of silenced genes, which recruit Pc and form Pc bodies in the nucleolus. These results suggest that the TAD organization of metazoan genomes is plastic and can be reconfigured quickly.

Targeting MLL1 H3K4 methyltransferase activity in mixed-lineage leukemia.

Here we report a comprehensive characterization of our recently developed inhibitor MM-401 that targets the MLL1 H3K4 methyltransferase activity. MM-401 is able to specifically inhibit MLL1 activity by blocking MLL1-WDR5 interaction and thus the complex assembly. This targeting strategy does not affect other mixed-lineage leukemia (MLL) family histone methyltransferases (HMTs), revealing a unique regulatory feature for the MLL1 complex. Using MM-401 and its enantiomer control MM-NC-401, we show that inhibiting MLL1 methyltransferase activity specifically blocks proliferation of MLL cells by inducing cell-cycle arrest, apoptosis, and myeloid differentiation without general toxicity to normal bone marrow cells or non-MLL cells. More importantly, transcriptome analyses show that MM-401 induces changes in gene expression similar to those of MLL1 deletion, supporting a predominant role of MLL1 activity in regulating MLL1-dependent leukemia transcription program. We envision broad applications for MM-401 in basic and translational research.

Regulatory annotation of genomic intervals based on tissue-specific expression QTLs.

MOTIVATION: Annotating a given genomic locus or a set of genomic loci is an important yet challenging task. This is especially true for the non-coding part of the genome which is enormous yet poorly understood. Since gene set enrichment analyses have demonstrated to be effective approach to annotate a set of genes, the same idea can be extended to explore the enrichment of functional elements or features in a set of genomic intervals to reveal potential functional connections. RESULTS: In this study, we describe a novel computational strategy named loci2path that takes advantage of the newly emerged, genome-wide and tissue-specific expression quantitative trait loci (eQTL) information to help annotate a set of genomic intervals in terms of transcription regulation. By checking the presence or the absence of millions of eQTLs in a set of input genomic intervals, combined with grouping eQTLs by the pathways or gene sets that their target genes belong to, loci2path build a bridge connecting genomic intervals to functional pathways and pre-defined biological-meaningful gene sets, revealing potential for regulatory connection. Our method enjoys two key advantages over existing methods: first, we no longer rely on proximity to link a locus to a gene which has shown to be unreliable; second, eQTL allows us to provide the regulatory annotation under the context of specific tissue types. To demonstrate its utilities, we apply loci2path on sets of genomic intervals harboring disease-associated variants as query. Using 1 702 612 eQTLs discovered by the Genotype-Tissue Expression (GTEx) project across 44 tissues and 6320 pathways or gene sets cataloged in MSigDB as annotation resource, our method successfully identifies highly relevant biological pathways and revealed disease mechanisms for psoriasis and other immune-related diseases. Tissue specificity analysis of associated eQTLs provide additional evidence of the distinct roles of different tissues played in the disease mechanisms. AVAILABILITY AND IMPLEMENTATION: loci2path is published as an open source Bioconductor package, and it is available at http://bioconductor.org/packages/release/bioc/html/loci2path.html. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Application of topic models to a compendium of ChIP-Seq datasets uncovers recurrent transcriptional regulatory modules.

MOTIVATION: The availability of thousands of genome-wide coupling chromatin immunoprecipitation (ChIP)-Seq datasets across hundreds of transcription factors (TFs) and cell lines provides an unprecedented opportunity to jointly analyze large-scale TF-binding in vivo, making possible the discovery of the potential interaction and cooperation among different TFs. The interacted and cooperated TFs can potentially form a transcriptional regulatory module (TRM) (e.g. co-binding TFs), which helps decipher the combinatorial regulatory mechanisms. RESULTS: We develop a computational method tfLDA to apply state-of-the-art topic models to multiple ChIP-Seq datasets to decipher the combinatorial binding events of multiple TFs. tfLDA is able to learn high-order combinatorial binding patterns of TFs from multiple ChIP-Seq profiles, interpret and visualize the combinatorial patterns. We apply the tfLDA to two cell lines with a rich collection of TFs and identify combinatorial binding patterns that show well-known TRMs and related TF co-binding events. AVAILABILITY AND IMPLEMENTATION: A software R package tfLDA is freely available at https://github.com/lichen-lab/tfLDA. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Characterization of the EZH2-MMSET histone methyltransferase regulatory axis in cancer.

Histone methyltransferases (HMTases), as chromatin modifiers, regulate the transcriptomic landscape in normal development as well in diseases such as cancer. Here, we molecularly order two HMTases, EZH2 and MMSET, that have established genetic links to oncogenesis. EZH2, which mediates histone H3K27 trimethylation and is associated with gene silencing, was shown to be coordinately expressed and function upstream of MMSET, which mediates H3K36 dimethylation and is associated with active transcription. We found that the EZH2-MMSET HMTase axis is coordinated by a microRNA network and that the oncogenic functions of EZH2 require MMSET activity. Together, these results suggest that the EZH2-MMSET HMTase axis coordinately functions as a master regulator of transcriptional repression, activation, and oncogenesis and may represent an attractive therapeutic target in cancer.

RT States: systematic annotation of the human genome using cell type-specific replication timing programs.

MOTIVATION: The replication timing (RT) program has been linked to many key biological processes including cell fate commitment, 3D chromatin organization and transcription regulation. Significant technology progress now allows to characterize the RT program in the entire human genome in a high-throughput and high-resolution fashion. These experiments suggest that RT changes dynamically during development in coordination with gene activity. Since RT is such a fundamental biological process, we believe that an effective quantitative profile of the local RT program from a diverse set of cell types in various developmental stages and lineages can provide crucial biological insights for a genomic locus. RESULTS: In this study, we explored recurrent and spatially coherent combinatorial profiles from 42 RT programs collected from multiple lineages at diverse differentiation states. We found that a Hidden Markov Model with 15 hidden states provide a good model to describe these genome-wide RT profiling data. Each of the hidden state represents a unique combination of RT profiles across different cell types which we refer to as 'RT states'. To understand the biological properties of these RT states, we inspected their relationship with chromatin states, gene expression, functional annotation and 3D chromosomal organization. We found that the newly defined RT states possess interesting genome-wide functional properties that add complementary information to the existing annotation of the human genome. AVAILABILITY AND IMPLEMENTATION: R scripts for inferring HMM models and Perl scripts for further analysis are available https://github.com/PouletAxel/script_HMM_Replication_timing. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Therapeutic targeting of BET bromodomain proteins in castration-resistant prostate cancer.

Men who develop metastatic castration-resistant prostate cancer (CRPC) invariably succumb to the disease. Progression to CRPC after androgen ablation therapy is predominantly driven by deregulated androgen receptor (AR) signalling. Despite the success of recently approved therapies targeting AR signalling, such as abiraterone and second-generation anti-androgens including MDV3100 (also known as enzalutamide), durable responses are limited, presumably owing to acquired resistance. Recently, JQ1 and I-BET762 two selective small-molecule inhibitors that target the amino-terminal bromodomains of BRD4, have been shown to exhibit anti-proliferative effects in a range of malignancies. Here we show that AR-signalling-competent human CRPC cell lines are preferentially sensitive to bromodomain and extraterminal (BET) inhibition. BRD4 physically interacts with the N-terminal domain of AR and can be disrupted by JQ1 (refs 11, 13). Like the direct AR antagonist MDV3100, JQ1 disrupted AR recruitment to target gene loci. By contrast with MDV3100, JQ1 functions downstream of AR, and more potently abrogated BRD4 localization to AR target loci and AR-mediated gene transcription, including induction of the TMPRSS2-ERG gene fusion and its oncogenic activity. In vivo, BET bromodomain inhibition was more efficacious than direct AR antagonism in CRPC xenograft mouse models. Taken together, these studies provide a novel epigenetic approach for the concerted blockade of oncogenic drivers in advanced prostate cancer.

Gene density, transcription, and insulators contribute to the partition of the Drosophila genome into physical domains.

The mechanisms responsible for the establishment of physical domains in metazoan chromosomes are poorly understood. Here we find that physical domains in Drosophila chromosomes are demarcated at regions of active transcription and high gene density that are enriched for transcription factors and specific combinations of insulator proteins. Physical domains contain different types of chromatin defined by the presence of specific proteins and epigenetic marks, with active chromatin preferentially located at the borders and silenced chromatin in the interior. Domain boundaries participate in long-range interactions that may contribute to the clustering of regions of active or silenced chromatin in the nucleus. Analysis of transgenes suggests that chromatin is more accessible and permissive to transcription at the borders than inside domains, independent of the presence of active or silencing histone modifications. These results suggest that the higher-order physical organization of chromatin may impose an additional level of regulation over classical epigenetic marks.

FOXP3 orchestrates H4K16 acetylation and H3K4 trimethylation for activation of multiple genes by recruiting MOF and causing displacement of PLU-1.

Both H4K16 acetylation and H3K4 trimethylation are required for gene activation. However, it is still largely unclear how these modifications are orchestrated by transcriptional factors. Here, we analyzed the mechanism of the transcriptional activation by FOXP3, an X-linked suppressor of autoimmune diseases and cancers. FOXP3 binds near transcriptional start sites of its target genes. By recruiting MOF and displacing histone H3K4 demethylase PLU-1, FOXP3 increases both H4K16 acetylation and H3K4 trimethylation at the FOXP3-associated chromatins of multiple FOXP3-activated genes. RNAi-mediated silencing of MOF reduced both gene activation and tumor suppression by FOXP3, while both somatic mutations in clinical cancer samples and targeted mutation of FOXP3 in mouse prostate epithelial cells disrupted nuclear localization of MOF. Our data demonstrate a pull-push model in which a single transcription factor orchestrates two epigenetic alterations necessary for gene activation and provide a mechanism for somatic inactivation of the FOXP3 protein function in cancer cells.

Omicseq: a web-based search engine for exploring omics datasets.

The development and application of high-throughput genomics technologies has resulted in massive quantities of diverse omics data that continue to accumulate rapidly. These rich datasets offer unprecedented and exciting opportunities to address long standing questions in biomedical research. However, our ability to explore and query the content of diverse omics data is very limited. Existing dataset search tools rely almost exclusively on the metadata. A text-based query for gene name(s) does not work well on datasets wherein the vast majority of their content is numeric. To overcome this barrier, we have developed Omicseq, a novel web-based platform that facilitates the easy interrogation of omics datasets holistically to improve 'findability' of relevant data. The core component of Omicseq is trackRank, a novel algorithm for ranking omics datasets that fully uses the numerical content of the dataset to determine relevance to the query entity. The Omicseq system is supported by a scalable and elastic, NoSQL database that hosts a large collection of processed omics datasets. In the front end, a simple, web-based interface allows users to enter queries and instantly receive search results as a list of ranked datasets deemed to be the most relevant. Omicseq is freely available at http://www.omicseq.org.

Heat Shock Protein Beta-1 Modifies Anterior to Posterior Purkinje Cell Vulnerability in a Mouse Model of Niemann-Pick Type C Disease.

Selective neuronal vulnerability is characteristic of most degenerative disorders of the CNS, yet mechanisms underlying this phenomenon remain poorly characterized. Many forms of cerebellar degeneration exhibit an anterior-to-posterior gradient of Purkinje cell loss including Niemann-Pick type C1 (NPC) disease, a lysosomal storage disorder characterized by progressive neurological deficits that often begin in childhood. Here, we sought to identify candidate genes underlying vulnerability of Purkinje cells in anterior cerebellar lobules using data freely available in the Allen Brain Atlas. This approach led to the identification of 16 candidate neuroprotective or susceptibility genes. We demonstrate that one candidate gene, heat shock protein beta-1 (HSPB1), promoted neuronal survival in cellular models of NPC disease through a mechanism that involved inhibition of apoptosis. Additionally, we show that over-expression of wild type HSPB1 or a phosphomimetic mutant in NPC mice slowed the progression of motor impairment and diminished cerebellar Purkinje cell loss. We confirmed the modulatory effect of Hspb1 on Purkinje cell degeneration in vivo, as knockdown by Hspb1 shRNA significantly enhanced neuron loss. These results suggest that strategies to promote HSPB1 activity may slow the rate of cerebellar degeneration in NPC disease and highlight the use of bioinformatics tools to uncover pathways leading to neuronal protection in neurodegenerative disorders.

traseR: an R package for performing trait-associated SNP enrichment analysis in genomic intervals.

UNLABELLED: Genome-wide association studies (GWASs) have successfully identified many sequence variants that are significantly associated with common diseases and traits. Tens of thousands of such trait-associated SNPs have already been cataloged, which we believe form a great resource for genomic research. Recent studies have demonstrated that the collection of trait-associated SNPs can be exploited to indicate whether a given genomic interval or intervals are likely to be functionally connected with certain phenotypes or diseases. Despite this importance, currently, there is no ready-to-use computational tool able to connect genomic intervals to phenotypes. Here, we present traseR, an easy-to-use R Bioconductor package that performs enrichment analyses of trait-associated SNPs in arbitrary genomic intervals with flexible options, including testing method, type of background and inclusion of SNPs in LD. AVAILABILITY AND IMPLEMENTATION: The traseR R package preloaded with up-to-date collection of trait-associated SNPs are freely available in Bioconductor CONTACT: zhaohui.qin@emory.edu SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Bayesian inference with historical data-based informative priors improves detection of differentially expressed genes.

MOTIVATION: Modern high-throughput biotechnologies such as microarray are capable of producing a massive amount of information for each sample. However, in a typical high-throughput experiment, only limited number of samples were assayed, thus the classical 'large p, small n' problem. On the other hand, rapid propagation of these high-throughput technologies has resulted in a substantial collection of data, often carried out on the same platform and using the same protocol. It is highly desirable to utilize the existing data when performing analysis and inference on a new dataset. RESULTS: Utilizing existing data can be carried out in a straightforward fashion under the Bayesian framework in which the repository of historical data can be exploited to build informative priors and used in new data analysis. In this work, using microarray data, we investigate the feasibility and effectiveness of deriving informative priors from historical data and using them in the problem of detecting differentially expressed genes. Through simulation and real data analysis, we show that the proposed strategy significantly outperforms existing methods including the popular and state-of-the-art Bayesian hierarchical model-based approaches. Our work illustrates the feasibility and benefits of exploiting the increasingly available genomics big data in statistical inference and presents a promising practical strategy for dealing with the 'large p, small n' problem. AVAILABILITY AND IMPLEMENTATION: Our method is implemented in R package IPBT, which is freely available from https://github.com/benliemory/IPBT CONTACT: yuzhu@purdue.edu; zhaohui.qin@emory.edu SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Base-resolution methylation patterns accurately predict transcription factor bindings in vivo.

Detecting in vivo transcription factor (TF) binding is important for understanding gene regulatory circuitries. ChIP-seq is a powerful technique to empirically define TF binding in vivo. However, the multitude of distinct TFs makes genome-wide profiling for them all labor-intensive and costly. Algorithms for in silico prediction of TF binding have been developed, based mostly on histone modification or DNase I hypersensitivity data in conjunction with DNA motif and other genomic features. However, technical limitations of these methods prevent them from being applied broadly, especially in clinical settings. We conducted a comprehensive survey involving multiple cell lines, TFs, and methylation types and found that there are intimate relationships between TF binding and methylation level changes around the binding sites. Exploiting the connection between DNA methylation and TF binding, we proposed a novel supervised learning approach to predict TF-DNA interaction using data from base-resolution whole-genome methylation sequencing experiments. We devised beta-binomial models to characterize methylation data around TF binding sites and the background. Along with other static genomic features, we adopted a random forest framework to predict TF-DNA interaction. After conducting comprehensive tests, we saw that the proposed method accurately predicts TF binding and performs favorably versus competing methods.

A novel statistical method for quantitative comparison of multiple ChIP-seq datasets.

MOTIVATION: ChIP-seq is a powerful technology to measure the protein binding or histone modification strength in the whole genome scale. Although there are a number of methods available for single ChIP-seq data analysis (e.g. 'peak detection'), rigorous statistical method for quantitative comparison of multiple ChIP-seq datasets with the considerations of data from control experiment, signal to noise ratios, biological variations and multiple-factor experimental designs is under-developed. RESULTS: In this work, we develop a statistical method to perform quantitative comparison of multiple ChIP-seq datasets and detect genomic regions showing differential protein binding or histone modification. We first detect peaks from all datasets and then union them to form a single set of candidate regions. The read counts from IP experiment at the candidate regions are assumed to follow Poisson distribution. The underlying Poisson rates are modeled as an experiment-specific function of artifacts and biological signals. We then obtain the estimated biological signals and compare them through the hypothesis testing procedure in a linear model framework. Simulations and real data analyses demonstrate that the proposed method provides more accurate and robust results compared with existing ones. AVAILABILITY AND IMPLEMENTATION: An R software package ChIPComp is freely available at http://web1.sph.emory.edu/users/hwu30/software/ChIPComp.html.

One Size Doesn't Fit All - RefEditor: Building Personalized Diploid Reference Genome to Improve Read Mapping and Genotype Calling in Next Generation Sequencing Studies.

With rapid decline of the sequencing cost, researchers today rush to embrace whole genome sequencing (WGS), or whole exome sequencing (WES) approach as the next powerful tool for relating genetic variants to human diseases and phenotypes. A fundamental step in analyzing WGS and WES data is mapping short sequencing reads back to the reference genome. This is an important issue because incorrectly mapped reads affect the downstream variant discovery, genotype calling and association analysis. Although many read mapping algorithms have been developed, the majority of them uses the universal reference genome and do not take sequence variants into consideration. Given that genetic variants are ubiquitous, it is highly desirable if they can be factored into the read mapping procedure. In this work, we developed a novel strategy that utilizes genotypes obtained a priori to customize the universal haploid reference genome into a personalized diploid reference genome. The new strategy is implemented in a program named RefEditor. When applying RefEditor to real data, we achieved encouraging improvements in read mapping, variant discovery and genotype calling. Compared to standard approaches, RefEditor can significantly increase genotype calling consistency (from 43% to 61% at 4X coverage; from 82% to 92% at 20X coverage) and reduce Mendelian inconsistency across various sequencing depths. Because many WGS and WES studies are conducted on cohorts that have been genotyped using array-based genotyping platforms previously or concurrently, we believe the proposed strategy will be of high value in practice, which can also be applied to the scenario where multiple NGS experiments are conducted on the same cohort. The RefEditor sources are available at https://github.com/superyuan/refeditor.