recount workflow: Accessing over 70,000 human RNA-seq samples with Bioconductor

Luminal epithelial cells integrate variable responses to aging into stereotypical changes that underlie breast cancer susceptibility

Effects from aging in single cells are heterogenous, whereas at the organ- and tissue-levels aging phenotypes tend to appear as stereotypical changes. The mammary epithelium is a bilayer of two major phenotypically and functionally distinct cell lineages: luminal epithelial and myoepithelial cells. Mammary luminal epithelia exhibit substantial stereotypical changes with age that merit attention because these cells are the putative cells-of-origin for breast cancers. We hypothesize that effects from aging that impinge upon maintenance of lineage fidelity increase susceptibility to cancer initiation. We generated and analyzed transcriptomes from primary luminal epithelial and myoepithelial cells from younger <30 (y)ears old and older >55 y women. In addition to age-dependent directional changes in gene expression, we observed increased transcriptional variance with age that contributed to genome-wide loss of lineage fidelity. Age-dependent variant responses were common to both lineages, whereas directional changes were almost exclusively detected in luminal epithelia and involved altered regulation of chromatin and genome organizers such as SATB1. Epithelial expression variance of gap junction protein GJB6 increased with age, and modulation of GJB6 expression in heterochronous co-cultures revealed that it provided a communication conduit from myoepithelial cells that drove directional change in luminal cells. Age-dependent luminal transcriptomes comprised a prominent signal that could be detected in bulk tissue during aging and transition into cancers. A machine learning classifier based on luminal-specific aging distinguished normal from cancer tissue and was highly predictive of breast cancer subtype. We speculate that luminal epithelia are the ultimate site of integration of the variant responses to aging in their surrounding tissue, and that their emergent phenotype both endows cells with the ability to become cancer-cells-of-origin and represents a biosensor that presages cancer susceptibility.

Genotype prediction of 336,463 samples from public expression data

Tens of thousands of RNA-sequencing experiments comprising hundreds of thousands of individual samples have now been performed. These data represent a broad range of experimental conditions, sequencing technologies, and hypotheses under study. The Recount project has aggregated and uniformly processed hundreds of thousands of publicly available RNA-seq samples. Most of these samples only include RNA expression measurements; genotype data for these same samples would enable a wide range of analyses including variant prioritization, eQTL analysis, and studies of allele specific expression. Here, we developed a statistical model based on the existing reference and alternative read counts from the RNA-seq experiments available through Recount3 to predict genotypes at autosomal biallelic loci in coding regions. We demonstrate the accuracy of our model using large-scale studies that measured both gene expression and genotype genome-wide. We show that our predictive model is highly accurate with 99.5% overall accuracy, 99.6% major allele accuracy, and 90.4% minor allele accuracy. Our model is robust to tissue and study effects, provided the coverage is high enough. We applied this model to genotype all the samples in Recount 3 and provide the largest ready-to-use expression repository containing genotype information. We illustrate that the predicted genotype from RNA-seq data is sufficient to unravel the underlying population structure of samples in Recount3 using Principal Component Analysis.

Comprehensive Analysis on Prognostic Signature Based on T Cell-Mediated Tumor Killing Related Genes in Gastric Cancer

Although immunotherapy is a valuable treatment for gastric cancer (GC), identifying the patients who would benefit most from this approach presents a challenge. In this study, GC patients were divided into two subtypes by consensus clustering according to T cell-mediated tumor killing related genes (TTKRGs), and there were significant differences in tumor-infiltrating immune cells, signaling pathways, and gene expression of immunomodulators and inhibitory immune checkpoints between the two subtypes. Then, we developed an individualized signature based on TTKRGs, and its clinical and predictive value in GC patients for chemotherapeutic and immunotherapeutic responses was assessed. We confirmed the expression levels of signature genes in GC tumor tissue using quantitative real-time polymerase chain reaction (qRT-PCR). Additionally, to improve the accuracy of GC prognosis predictions, we established a nomogram. We further identified some compounds as sensitive drugs targeting GC risk groups. The signature showed significant predictive ability across RNA-seq, microarray, and qRT-PCR cohorts, which could assist in predicting survival, immunotherapeutic and chemotherapeutic outcomes in GC patients.

Deep Learning for Predicting Gene Regulatory Networks: A Step-by-Step Protocol in R

Deep learning has emerged as a powerful tool for solving complex problems, including reconstruction of gene regulatory networks within the realm of biology. These networks consist of transcription factors and their associations with genes they regulate. Despite the utility of deep learning methods in studying gene expression and regulation, their accessibility remains limited for biologists,  mainly due to the prerequisites of programming skills and a nuanced grasp of the underlying algorithms. This chapter presents a deep learning protocol that utilize TensorFlow and the Keras API in R/RStudio, with the aim of making deep learning accessible for individuals without specialized expertise. The protocol focuses on the genome-wide prediction of regulatory interactions between transcription factors and genes, leveraging publicly available gene expression data in conjunction with well-established benchmarks. The protocol encompasses pivotal phases including data preprocessing, conceptualization of neural network architectures, iterative processes of model training and validation, as well as forecasting of novel regulatory associations. Furthermore, it provides insights into parameter tuning for deep learning models. By adhering to this protocol, researchers are expected to gain a comprehensive understanding of applying deep learning techniques to predict regulatory interactions. This protocol can be readily modifiable to serve diverse research problems, thereby empowering scientists to effectively harness the capabilities of deep learning in their investigations.

Heterogeneity in the gene regulatory landscape of leiomyosarcoma

Characterizing inter-tumor heterogeneity is crucial for selecting suitable cancer therapy, as the presence of diverse molecular subgroups of patients can be associated with disease outcome or response to treatment. While cancer subtypes are often characterized by differences in gene expression, the mechanisms driving these differences are generally unknown. We set out to model the regulatory mechanisms driving sarcoma heterogeneity based on patient-specific, genome-wide gene regulatory networks. We developed a new computational framework, PORCUPINE, which combines knowledge on biological pathways with permutation-based network analysis to identify pathways that exhibit significant regulatory heterogeneity across a patient population. We applied PORCUPINE to patient-specific leiomyosarcoma networks modeled on data from The Cancer Genome Atlas and validated our results in an independent dataset from the German Cancer Research Center. PORCUPINE identified 37 heterogeneously regulated pathways, including pathways representing potential targets for treatment of subgroups of leiomyosarcoma patients, such as FGFR and CTLA4 inhibitory signaling. We validated the detected regulatory heterogeneity through analysis of networks and chromatin states in leiomyosarcoma cell lines. We showed that the heterogeneity identified with PORCUPINE is not associated with methylation profiles or clinical features, thereby suggesting an independent mechanism of patient heterogeneity driven by the complex landscape of gene regulatory interactions.

Comprehensive analysis of ubiquitously expressed genes in human, from a data-driven perspective

Comprehensive characterization of spatial and temporal gene expression patterns in humans is critical for uncovering the regulatory codes of the human genome and understanding the molecular mechanism of human disease. The ubiquitously expressed genes (UEGs) refer to those genes expressed across a majority, if not all, phenotypic and physiological conditions of an organism. It is known that many human genes are broadly expressed across tissues. However, most previous UEG studies have only focused on providing a list of UEGs without capturing their global expression patterns, thus limiting the potential use of UEG information. In this article, we proposed a novel data-driven framework to leverage the extensive collection of ∼40,000 human transcriptomes to derive a list of UEGs and their corresponding global expression patterns, which offers a valuable resource to further characterize human transcriptome. Our results suggest that about half (12,234; 49.01%) of the human genes are expressed in at least 80% of human transcriptomes, and the median size of the human transcriptome is 16,342 (65.44%). Through gene clustering, we identified a set of UEGs, named LoVarUEGs, that have stable expression across human transcriptomes and can be used as internal reference genes for expression measurement. To further demonstrate the usefulness of this resource, we evaluated the global expression patterns for 16 previously predicted disallowed genes in islets beta cells and found that seven of these genes showed relatively more varied expression patterns, suggesting that the repression of these genes may not be unique to islets beta cells.

Novel integrated workflow allows production and in-depth quality assessment of multifactorial reprogrammed skeletal muscle cells from human stem cells

Skeletal muscle tissue engineering aims at generating biological substitutes that restore, maintain or improve normal muscle function; however, the quality of cells produced by current protocols remains insufficient. Here, we developed a multifactor-based protocol that combines adenovector (AdV)-mediated MYOD expression, small molecule inhibitor and growth factor treatment, and electrical pulse stimulation (EPS) to efficiently reprogram different types of human-derived multipotent stem cells into physiologically functional skeletal muscle cells (SMCs). The protocol was complemented through a novel in silico workflow that allows for in-depth estimation and potentially optimization of the quality of generated muscle tissue, based on the transcriptomes of transdifferentiated cells. We additionally patch-clamped phenotypic SMCs to associate their bioelectrical characteristics with their transcriptome reprogramming. Overall, we set up a comprehensive and dynamic approach at the nexus of viral vector-based technology, bioinformatics, and electrophysiology that facilitates production of high-quality skeletal muscle cells and can guide iterative cycles to improve myo-differentiation protocols.

Comprehensive Analysis of an Individualized Immune-Related lncRNA Pair Signature in Gastric Cancer

Long noncoding RNAs (lncRNAs) have diverse functions, including immune regulation. Increasing studies have reported immune-related lncRNAs in the prognosis of multiple cancers. In this study, we developed an individualized signature containing 13 immune-related lncRNA pairs (IRLPs) which could predict the overall survival, disease-free survival, progression-free survival, and disease-specific survival of gastric cancer (GC) patients in The Cancer Genome Atlas (TCGA) cohort, and internal and external validations, signature comparisons, and subgroup analyses further confirmed its superiority, stability, and generalizability. Notably, this signature also showed good applicability in discriminating the prognosis of pan-cancer patients. Then, we constructed and validated a nomogram for overall survival based on the signature and clinical factors, which allowed more accurate predictions of GC prognosis. In addition, we revealed that the low survival rate of patients with high-risk scores may be due to their aggressive clinical features, enriched cancer-related signaling pathways, the infiltration of specific immunosuppressive cells, and low tumor mutation burden. We further predicted obviously worse immunotherapeutic responses in the high-risk groups and identified some candidate compounds targeting GC risk group differentiation. This signature based on the IRLPs may be promising for predicting the survival outcomes and immunotherapeutic responses of GC patients in clinical practice.

A Multiparametric Pharmacogenomic Strategy for Drug Repositioning predicts Therapeutic Efficacy for Glioblastoma Cell Lines

Background

Poor prognosis of glioblastoma patients and the extensive heterogeneity of glioblastoma at both the molecular and cellular level necessitates developing novel individualized treatment modalities via genomics-driven approaches.

Methods

This study leverages numerous pharmacogenomic and tissue databases to examine drug repositioning for glioblastoma. RNAseq of glioblastoma tumor samples from The Cancer Genome Atlas (TCGA, n=117) were compared to “normal” frontal lobe samples from Genotype-Tissue Expression Portal (GTEX, n=120) to find differentially expressed genes (DEGs). Using compound-gene expression data and drug activity data from the Library of Integrated Network-Based Cellular Signatures (LINCS, n=66,512 compounds) CCLE (71 glioma cell lines), and Chemical European Molecular Biology Laboratory (ChEMBL) platforms, we employed a summarized reversal gene expression metric (sRGES) to “reverse” the resultant disease signature for GBM and its subtypes. A multi-parametric strategy was employed to stratify compounds capable of blood brain barrier penetrance with a favorable pharmacokinetic profile (CNS-MPO).

Results

Significant correlations were identified between sRGES and drug efficacy in GBM cell lines in both ChEMBL(r=0.37,p&lt;.001) and Cancer Therapeutic Response Portal (CTRP) databases (r=0.35, p&lt;0.001). Our multiparametric algorithm identified two classes of drugs with highest sRGES and CNS-MPO: HDAC inhibitors (vorinostat and entinostat) and topoisomerase inhibitors suitable for drug repurposing.

Conclusions

Our studies suggest that reversal of glioblastoma disease signature correlates with drug potency for various GBM subtypes. This multiparametric approach may set the foundation for an early-phase personalized -omics clinical trial for glioblastoma by effectively identifying drugs that are capable of reversing the disease signature and have favorable pharmacokinetic and safety profiles.

A survey of RNA editing at single-cell resolution links interneurons to schizophrenia and autism

Conversion of adenosine to inosine in RNA by ADAR enzymes, termed "RNA editing," is essential for healthy brain development. Editing is dysregulated in neuropsychiatric diseases, but has not yet been investigated at scale at the level of individual neurons. We quantified RNA editing sites in nuclear transcriptomes of 3055 neurons from six cortical regions of a neurotypical female donor, and found 41,930 sites present in at least ten nuclei. Most sites were located within Alu repeats in introns or 3' UTRs, and approximately 80% were cataloged in public RNA editing databases. We identified 9285 putative novel editing sites, 29% of which were also detectable in unrelated donors. Intersection with results from bulk RNA-seq studies provided cell-type and spatial context for 1730 sites that are differentially edited in schizophrenic brain donors, and 910 such sites in autistic donors. Autism-related genes were also enriched with editing sites predicted to modify RNA structure. Inhibitory neurons showed higher overall transcriptome editing than excitatory neurons, and the highest editing rates were observed in the frontal cortex. We used generalized linear models to identify differentially edited sites and genes between cell types. Twenty nine genes were preferentially edited in excitatory neurons, and 43 genes were edited more heavily in inhibitory neurons, including RBFOX1, its target genes, and genes in the autism-associated Prader-Willi locus (15q11). The abundance of SNORD115/116 genes from locus 15q11 was positively associated with editing activity across the transcriptome. We contend that insufficient editing of autism-related genes in inhibitory neurons may contribute to the specific perturbation of those cells in autism.

recount3: summaries and queries for large-scale RNA-seq expression and splicing

We present recount3, a resource consisting of over 750,000 publicly available human and mouse RNA sequencing (RNA-seq) samples uniformly processed by our new Monorail analysis pipeline. To facilitate access to the data, we provide the recount3 and snapcount R/Bioconductor packages as well as complementary web resources. Using these tools, data can be downloaded as study-level summaries or queried for specific exon-exon junctions, genes, samples, or other features. Monorail can be used to process local and/or private data, allowing results to be directly compared to any study in recount3. Taken together, our tools help biologists maximize the utility of publicly available RNA-seq data, especially to improve their understanding of newly collected data. recount3 is available from http://rna.recount.bio .

In Silico Methods for the Identification of Diagnostic and Favorable Prognostic Markers in Acute Myeloid Leukemia

Acute myeloid leukemia (AML), the most common type of acute leukemia in adults, is mainly asymptomatic at early stages and progresses/recurs rapidly and frequently. These attributes necessitate the identification of biomarkers for timely diagnosis and accurate prognosis. In this study, differential gene expression analysis was performed on large-scale transcriptomics data of AML patients versus corresponding normal tissue. Weighted gene co-expression network analysis was conducted to construct networks of co-expressed genes, and detect gene modules. Finally, hub genes were identified from selected modules by applying network-based methods. This robust and integrative bioinformatics approach revealed a set of twenty-four genes, mainly related to cell cycle and immune response, the diagnostic significance of which was subsequently compared against two independent gene expression datasets. Furthermore, based on a recent notion suggesting that molecular characteristics of a few, unusual patients with exceptionally favorable survival can provide insights for improving the outcome of individuals with more typical disease trajectories, we defined groups of long-term survivors in AML patient cohorts and compared their transcriptomes versus the general population to infer favorable prognostic signatures. These findings could have potential applications in the clinical setting, in particular, in diagnosis and prognosis of AML.

clusterProfiler 4.0: A universal enrichment tool for interpreting omics data

Functional enrichment analysis is pivotal for interpreting high-throughput omics data in life science. It is crucial for this type of tool to use the latest annotation databases for as many organisms as possible. To meet these requirements, we present here an updated version of our popular Bioconductor package, clusterProfiler 4.0. This package has been enhanced considerably compared with its original version published 9 years ago. The new version provides a universal interface for functional enrichment analysis in thousands of organisms based on internally supported ontologies and pathways as well as annotation data provided by users or derived from online databases. It also extends the dplyr and ggplot2 packages to offer tidy interfaces for data operation and visualization. Other new features include gene set enrichment analysis and comparison of enrichment results from multiple gene lists. We anticipate that clusterProfiler 4.0 will be applied to a wide range of scenarios across diverse organisms.

•

clusterProfiler supports exploring functional characteristics of both coding and non-coding genomics data for thousands of species with up-to-date gene annotation

•

It provides a universal interface for gene functional annotation from a variety of sources and thus can be applied in diverse scenarios

•

It provides a tidy interface to access, manipulate, and visualize enrichment results to help users achieve efficient data interpretation

•

Datasets obtained from multiple treatments and time points can be analyzed and compared in a single run, easily revealing functional consensus and differences among distinct conditions

Personalized single-cell networks: a framework to predict the response of any gene to any drug for any patient

BACKGROUND

The last decade has seen a major increase in the availability of genomic data. This includes expert-curated databases that describe the biological activity of genes, as well as high-throughput assays that measure gene expression in bulk tissue and single cells. Integrating these heterogeneous data sources can generate new hypotheses about biological systems. Our primary objective is to combine population-level drug-response data with patient-level single-cell expression data to predict how any gene will respond to any drug for any patient.

METHODS

We take 2 approaches to benchmarking a "dual-channel" random walk with restart (RWR) for data integration. First, we evaluate how well RWR can predict known gene functions from single-cell gene co-expression networks. Second, we evaluate how well RWR can predict known drug responses from individual cell networks. We then present two exploratory applications. In the first application, we combine the Gene Ontology database with glioblastoma single cells from 5 individual patients to identify genes whose functions differ between cancers. In the second application, we combine the LINCS drug-response database with the same glioblastoma data to identify genes that may exhibit patient-specific drug responses.

CONCLUSIONS

Our manuscript introduces two innovations to the integration of heterogeneous biological data. First, we use a "dual-channel" method to predict up-regulation and down-regulation separately. Second, we use individualized single-cell gene co-expression networks to make personalized predictions. These innovations let us predict gene function and drug response for individual patients. Taken together, our work shows promise that single-cell co-expression data could be combined in heterogeneous information networks to facilitate precision medicine.

FGFR1 overexpression in non-small cell lung cancer is mediated by genetic and epigenetic mechanisms and is a determinant of FGFR1 inhibitor response

Amplification of fibroblast growth factor receptor 1 (FGFR1) in non-small cell lung cancer (NSCLC) has been considered as an actionable drug target. However, pan-FGFR tyrosine kinase inhibitors did not demonstrate convincing clinical efficacy in FGFR1-amplified NSCLC patients. This study aimed to characterise the molecular context of FGFR1 expression and to define biomarkers predictive of FGFR1 inhibitor response. In this study, 635 NSCLC samples were characterised for FGFR1 protein expression by immunohistochemistry and copy number gain (CNG) by in situ hybridisation (n = 298) or DNA microarray (n = 189). FGFR1 gene expression (n = 369) and immune cell profiles (n = 309) were also examined. Furthermore, gene expression, methylation and microRNA data from The Cancer Genome Atlas (TCGA) were compared. A panel of FGFR1-amplified NSCLC patient-derived xenograft (PDX) models were tested for response to the selective FGFR1 antagonist M6123. A minority of patients demonstrated FGFR1 CNG (10.5%) or increased FGFR1 mRNA (8.7%) and protein expression (4.4%). FGFR1 CNG correlated weakly with FGFR1 gene and protein expression. Tumours overexpressing FGFR1 protein were typically devoid of driver alterations (e.g. EGFR, KRAS) and showed reduced infiltration of T-lymphocytes and lower PD-L1 expression. Promoter methylation and microRNA were identified as regulators of FGFR1 expression in NSCLC and other cancers. Finally, NSCLC PDX models demonstrating FGFR1 amplification and FGFR1 protein overexpression were sensitive to M6123. The unique molecular and immune features of tumours with high FGFR1 expression provide a rationale to stratify patients in future clinical trials of FGFR1 pathway-targeting agents.

Opportunities and Challenges in Democratizing Immunology Datasets

The field of immunology is rapidly progressing toward a systems-level understanding of immunity to tackle complex infectious diseases, autoimmune conditions, cancer, and beyond. In the last couple of decades, advancements in data acquisition techniques have presented opportunities to explore untapped areas of immunological research. Broad initiatives are launched to disseminate the datasets siloed in the global, federated, or private repositories, facilitating interoperability across various research domains. Concurrently, the application of computational methods, such as network analysis, meta-analysis, and machine learning have propelled the field forward by providing insight into salient features that influence the immunological response, which was otherwise left unexplored. Here, we review the opportunities and challenges in democratizing datasets, repositories, and community-wide knowledge sharing tools. We present use cases for repurposing open-access immunology datasets with advanced machine learning applications and more.

Comparative epigenome-wide analysis highlights placenta-specific differentially methylated regions

Aim: The placenta undergoes DNA methylation (DNAm) programming that is unique compared with all other fetal tissues. We aim to decipher some of the physiologic roles of the placenta by comparing its DNAm profile with that of another fetal tissue. Materials & methods: We performed a comparative analysis of genome-wide DNAm of 444 placentas paired with cord blood samples collected at birth. Gene ontology term analyses were conducted on the resulting differentially methylated regions. Results: Genomic regions upstream of transcription start sites showing lower DNAm in the placenta were enriched with terms related to miRNA functions and genes encoding G-protein-coupled receptors. Conclusion: These results highlight genomic regions that are differentially methylated in the placenta in contrast to fetal blood.

Style transfer with variational autoencoders is a promising approach to RNA-Seq data harmonization and analysis

Motivation

The transcriptomic data are being frequently used in the research of biomarker genes of different diseases and biological states. The most common tasks there are the data harmonization and treatment outcome prediction. Both of them can be addressed via the style transfer approach. Either technical factors or any biological details about the samples which we would like to control (gender, biological state, treatment, etc.) can be used as style components.

Results

The proposed style transfer solution is based on Conditional Variational Autoencoders, Y-Autoencoders and adversarial feature decomposition. To quantitatively measure the quality of the style transfer, neural network classifiers which predict the style and semantics after training on real expression were used. Comparison with several existing style-transfer based approaches shows that proposed model has the highest style prediction accuracy on all considered datasets while having comparable or the best semantics prediction accuracy.

Availability and implementation

https://github.com/NRshka/stvae-source.

Supplementary information

Supplementary data are available at Bioinformatics online.

Maximizing the reusability of gene expression data by predicting missing metadata

Reusability is part of the FAIR data principle, which aims to make data Findable, Accessible, Interoperable, and Reusable. One of the current efforts to increase the reusability of public genomics data has been to focus on the inclusion of quality metadata associated with the data. When necessary metadata are missing, most researchers will consider the data useless. In this study, we developed a framework to predict the missing metadata of gene expression datasets to maximize their reusability. We found that when using predicted data to conduct other analyses, it is not optimal to use all the predicted data. Instead, one should only use the subset of data, which can be predicted accurately. We proposed a new metric called Proportion of Cases Accurately Predicted (PCAP), which is optimized in our specifically-designed machine learning pipeline. The new approach performed better than pipelines using commonly used metrics such as F1-score in terms of maximizing the reusability of data with missing values. We also found that different variables might need to be predicted using different machine learning methods and/or different data processing protocols. Using differential gene expression analysis as an example, we showed that when missing variables are accurately predicted, the corresponding gene expression data can be reliably used in downstream analyses.

Low Entropy Sub-Networks Prevent the Integration of Metabolomic and Transcriptomic Data

The constantly and rapidly increasing amount of the biological data gained from many different high-throughput experiments opens up new possibilities for data- and model-driven inference. Yet, alongside, emerges a problem of risks related to data integration techniques. The latter are not so widely taken account of. Especially, the approaches based on the flux balance analysis (FBA) are sensitive to the structure of a metabolic network for which the low-entropy clusters can prevent the inference from the activity of the metabolic reactions. In the following article, we set forth problems that may arise during the integration of metabolomic data with gene expression datasets. We analyze common pitfalls, provide their possible solutions, and exemplify them by a case study of the renal cell carcinoma (RCC). Using the proposed approach we provide a metabolic description of the known morphological RCC subtypes and suggest a possible existence of the poor-prognosis cluster of patients, which are commonly characterized by the low activity of the drug transporting enzymes crucial in the chemotherapy. This discovery suits and extends the already known poor-prognosis characteristics of RCC. Finally, the goal of this work is also to point out the problem that arises from the integration of high-throughput data with the inherently nonuniform, manually curated low-throughput data. In such cases, the over-represented information may potentially overshadow the non-trivial discoveries.

PyGNA: a unified framework for geneset network analysis

BACKGROUND

Gene and protein interaction experiments provide unique opportunities to study the molecular wiring of a cell. Integrating high-throughput functional genomics data with this information can help identifying networks associated with complex diseases and phenotypes.

RESULTS

Here we introduce an integrated statistical framework to test network properties of single and multiple genesets under different interaction models. We implemented this framework as an open-source software, called Python Geneset Network Analysis (PyGNA). Our software is designed for easy integration into existing analysis pipelines and to generate high quality figures and reports. We also developed PyGNA to take advantage of multi-core systems to generate calibrated null distributions on large datasets. We then present the results of extensive benchmarking of the tests implemented in PyGNA and a use case inspired by RNA sequencing data analysis, showing how PyGNA can be easily integrated to study biological networks. PyGNA is available at http://github.com/stracquadaniolab/pygna and can be easily installed using the PyPi or Anaconda package managers, and Docker.

CONCLUSIONS

We present a tool for network-aware geneset analysis. PyGNA can either be readily used and easily integrated into existing high-performance data analysis pipelines or as a Python package to implement new tests and analyses. With the increasing availability of population-scale omic data, PyGNA provides a viable approach for large scale geneset network analysis.

Epigenomic and Transcriptomic Dynamics During Human Heart Organogenesis

RATIONALE

There is growing evidence that common variants and rare sequence alterations in regulatory sequences can result in birth defects or predisposition to disease. Congenital heart defects are the most common birth defect and have a clear genetic component, yet only a third of cases can be attributed to structural variation in the genome or a mutation in a gene. The remaining unknown cases could be caused by alterations in regulatory sequences.

OBJECTIVE

Identify regulatory sequences and gene expression networks that are active during organogenesis of the human heart. Determine whether these sites and networks are enriched for disease-relevant genes and associated genetic variation.

METHODS AND RESULTS

We characterized ChromHMM (chromatin state) and gene expression dynamics during human heart organogenesis. We profiled 7 histone modifications in embryonic hearts from each of 9 distinct Carnegie stages (13-14, 16-21, and 23), annotated chromatin states, and compared these maps to over 100 human tissues and cell types. We also generated RNA-sequencing data, performed differential expression, and constructed weighted gene coexpression networks. We identified 177 412 heart enhancers; 12 395 had not been previously annotated as strong enhancers. We identified 92% of all functionally validated heart-positive enhancers (n=281; 7.5× enrichment; P<2.2×10-16). Integration of these data demonstrated novel heart enhancers are enriched near genes expressed more strongly in cardiac tissue and are enriched for variants associated with ECG measures and atrial fibrillation. Our gene expression network analysis identified gene modules strongly enriched for heart-related functions, regulatory control by heart-specific enhancers, and putative disease genes.

CONCLUSIONS

Well-connected hub genes with heart-specific expression targeted by embryonic heart-specific enhancers are likely disease candidates. Our functional annotations will allow for better interpretation of whole genome sequencing data in the large number of patients affected by congenital heart defects.

A pan-cancer assessment of alterations of the kinase domain of ULK1, an upstream regulator of autophagy

Autophagy is a key clearance process to recycle damaged cellular components. One important upstream regulator of autophagy is ULK1 kinase. Several three-dimensional structures of the ULK1 catalytic domain are available, but a comprehensive study, including molecular dynamics, is missing. Also, an exhaustive description of ULK1 alterations found in cancer samples is presently lacking. We here applied a framework which links -omics data to structural protein ensembles to study ULK1 alterations from genomics data available for more than 30 cancer types. We predicted the effects of mutations on ULK1 function and structural stability, accounting for protein dynamics, and the different layers of changes that a mutation can induce in a protein at the functional and structural level. ULK1 is down-regulated in gynecological tumors. In other cancer types, ULK2 could compensate for ULK1 downregulation and, in the majority of the cases, no marked changes in expression have been found. 36 missense mutations of ULK1, not limited to the catalytic domain, are co-occurring with mutations in a large number of ULK1 interactors or substrates, suggesting a pronounced effect of the upstream steps of autophagy in many cancer types. Moreover, our results pinpoint that more than 50% of the mutations in the kinase domain of ULK1, here investigated, are predicted to affect protein stability. Three mutations (S184F, D102N, and A28V) are predicted with only impact on kinase activity, either modifying the functional dynamics or the capability to exert effects from distal sites to the functional and catalytic regions. The framework here applied could be extended to other protein targets to aid the classification of missense mutations from cancer genomics studies, as well as to prioritize variants for experimental validation, or to select the appropriate biological readouts for experiments.

SliceIt: A genome-wide resource and visualization tool to design CRISPR/Cas9 screens for editing protein-RNA interaction sites in the human genome

Several protein-RNA cross linking protocols have been established in recent years to delineate the molecular interaction of an RNA Binding Protein (RBP) and its target RNAs. However, functional dissection of the role of the RBP binding sites in modulating the post-transcriptional fate of the target RNA remains challenging. CRISPR/Cas9 genome editing system is being commonly employed to perturb both coding and noncoding regions in the genome. With the advancements in genome-scale CRISPR/Cas9 screens, it is now possible to not only perturb specific binding sites but also probe the global impact of protein-RNA interaction sites across cell types. Here, we present SliceIt (http://sliceit.soic.iupui.edu/), a database of in silico sgRNA (single guide RNA) library to facilitate conducting such high throughput screens. SliceIt comprises of ~4.8 million unique sgRNAs with an estimated range of 2-8 sgRNAs designed per RBP binding site, for eCLIP experiments of >100 RBPs in HepG2 and K562 cell lines from the ENCODE project. SliceIt provides a user friendly environment, developed using advanced search engine framework, Elasticsearch. It is available in both table and genome browser views facilitating the easy navigation of RBP binding sites, designed sgRNAs, exon expression levels across 53 human tissues along with prevalence of SNPs and GWAS hits on binding sites. Exon expression profiles enable examination of locus specific changes proximal to the binding sites. Users can also upload custom tracks of various file formats directly onto genome browser, to navigate additional genomic features in the genome and compare with other types of omics profiles. All the binding site-centric information is dynamically accessible via "search by gene", "search by coordinates" and "search by RBP" options and readily available to download. Validation of the sgRNA library in SliceIt was performed by selecting RBP binding sites in Lipt1 gene and designing sgRNAs. Effect of CRISPR/Cas9 perturbations on the selected binding sites in HepG2 cell line, was confirmed based on altered proximal exon expression levels using qPCR, further supporting the utility of the resource to design experiments for perturbing protein-RNA interaction networks. Thus, SliceIt provides a one-stop repertoire of guide RNA library to perturb RBP binding sites, along with several layers of functional information to design both low and high throughput CRISPR/Cas9 screens, for studying the phenotypes and diseases associated with RBP binding sites.

Investigating the effect of dependence between conditions with Bayesian Linear Mixed Models for motif activity analysis

MOTIVATION

Cellular identity and behavior is controlled by complex gene regulatory networks. Transcription factors (TFs) bind to specific DNA sequences to regulate the transcription of their target genes. On the basis of these TF motifs in cis-regulatory elements we can model the influence of TFs on gene expression. In such models of TF motif activity the data is usually modeled assuming a linear relationship between the motif activity and the gene expression level. A commonly used method to model motif influence is based on Ridge Regression. One important assumption of linear regression is the independence between samples. However, if samples are generated from the same cell line, tissue, or other biological source, this assumption may be invalid. This same assumption of independence is also applied to different yet similar experimental conditions, which may also be inappropriate. In theory, the independence assumption between samples could lead to loss in signal detection. Here we investigate whether a Bayesian model that allows for correlations results in more accurate inference of motif activities.

RESULTS

We extend the Ridge Regression to a Bayesian Linear Mixed Model, which allows us to model dependence between different samples. In a simulation study, we investigate the differences between the two model assumptions. We show that our Bayesian Linear Mixed Model implementation outperforms Ridge Regression in a simulation scenario where the noise, which is the signal that can not be explained by TF motifs, is uncorrelated. However, we demonstrate that there is no such gain in performance if the noise has a similar covariance structure over samples as the signal that can be explained by motifs. We give a mathematical explanation to why this is the case. Using four representative real datasets we show that at most ∼​40% of the signal is explained by motifs using the linear model. With these data there is no advantage to using the Bayesian Linear Mixed Model, due to the similarity of the covariance structure.

AVAILABILITY & IMPLEMENTATION

The project implementation is available at https://github.com/Sim19/SimGEXPwMotifs.

Consistent RNA sequencing contamination in GTEx and other data sets

A challenge of next generation sequencing is read contamination. We use Genotype-Tissue Expression (GTEx) datasets and technical metadata along with RNA-seq datasets from other studies to understand factors that contribute to contamination. Here we report, of 48 analyzed tissues in GTEx, 26 have variant co-expression clusters of four highly expressed and pancreas-enriched genes (PRSS1, PNLIP, CLPS, and/or CELA3A). Fourteen additional highly expressed genes from other tissues also indicate contamination. Sample contamination is strongly associated with a sample being sequenced on the same day as a tissue that natively expresses those genes. Discrepant SNPs across four contaminating genes validate the contamination. Low-level contamination affects ~40% of samples and leads to numerous eQTL assignments in inappropriate tissues among these 18 genes. This type of contamination occurs widely, impacting bulk and single cell (scRNA-seq) data set analysis. In conclusion, highly expressed, tissue-enriched genes basally contaminate GTEx and other datasets impacting analyses.

Putatively cancer-specific exon-exon junctions are shared across patients and present in developmental and other non-cancer cells

This study probes the distribution of putatively cancer-specific junctions across a broad set of publicly available non-cancer human RNA sequencing (RNA-seq) datasets. We compared cancer and non-cancer RNA-seq data from The Cancer Genome Atlas (TCGA), the Genotype-Tissue Expression (GTEx) Project and the Sequence Read Archive. We found that (i) averaging across cancer types, 80.6% of exon-exon junctions thought to be cancer-specific based on comparison with tissue-matched samples (σ = 13.0%) are in fact present in other adult non-cancer tissues throughout the body; (ii) 30.8% of junctions not present in any GTEx or TCGA normal tissues are shared by multiple samples within at least one cancer type cohort, and 87.4% of these distinguish between different cancer types; and (iii) many of these junctions not found in GTEx or TCGA normal tissues (15.4% on average, σ = 2.4%) are also found in embryological and other developmentally associated cells. These findings refine the meaning of RNA splicing event novelty, particularly with respect to the human neoepitope repertoire. Ultimately, cancer-specific exon-exon junctions may have a substantial causal relationship with the biology of disease.

Recounting the FANTOM CAGE-Associated Transcriptome

Long noncoding RNAs (lncRNAs) have emerged as key coordinators of biological and cellular processes. Characterizing lncRNA expression across cells and tissues is key to understanding their role in determining phenotypes, including human diseases. We present here FC-R2, a comprehensive expression atlas across a broadly defined human transcriptome, inclusive of over 109,000 coding and noncoding genes, as described in the FANTOM CAGE-Associated Transcriptome (FANTOM-CAT) study. This atlas greatly extends the gene annotation used in the original recount2 resource. We demonstrate the utility of the FC-R2 atlas by reproducing key findings from published large studies and by generating new results across normal and diseased human samples. In particular, we (a) identify tissue-specific transcription profiles for distinct classes of coding and noncoding genes, (b) perform differential expression analysis across thirteen cancer types, identifying novel noncoding genes potentially involved in tumor pathogenesis and progression, and (c) confirm the prognostic value for several enhancer lncRNAs expression in cancer. Our resource is instrumental for the systematic molecular characterization of lncRNA by the FANTOM6 Consortium. In conclusion, comprised of over 70,000 samples, the FC-R2 atlas will empower other researchers to investigate functions and biological roles of both known coding genes and novel lncRNAs.

Using R and Bioconductor in Clinical Genomics and Transcriptomics

Bioinformatics pipelines are essential in the analysis of genomic and transcriptomic data generated by next-generation sequencing (NGS). Recent guidelines emphasize the need for rigorous validation and assessment of robustness, reproducibility, and quality of NGS analytic pipelines intended for clinical use. Software tools written in the R statistical language and, in particular, the set of tools available in the Bioconductor repository are widely used in research bioinformatics; and these frameworks offer several advantages for use in clinical bioinformatics, including the breath of available tools, modular nature of software packages, ease of installation, enforcement of interoperability, version control, and short learning curve. This review provides an introduction to R and Bioconductor software, its advantages and limitations for clinical bioinformatics, and illustrative examples of tools that can be used in various steps of NGS analysis.

Improving the diagnostic yield of exome- sequencing by predicting gene-phenotype associations using large-scale gene expression analysis

The diagnostic yield of exome and genome sequencing remains low (8-70%), due to incomplete knowledge on the genes that cause disease. To improve this, we use RNA-seq data from 31,499 samples to predict which genes cause specific disease phenotypes, and develop GeneNetwork Assisted Diagnostic Optimization (GADO). We show that this unbiased method, which does not rely upon specific knowledge on individual genes, is effective in both identifying previously unknown disease gene associations, and flagging genes that have previously been incorrectly implicated in disease. GADO can be run on www.genenetwork.nl by supplying HPO-terms and a list of genes that contain candidate variants. Finally, applying GADO to a cohort of 61 patients for whom exome-sequencing analysis had not resulted in a genetic diagnosis, yields likely causative genes for ten cases.

New functionalities in the TCGAbiolinks package for the study and integration of cancer data from GDC and GTEx

The advent of Next-Generation Sequencing (NGS) technologies has opened new perspectives in deciphering the genetic mechanisms underlying complex diseases. Nowadays, the amount of genomic data is massive and substantial efforts and new tools are required to unveil the information hidden in the data. The Genomic Data Commons (GDC) Data Portal is a platform that contains different genomic studies including the ones from The Cancer Genome Atlas (TCGA) and the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) initiatives, accounting for more than 40 tumor types originating from nearly 30000 patients. Such platforms, although very attractive, must make sure the stored data are easily accessible and adequately harmonized. Moreover, they have the primary focus on the data storage in a unique place, and they do not provide a comprehensive toolkit for analyses and interpretation of the data. To fulfill this urgent need, comprehensive but easily accessible computational methods for integrative analyses of genomic data that do not renounce a robust statistical and theoretical framework are required. In this context, the R/Bioconductor package TCGAbiolinks was developed, offering a variety of bioinformatics functionalities. Here we introduce new features and enhancements of TCGAbiolinks in terms of i) more accurate and flexible pipelines for differential expression analyses, ii) different methods for tumor purity estimation and filtering, iii) integration of normal samples from other platforms iv) support for other genomics datasets, exemplified here by the TARGET data. Evidence has shown that accounting for tumor purity is essential in the study of tumorigenesis, as these factors promote confounding behavior regarding differential expression analysis. With this in mind, we implemented these filtering procedures in TCGAbiolinks. Moreover, a limitation of some of the TCGA datasets is the unavailability or paucity of corresponding normal samples. We thus integrated into TCGAbiolinks the possibility to use normal samples from the Genotype-Tissue Expression (GTEx) project, which is another large-scale repository cataloging gene expression from healthy individuals. The new functionalities are available in the TCGAbiolinks version 2.8 and higher released in Bioconductor version 3.7.

The advent of Next-Generation Sequencing (NGS) technologies has been generating a massive amount of data which require continuous efforts in developing and maintain computational tool for data analyses. The Genomic Data Commons (GDC) Data Portal is a platform that contains different cancer genomic studies. Such platforms have often the primary focus on the data storage and they do not provide a comprehensive toolkit for analyses. To fulfil this urgent need, comprehensive but accessible computational protocols that do not renounce a robust statistical framework are thus required. In this context, we here present the new functions of the R/Bioconductor package TCGAbiolinks to improve the discovery of differentially expressed genes in cancer and tumor (sub)types, include the estimate of tumor purity and tumor infiltrations, use normal samples from other platforms and support more broadly other genomics datasets.