Publications

Background Visium is a widely-used spatially-resolved transcriptomics assay available from 10x Genomics. Standard Visium capture areas (6.5mm by 6.5mm) limit the survey of larger tissue structures, but combining overlapping images and associated gene expression data allow for more complex study designs. Current software can handle nested or partial image overlaps, but is designed for merging up to two capture areas, and cannot account for some technical scenarios related to capture area alignment. Results We generated Visium data from a postmortem human tissue sample such that two capture areas were partially overlapping and a third one was adjacent. We developed the R/Bioconductor package visiumStitched, which facilitates stitching the images together with Fiji (ImageJ), and constructing SpatialExperiment R objects with the stitched images and gene expression data. visiumStitched constructs an artificial hexagonal array grid which allows seamless downstream analyses such as spatially-aware clustering without discarding data from overlapping spots. Data stitched with visiumStitched can then be interactively visualized with spatialLIBD. Conclusions visiumStitched provides a simple, but flexible framework to handle various multi-capture area study design scenarios. Specifically, it resolves a data processing step without disrupting analysis workflows and without discarding data from overlapping spots. visiumStitched relies on affine transformations by Fiji, which have limitations and are less accurate when aligning against an atlas or other situations. visiumStitched provides an easy-to-use solution which expands possibilities for designing multi-capture area study designs.

Maternal smoking during pregnancy (MSDP) is associated with significant cognitive and behavioral effects on offspring. While neurodevelopmental outcomes have been studied for prenatal exposure to nicotine, the main psychoactive component of cigarette smoke, its contribution to MSDP effects has never been explored. Comparing the effects of these substances on molecular signaling in the prenatal and adult brain may provide insights into nicotinic and broader tobacco consequences that are developmental-stage specific or age-independent. Pregnant mice were administered nicotine or exposed to chronic cigarette smoke, and RNA-sequencing was performed on frontal cortices of postnatal day 0 pups born to these mice, as well as on frontal cortices and blood of the adult dams. We identified 1,010 and 4,165 differentially expressed genes (DEGs) in nicotine and smoking-exposed pup brains, respectively (FDR<0.05, Ns = 19 nicotine-exposed vs 23 vehicle-exposed; 46 smoking-exposed vs 49 controls). Prenatal nicotine exposure (PNE) alone was related to dopaminergic synapses and long-term synaptic depression, whereas MSDP was associated with the SNARE complex and vesicle transport. Both substances affected SMN-Sm protein complexes and postsynaptic endosomes. Analyses at the transcript, exon, and exon-exon junction levels supported gene level results and revealed additional smoking-affected processes. No DEGs at FDR<0.05 were found in adult mouse brain for any substance (12 nicotine-administered vs 11 vehicle-administered; 12 smoking-exposed vs 12 controls), nor in adult blood (12 smoking-exposed vs 12 controls), and only 3% and 6.41% of the DEGs in smoking-exposed pup brain replicated in smoking-exposed blood and human prenatal brain, respectively. Together, these results demonstrate variable but overlapping molecular effects of PNE and MSDP on the developing brain, and attenuated effects of both smoking and nicotine on adult versus fetal brain.

Neurodegenerative disorders have overlapping symptoms and have high comorbidity rates, but this is not reflected in overlaps of risk genes. We have investigated whether ligand-receptor interactions (LRIs) are a mechanism by which distinct genes associated with disease risk can impact overlapping outcomes. We found that LRIs are likely disrupted in neurological disease and that the ligand-receptor networks associated with neurological diseases have substantial overlaps. Specifically, 96.8% of LRIs associated with disease risk are interconnected in a single LR network. These ligands and receptors are enriched for roles in inflammatory pathways and highlight the role of glia in cross-disease risk. Disruption to this LR network due to disease-associated processes (e.g. differential transcript use, protein misfolding) is likely to contribute to disease progression and risk of comorbidity. Our findings have implications for drug development, as they highlight the potential benefits and risks of pursuing cross-disease drug targets.

The molecular pathology of stress-related disorders remains elusive. Our brain multiregion, multiomic study of posttraumatic stress disorder (PTSD) and major depressive disorder (MDD) included the central nucleus of the amygdala, hippocampal dentate gyrus, and medial prefrontal cortex (mPFC). Genes and exons within the mPFC carried most disease signals replicated across two independent cohorts. Pathways pointed to immune function, neuronal and synaptic regulation, and stress hormones. Multiomic factor and gene network analyses provided the underlying genomic structure. Single nucleus RNA sequencing in dorsolateral PFC revealed dysregulated (stress-related) signals in neuronal and non-neuronal cell types. Analyses of brain-blood intersections in >50,000 UK Biobank participants were conducted along with fine-mapping of the results of PTSD and MDD genome-wide association studies to distinguish risk from disease processes. Our data suggest shared and distinct molecular pathology in both disorders and propose potential therapeutic targets and biomarkers.

The molecular organization of the human neocortex historically has been studied in the context of its histological layers. However, emerging spatial transcriptomic technologies have enabled unbiased identification of transcriptionally defined spatial domains that move beyond classic cytoarchitecture. We used the Visium spatial gene expression platform to generate a data-driven molecular neuroanatomical atlas across the anterior-posterior axis of the human dorsolateral prefrontal cortex. Integration with paired single-nucleus RNA-sequencing data revealed distinct cell type compositions and cell-cell interactions across spatial domains. Using PsychENCODE and publicly available data, we mapped the enrichment of cell types and genes associated with neuropsychiatric disorders to discrete spatial domains.

Ancestral differences in genomic variation affect the regulation of gene expression; however, most gene expression studies have been limited to European ancestry samples or adjusted to identify ancestry-independent associations. Here, we instead examined the impact of genetic ancestry on gene expression and DNA methylation in the postmortem brain tissue of admixed Black American neurotypical individuals to identify ancestry-dependent and ancestry-independent contributions. Ancestry-associated differentially expressed genes (DEGs), transcripts and gene networks, while notably not implicating neurons, are enriched for genes related to the immune response and vascular tissue and explain up to 26% of heritability for ischemic stroke, 27% of heritability for Parkinson disease and 30% of heritability for Alzheimer’s disease. Ancestry-associated DEGs also show general enrichment for the heritability of diverse immune-related traits but depletion for psychiatric-related traits. We also compared Black and non-Hispanic white Americans, confirming most ancestry-associated DEGs. Our results delineate the extent to which genetic ancestry affects differences in gene expression in the human brain and the implications for brain illness risk.

Schizophrenia is a complex neuropsychiatric disorder with sexually dimorphic features, including differential symptomatology, drug responsiveness, and male incidence rate. Prior large-scale transcriptome analyses for sex differences in schizophrenia have focused on the prefrontal cortex. Analyzing BrainSeq Consortium data (caudate nucleus: n = 399, dorsolateral prefrontal cortex: n = 377, and hippocampus: n = 394), we identified 831 unique genes that exhibit sex differences across brain regions, enriched for immune-related pathways. We observed X-chromosome dosage reduction in the hippocampus of male individuals with schizophrenia. Our sex interaction model revealed 148 junctions dysregulated in a sex-specific manner in schizophrenia. Sex-specific schizophrenia analysis identified dozens of differentially expressed genes, notably enriched in immune-related pathways. Finally, our sex-interacting expression quantitative trait loci analysis revealed 704 unique genes, nine associated with schizophrenia risk. These findings emphasize the importance of sex-informed analysis of sexually dimorphic traits, inform personalized therapeutic strategies in schizophrenia, and highlight the need for increased female samples for schizophrenia analyses.

Pathophysiology of many neuropsychiatric disorders, including schizophrenia (SCZD), is linked to habenula (Hb) function. While pharmacotherapies and deep brain stimulation targeting the Hb are emerging as promising therapeutic treatments, little is known about the cell type-specific transcriptomic organization of the human Hb or how it is altered in SCZD. Here we define the molecular neuroanatomy of the human Hb and identify transcriptomic changes in individuals with SCZD compared to neurotypical controls. Utilizing Hb-enriched postmortem human brain tissue, we performed single nucleus RNA-sequencing (snRNA-seq; n=7 neurotypical donors) and identified 17 molecularly defined Hb cell types across 16,437 nuclei, including 3 medial and 7 lateral Hb populations, several of which were conserved between rodents and humans. Single molecule fluorescent in situ hybridization (smFISH; n=3 neurotypical donors) validated snRNA-seq Hb cell types and mapped their spatial locations. Bulk RNA-sequencing and cell type deconvolution in Hb-enriched tissue from 35 individuals with SCZD and 33 neurotypical controls yielded 45 SCZD-associated differentially expressed genes (DEGs, FDR < 0.05), with 32 (71%) unique to Hb-enriched tissue. eQTL analysis identified 717 independent SNP-gene pairs (FDR < 0.05), where either the SNP is a SCZD risk variant (16 pairs) or the gene is a SCZD DEG (7 pairs). eQTL and SCZD risk colocalization analysis identified 16 colocalized genes. These results identify topographically organized cell types with distinct molecular signatures in the human Hb and demonstrate unique genetic changes associated with SCZD, thereby providing novel molecular insights into the role of Hb in neuropsychiatric disorders.

Background Cellular deconvolution of bulk RNA-sequencing (RNA-seq) data using single cell or nuclei RNA-seq (sc/snRNA-seq) reference data is an important strategy for estimating cell type composition in heterogeneous tissues, such as human brain. Computational methods for deconvolution have been developed and benchmarked against simulated data, pseudobulked sc/snRNA-seq data, or immunohistochemistry reference data. A major limitation in developing improved deconvolution algorithms has been the lack of integrated datasets with orthogonal measurements of gene expression and estimates of cell type proportions on the same tissue sample. Deconvolution algorithm performance has not yet been evaluated across different RNA extraction methods (cytosolic, nuclear, or whole cell RNA), different library preparation types (mRNA enrichment vs. ribosomal RNA depletion), or with matched single cell reference datasets. Results A rich multi-assay dataset was generated in postmortem human dorsolateral prefrontal cortex (DLPFC) from 22 tissue blocks. Assays included spatially-resolved transcriptomics, snRNA-seq, bulk RNA-seq (across six library/extraction RNA-seq combinations), and RNAScope/Immunofluorescence (RNAScope/IF) for six broad cell types. The Mean Ratio method, implemented in the DeconvoBuddies R package, was developed for selecting cell type marker genes. Six computational deconvolution algorithms were evaluated in DLPFC and predicted cell type proportions were compared to orthogonal RNAScope/IF measurements. Conclusions Bisque and hspe were the most accurate methods, were robust to differences in RNA library types and extractions. This multi-assay dataset showed that cell size differences, marker genes differentially quantified across RNA libraries, and cell composition variability in reference snRNA-seq impact the accuracy of current deconvolution methods.

Norepinephrine (NE) neurons in the locus coeruleus (LC) make long-range projections throughout the central nervous system, playing critical roles in arousal and mood, as well as various components of cognition including attention, learning, and memory. The LC-NE system is also implicated in multiple neurological and neuropsychiatric disorders. Importantly, LC-NE neurons are highly sensitive to degeneration in both Alzheimer’s and Parkinson’s disease. Despite the clinical importance of the brain region and the prominent role of LC-NE neurons in a variety of brain and behavioral functions, a detailed molecular characterization of the LC is lacking. Here, we used a combination of spatially-resolved transcriptomics and single-nucleus RNA-sequencing to characterize the molecular landscape of the LC region and the transcriptomic profile of LC-NE neurons in the human brain. We provide a freely accessible resource of these data in web-accessible and downloadable formats.

The lateral septum (LS), a GABAergic structure located in the basal forebrain, is implicated in social behavior, learning, and memory. We previously demonstrated that expression of tropomyosin kinase receptor B (TrkB) in LS neurons is required for social novelty recognition. To better understand molecular mechanisms by which TrkB signaling controls behavior, we locally knocked down TrkB in LS and used bulk RNA-sequencing to identify changes in gene expression downstream of TrkB. TrkB knockdown induces upregulation of genes associated with inflammation and immune responses, and downregulation of genes associated with synaptic signaling and plasticity. Next, we generated one of the first atlases of molecular profiles for LS cell types using single nucleus RNA-sequencing (snRNA-seq). We identified markers for the septum broadly, and the LS specifically, as well as for all neuronal cell types. We then investigated whether the differentially expressed genes (DEGs) induced by TrkB knockdown map to specific LS cell types. Enrichment testing identified that downregulated DEGs are broadly expressed across neuronal clusters. Enrichment analyses of these DEGs demonstrated that downregulated genes are uniquely expressed in the LS, and associated with either synaptic plasticity or neurodevelopmental disorders. Upregulated genes are enriched in LS microglia, associated with immune response and inflammation, and linked to both neurodegenerative disease and neuropsychiatric disorders. In addition, many of these genes are implicated in regulating social behaviors. In summary, the findings implicate TrkB signaling in the LS as a critical regulator of gene networks associated with psychiatric disorders that display social deficits, including schizophrenia and autism, and with neurodegenerative diseases, including Alzheimer’s.

Deconvolution of cell mixtures in “bulk” transcriptomic samples from homogenate human tissue is important for understanding disease pathologies. However, several experimental and computational challenges impede transcriptomics-based deconvolution approaches using single-cell/nucleus RNA-seq reference atlases. Cells from the brain and blood have substantially different sizes, total mRNA, and transcriptional activities, and existing approaches may quantify total mRNA instead of cell type proportions. Further, standards are lacking for the use of cell reference atlases and integrative analyses of single-cell and spatial transcriptomics data. We discuss how to approach these key challenges with orthogonal “gold standard” datasets for evaluating deconvolution methods.

The creation of effective visualizations is a fundamental component of data analysis. In biomedical research, new challenges are emerging to visualize multi-dimensional data in a 2D space, but current data visualization tools have limited capabilities. To address this problem, we leverage Gestalt principles to improve the design and interpretability of multi-dimensional data in 2D data visualizations, layering aesthetics to display multiple variables. The proposed visualization can be applied to spatially-resolved transcriptomics data, but also broadly to data visualized in 2D space, such as embedding visualizations. We provide an open source R package escheR, which is built off of the state-of-the-art ggplot2 visualization framework and can be seamlessly integrated into genomics toolboxes and workflows.

The dentate gyrus of the hippocampus is important for many cognitive functions, including learning, memory, and mood. Here, we investigated age-associated changes in transcriptome-wide spatial gene expression in the human dentate gyrus across the lifespan. Genes associated with neurogenesis and the extracellular matrix were enriched in infants, while gene markers of inhibitory neurons and cell proliferation showed increases and decreases in post-infancy, respectively. While we did not find evidence for neural proliferation post-infancy, we did identify molecular signatures supporting protracted maturation of granule cells. We also identified a wide-spread hippocampal aging signature and an age-associated increase in genes related to neuroinflammation. Our findings suggest major changes to the putative neurogenic niche after infancy and identify molecular foci of brain aging in glial and neuropil enriched tissue.

Molecular mechanisms of neuropsychiatric disorders are challenging to study in human brain. For decades, the preferred model has been to study postmortem human brain samples despite the limitations they entail. A recent study generated RNA sequencing data from biopsies of prefrontal cortex from living patients with Parkinson’s Disease and compared gene expression to postmortem tissue samples, from which they found vast differences between the two. This led the authors to question the utility of postmortem human brain studies. Through re-analysis of the same data, we unexpectedly found that the living brain tissue samples were of much lower quality than the postmortem samples across multiple standard metrics. We also performed simulations that illustrate the effects of ignoring RNA degradation in differential gene expression analyses, showing the effects can be substantial and of similar magnitude to what the authors find. For these reasons, we believe the authors’ conclusions are unjustified. To the contrary, while opportunities to study gene expression in the living brain are welcome, evidence that this eclipses the value of postmortem analyses is not apparent.

Neuropathological lesions in the brains of individuals affected with neurodegenerative disorders are hypothesized to trigger molecular and cellular processes that disturb the homeostasis of local microenvironments. Here, we applied the 10x Genomics Visium Spatial Proteogenomics (Visium-SPG) platform, which couples spatial gene expression with immunofluorescence (IF) protein co-detection, to evaluate its ability to quantify changes in spatial gene expression with respect to amyloid-beta (Aβ) and hyperphosphorylated tau (pTau) pathology in post-mortem human brain tissue from individuals with Alzheimer’s disease (AD). We identified transcriptomic signatures associated with proximity to Aβ in the human inferior temporal cortex during late-stage AD, which we further investigated at cellular resolution with combined IF and single-molecule fluorescent in situ hybridization (smFISH). The study provides a data analysis workflow for Visium-SPG, and the data represent a proof-of-principle for the power of multi-omic profiling in identifying changes in molecular dynamics that are spatially associated with pathology in the human brain. We provide the scientific community with web-based, interactive resources to access the datasets of the spatially resolved AD-related transcriptomes.

We define and identify a new class of control genes for next-generation sequencing called total RNA expression genes (TREGs), which correlate with total RNA abundance in cell types of different sizes and transcriptional activity. We provide a data-driven method to identify TREGs from single-cell RNA sequencing data, allowing the estimation of total amount of RNA when restricted to quantifying a limited number of genes. We demonstrate our method in postmortem human brain using multiplex single-molecule fluorescent in situ hybridization and compare candidate TREGs against classic housekeeping genes. We identify AKT3 as a top TREG across five brain regions.

Activity-regulated gene (ARG) expression patterns in the hippocampus (HPC) regulate synaptic plasticity, learning, and memory, and are linked to both risk and treatment responses for many neuropsychiatric disorders. The HPC contains discrete classes of neurons with specialized functions, but cell type-specific activity-regulated transcriptional programs are not well characterized. Here, we used single-nucleus RNA-sequencing (snRNA-seq) in a mouse model of acute electroconvulsive seizures (ECS) to identify cell type-specific molecular signatures associated with induced activity in HPC neurons. We used unsupervised clustering and a priori marker genes to computationally annotate 15,990 high-quality HPC neuronal nuclei from N = 4 mice across all major HPC subregions and neuron types. Activity-induced transcriptomic responses were divergent across neuron populations, with dentate granule cells being particularly responsive to activity. Differential expression analysis identified both upregulated and downregulated cell type-specific gene sets in neurons following ECS. Within these gene sets, we identified enrichment of pathways associated with varying biological processes such as synapse organization, cellular signaling, and transcriptional regulation. Finally, we used matrix factorization to reveal continuous gene expression patterns differentially associated with cell type, ECS, and biological processes. This work provides a rich resource for interrogating activity-regulated transcriptional responses in HPC neurons at single-nuclei resolution in the context of ECS, which can provide biological insight into the roles of defined neuronal subtypes in HPC function.

Alternative splicing impacts most multi-exonic human genes. Inaccuracies during this process may have an important role in ageing and disease. Here, we investigated mis-splicing using RNA-sequencing data from ~14K control samples and 42 human body sites, focusing on split reads partially mapping to known transcripts in annotation. We show that mis-splicing occurs at different rates across introns and tissues and that these splicing inaccuracies are primarily affected by the abundance of core components of the spliceosome assembly and its regulators. Using publicly available data on short-hairpin RNA-knockdowns of numerous spliceosomal components and related regulators, we found support for the importance of RNA-binding proteins in mis-splicing. We also demonstrated that age is positively correlated with mis-splicing, and it affects genes implicated in neurodegenerative diseases. This in-depth characterisation of mis-splicing can have important implications for our understanding of the role of splicing inaccuracies in human disease and the interpretation of long-read RNA-sequencing data.

Most studies of gene expression in the brains of individuals with schizophrenia have focused on cortical regions, but subcortical nuclei such as the striatum are prominently implicated in the disease, and current antipsychotic drugs target the striatum’s dense dopaminergic innervation. Here, we performed a comprehensive analysis of the genetic and transcriptional landscape of schizophrenia in the postmortem caudate nucleus of the striatum of 443 individuals (245 neurotypical individuals, 154 individuals with schizophrenia and 44 individuals with bipolar disorder), 210 from African and 233 from European ancestries. Integrating expression quantitative trait loci analysis, Mendelian randomization with the latest schizophrenia genome-wide association study, transcriptome-wide association study and differential expression analysis, we identified many genes associated with schizophrenia risk, including potentially the dopamine D2 receptor short isoform. We found that antipsychotic medication has an extensive influence on caudate gene expression. We constructed caudate nucleus gene expression networks that highlight interactions involving schizophrenia risk. These analyses provide a resource for the study of schizophrenia and insights into risk mechanisms and potential therapeutic targets.

Following the recent review by Merikangas et al. summarizing differential gene expression in schizophrenia, we wanted to provide readers with our perspective on the state of the field in order to highlight current resources and ongoing challenges. We have led data collection and analysis of post mortem brains of control and schizophrenia samples across multiple large-scale cohorts, and we have provided data and results resources for each of these studies (Table 1).

Poster presented at BoG2022

Summary SpatialExperiment is a new data infrastructure for storing and accessing spatially resolved transcriptomics data, implemented within the R/Bioconductor framework, which provides advantages of modularity, interoperability, standardized operations, and comprehensive documentation. Here, we demonstrate the structure and user interface with examples from the 10x Genomics Visium and seqFISH platforms, and provide access to example datasets and visualization tools in the STexampleData, TENxVisiumData, and ggspavis packages. Availability and Implementation The SpatialExperiment, STexampleData, TENxVisiumData, and ggspavis packages are available from Bioconductor. The package versions described in this manuscript are available in Bioconductor version 3.15 onwards. Supplementary information Supplementary tables and figures are available at Bioinformatics online.

Objective: The authors sought to study the transcriptomic and genomic features of completed suicide by parsing the method chosen, to capture molecular correlates of the distinctive frame of mind of individuals who die by suicide, while reducing heterogeneity. Methods: The authors analyzed gene expression (RNA sequencing) from postmortem dorsolateral prefrontal cortex of patients who died by suicide with violent compared with nonviolent means, nonsuicide patients with the same psychiatric disorders, and a neurotypical group (total N=329). They then examined genomic risk scores (GRSs) for each psychiatric disorder included, and GRSs for cognition (IQ) and for suicide attempt, testing how they predict diagnosis or traits (total N=888). Results: Patients who died by suicide by violent means showed a transcriptomic pattern remarkably divergent from each of the other patient groups but less from the neurotypical group; consistently, their genomic profile of risk was relatively low for their diagnosed illness as well as for suicide attempt, and relatively high for IQ: they were more similar to the neurotypical group than to other patients. Differentially expressed genes (DEGs) associated with patients who died by suicide by violent means pointed to purinergic signaling in microglia, showing similarities to a genome-wide association study of Drosophila aggression. Weighted gene coexpression network analysis revealed that these DEGs were coexpressed in a context of mitochondrial metabolic activation unique to suicide by violent means. Conclusions: These findings suggest that patients who die by suicide by violent means are in part biologically separable from other patients with the same diagnoses, and their behavioral outcome may be less dependent on genetic risk for conventional psychiatric disorders and be associated with an alteration of purinergic signaling and mitochondrial metabolism.

Recent genetic studies have identified variants associated with bipolar disorder (BD), but it remains unclear how brain gene expression is altered in BD and how genetic risk for BD may contribute to these alterations. Here, we obtained transcriptomes from subgenual anterior cingulate cortex and amygdala samples from post-mortem brains of individuals with BD and neurotypical controls, including 511 total samples from 295 unique donors. We examined differential gene expression between cases and controls and the transcriptional effects of BD-associated genetic variants. We found two coexpressed modules that were associated with transcriptional changes in BD: one enriched for immune and inflammatory genes and the other with genes related to the postsynaptic membrane. Over 50% of BD genome-wide significant loci contained significant expression quantitative trait loci (QTL) (eQTL), and these data converged on several individual genes, including SCN2A and GRIN2A. Thus, these data implicate specific genes and pathways that may contribute to the pathology of BD.

Single-cell gene expression technologies are powerful tools to study cell types in the human brain, but efforts have largely focused on cortical brain regions. We therefore created a single-nucleus RNA-sequencing resource of 70,615 high-quality nuclei to generate a molecular taxonomy of cell types across five human brain regions that serve as key nodes of the human brain reward circuitry: nucleus accumbens, amygdala, subgenual anterior cingulate cortex, hippocampus, and dorsolateral prefrontal cortex. We first identified novel subpopulations of interneurons and medium spiny neurons (MSNs) in the nucleus accumbens and further characterized robust GABAergic inhibitory cell populations in the amygdala. Joint analyses across the 107 reported cell classes revealed cell-type substructure and unique patterns of transcriptomic dynamics. We identified discrete subpopulations of D1- and D2-expressing MSNs in the nucleus accumbens to which we mapped cell-type-specific enrichment for genetic risk associated with both psychiatric disease and addiction.

DNA methylation (DNAm) is an epigenetic regulator of gene expression and a hallmark of gene-environment interaction. Using whole-genome bisulfite sequencing, we have surveyed DNAm in 344 samples of human postmortem brain tissue from neurotypical subjects and individuals with schizophrenia. We identify genetic influence on local methylation levels throughout the genome, both at CpG sites and CpH sites, with 86% of SNPs and 55% of CpGs being part of methylation quantitative trait loci (meQTLs). These associations can further be clustered into regions that are differentially methylated by a given SNP, highlighting the genes and regions with which these loci are epigenetically associated. These findings can be used to better characterize schizophrenia GWAS-identified variants as epigenetic risk variants. Regions differentially methylated by schizophrenia risk-SNPs explain much of the heritability associated with risk loci, despite covering only a fraction of the genomic space. We provide a comprehensive, single base resolution view of association between genetic variation and genomic methylation, and implicate schizophrenia GWAS-associated variants as influencing the epigenetic plasticity of the brain.

Although next-generation sequencing technologies have accelerated the discovery of novel gene-to-disease associations, many patients with suspected Mendelian diseases still leave the clinic without a genetic diagnosis. An estimated one third of these patients will have disorders caused by mutations impacting splicing. RNA-sequencing has been shown to be a promising diagnostic tool, however few methods have been developed to integrate RNA-sequencing data into the diagnostic pipeline. Here, we introduce dasper, an R/Bioconductor package that improves upon existing tools for detecting aberrant splicing by using machine learning to incorporate disruptions in exon-exon junction counts as well as coverage. dasper is designed for diagnostics, providing a rank-based report of how aberrant each splicing event looks, as well as including visualization functionality to facilitate interpretation. We validate dasper using 16 patient-derived fibroblast cell lines harbouring pathogenic variants known to impact splicing. We find that dasper is able to detect pathogenic splicing events with greater accuracy than existing LeafCutterMD or z-score approaches. Furthermore, by only applying a broad OMIM gene filter (without any variant-level filters), dasper is able to detect pathogenic splicing events within the top 10 most aberrant identified for each patient. Since using publicly available control data minimises costs associated with incorporating RNA-sequencing into diagnostic pipelines, we also investigate the use of 504 GTEx fibroblast samples as controls. We find that dasper leverages publicly available data effectively, ranking pathogenic splicing events in the top 25. Thus, we believe dasper can increase diagnostic yield for a pathogenic splicing variants and enable the efficient implementation of RNA-sequencing for diagnostics in clinical laboratories.

Background. Recent breakthroughs in psychiatric genetics have implicated biological pathways onto which genetic risk for psychiatric disorders converges. However, these studies do not reveal the developmental time point(s) at which these pathways are relevant. Methods. We aimed to determine the relationship between psychiatric risk and developmental gene expression relating to discrete biological pathways. We used postmortem RNA sequencing data (BrainSeq and BrainSpan) from brain tissue at multiple prenatal and postnatal time points, with summary statistics from recent genome-wide association studies of schizophrenia, bipolar disorder, and major depressive disorder. We prioritized gene sets for overall enrichment of association with each disorder and then tested the relationship between the association of their constituent genes with their relative expression at each developmental stage. Results. We observed relationships between the expression of genes involved in voltage-gated cation channel activity during early midfetal, adolescence, and early adulthood time points and association with schizophrenia and bipolar disorder, such that genes more strongly associated with these disorders had relatively low expression during early midfetal development and higher expression during adolescence and early adulthood. The relationship with schizophrenia was strongest for the subset of genes related to calcium channel activity, while for bipolar disorder, the relationship was distributed between calcium and potassium channel activity genes. Conclusions. Our results indicate periods during development when biological pathways related to the activity of calcium and potassium channels may be most vulnerable to the effects of genetic variants conferring risk for psychiatric disorders. Furthermore, they indicate key time points and potential targets for disorder-specific therapeutic interventions.

Poster on the HumanPilot project for BioC2020 using spatialLIBD.

DNA methylation (DNAm) is a key epigenetic regulator of gene expression across development. The developing prenatal brain is a highly dynamic tissue, but our understanding of key drivers of epigenetic variability across development is limited. We, therefore, assessed genomic methylation at over 39 million sites in the prenatal cortex using whole-genome bisulfite sequencing and found loci and regions in which methylation levels are dynamic across development. We saw that DNAm at these loci was associated with nearby gene expression and enriched for enhancer chromatin states in prenatal brain tissue. Additionally, these loci were enriched for genes associated with neuropsychiatric disorders and genes involved with neurogenesis. We also found autosomal differences in DNAm between the sexes during prenatal development, though these have less clear functional consequences. We lastly confirmed that the dynamic methylation at this critical period is specifically CpG methylation, with generally low levels of CpH methylation. Our findings provide detailed insight into prenatal brain development as well as clues to the pathogenesis of psychiatric traits seen later in life.

SynGAP is a synaptic Ras GTPase-activating protein (GAP) with four C-terminal splice variants: α1, α2, β, and γ. Although studies have implicated SYNGAP1 in several cognitive disorders, it is not clear which SynGAP isoforms contribute to disease. Here, we demonstrate that SynGAP isoforms exhibit unique spatiotemporal expression patterns and play distinct roles in neuronal and synaptic development in mouse neurons. SynGAP-α1, which undergoes liquid-liquid phase separation with PSD-95, is highly enriched in synapses and is required for LTP. In contrast, SynGAP-β, which does not bind PSD-95 PDZ domains, is less synaptically targeted and promotes dendritic arborization. A mutation in SynGAP-α1 that disrupts phase separation and synaptic targeting abolishes its ability to regulate plasticity and instead causes it to drive dendritic development like SynGAP-β. These results demonstrate that distinct intrinsic biochemical properties of SynGAP isoforms determine their function, and individual isoforms may differentially contribute to the pathogenesis of SYNGAP1-related cognitive disorders.

RegulonDB has collected, harmonized and centralized data from hundreds of experiments for nearly two decades and is considered a point of reference for transcriptional regulation in Escherichia coli K12. Here, we present the regutools R package to facilitate programmatic access to RegulonDB data in computational biology. regutools gives researchers the possibility of writing reproducible workflows with automated queries to RegulonDB. The regutools package serves as a bridge between RegulonDB data and the Bioconductor ecosystem by reusing the data structures and statistical methods powered by other Bioconductor packages. We demonstrate the integration of regutools with Bioconductor by analyzing transcription factor DNA binding sites and transcriptional regulatory networks from RegulonDB. We anticipate that regutools will serve as a useful building block in our progress to further our understanding of gene regulatory networks.

Jaffe et al. profile the granule cell layer of the human hippocampus and find unique molecular associations for aging and genetic variation, as well as diagnosis with schizophrenia and its genetic risk, that were previously undiscovered in homogenate tissue.

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 enhancers 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.

Human induced pluripotent stem cells (hiPSCs) are a powerful model of neural differentiation and maturation. We present a hiPSC transcriptomics resource on corticogenesis from 5 iPSC donor and 13 subclonal lines across 9 time points over 5 broad conditions: self-renewal, early neuronal differentiation, neural precursor cells (NPCs), assembled rosettes, and differentiated neuronal cells. We identify widespread changes in the expression of both individual features and global patterns of transcription. We next demonstrate that co-culturing human NPCs with rodent astrocytes results in mutually synergistic maturation, and that cell type-specific expression data can be extracted using only sequencing read alignments without cell sorting. We lastly adapt a previously generated RNA deconvolution approach to single-cell expression data to estimate the relative neuronal maturity of iPSC-derived neuronal cultures and human brain tissue. Using many public datasets, we demonstrate neuronal cultures are maturationally heterogeneous but contain subsets of neurons more mature than previously observed.

Background: DNA methylation (DNAm) is a critical regulator of both development and cellular identity and shows unique patterns in neurons. To better characterize maturational changes in DNAm patterns in these cells, we profile the DNAm landscape at single-base resolution across the first two decades of human neocortical development in NeuN+ neurons using whole-genome bisulfite sequencing and compare them to non-neurons (primarily glia) and prenatal homogenate cortex. Results: We show that DNAm changes more dramatically during the first 5 years of postnatal life than during the entire remaining period. We further refine global patterns of increasingly divergent neuronal CpG and CpH methylation (mCpG and mCpH) into six developmental trajectories and find that in contrast to genome-wide patterns, neighboring mCpG and mCpH levels within these regions are highly correlated. We integrate paired RNAseq data and identify putative regulation of hundreds of transcripts and their splicing events exclusively by mCpH levels, independently from mCpG levels, across this period. We finally explore the relationship between DNAm patterns and development of brain-related phenotypes and find enriched heritability for many phenotypes within identified DNAm features. Conclusions: By profiling DNAm changes in NeuN-sorted neurons over the span of human cortical development, we identify novel, dynamic regions of DNAm that would be masked in homogenate DNAm data; expand on the relationship between CpG methylation, CpH methylation, and gene expression; and find enrichment particularly for neuropsychiatric diseases in genomic regions with cell type-specific, developmentally dynamic DNAm patterns. Keywords: DNA methylation, Neurodevelopment, Gene expression, Non-CpG methylation.

Background: RNA sequencing offers advantages over other quantification methods for microRNA (miRNA), yet numerous biases make reliable quantification challenging. Previous evaluations of these biases have focused on adapter ligation bias with limited evaluation of reverse transcription bias or amplification bias. Furthermore, evaluations of the quantification of isomiRs (miRNA isoforms) or the influence of starting amount on performance have been very limited. No study had yet evaluated the quantification of isomiRs of altered length or compared the consistency of results derived from multiple moderate starting inputs. We therefore evaluated quantifications of miRNA and isomiRs using four library preparation kits, with various starting amounts, as well as quantifications following removal of duplicate reads using unique molecular identifiers (UMIs) to mitigate reverse transcription and amplification biases. Results: All methods resulted in false isomiR detection; however, the adapter-free method tested was especially prone to false isomiR detection. We demonstrate that using UMIs improves accuracy and we provide a guide for input amounts to improve consistency. Conclusions: Our data show differences and limitations of current methods, thus raising concerns about the validity of quantification of miRNA and isomiRs across studies. We advocate for the use of UMIs to improve accuracy and reliability of miRNA quantifications.

The hippocampus formation, although prominently implicated in schizophrenia pathogenesis, has been overlooked in large-scale genomics efforts in the schizophrenic brain. We performed RNA-seq in hippocampi and dorsolateral prefrontal cortices (DLPFCs) from 551 individuals (286 with schizophrenia). We identified substantial regional differences in gene expression and found widespread developmental differences that were independent of cellular composition. We identified 48 and 245 differentially expressed genes (DEGs) associated with schizophrenia within the hippocampus and DLPFC, with little overlap between the brain regions. 124 of 163 (76.6%) of schizophrenia GWAS risk loci contained eQTLs in any region. Transcriptome-wide association studies in each region identified many novel schizophrenia risk features that were brain region-specific. Last, we identified potential molecular correlates of in vivo evidence of altered prefrontal-hippocampal functional coherence in schizophrenia. These results underscore the complexity and regional heterogeneity of the transcriptional correlates of schizophrenia and offer new insights into potentially causative biology.

Poster presented at BoG2019

The usability of publicly-available gene expression data is often limited by the availability of high-quality, standardized biological phenotype and experimental condition information (“metadata”). We released the recount2 project, which involved re-processing ∼70,000 samples in the Sequencing Read Archive (SRA), Genotype-Tissue Expression (GTEx), and The Cancer Genome Atlas (TCGA) projects. While samples from the latter two projects are well-characterized with extensive metadata, the ∼50,000 RNA-seq samples from SRA in recount2 are inconsistently annotated with metadata. Tissue type, sex, and library type can be estimated from the RNA sequencing (RNA-seq) data itself. However, more detailed and harder to predict metadata, like age and diagnosis, must ideally be provided by labs that deposit the data. To facilitate more analyses within human brain tissue data, we have complemented phenotype predictions by manually constructing a uniformly-curated database of public RNA-seq samples present in SRA and recount2. We describe the reproducible curation process for constructing recount-brain that involves systematic review of the primary manuscript, which can serve as a guide to annotate other studies and tissues. We further expanded recount-brain by merging it with GTEx and TCGA brain samples as well as linking to controlled vocabulary terms for tissue, Brodmann area and disease. Furthermore, we illustrate how to integrate the sample metadata in recount-brain with the gene expression data in recount2 to perform differential expression analysis. We then provide three analysis examples involving modeling postmortem interval, glioblastoma, and meta-analyses across GTEx and TCGA. Overall, recount-brain facilitates expression analyses and improves their reproducibility as individual researchers do not have to manually curate the sample metadata. recount-brain is available via the add_metadata() function from the recount Bioconductor package at bioconductor.org/packages/recount.

Understanding the molecular profile of every human cell type is essential for understanding its role in normal physi-ology and disease. Technological advancements in DNA sequencing, mass spectrometry, and computational methodsallow us to carry out multiomics analyses although such approaches are not routine yet. Human umbilical vein endothe-lial cells (HUVECs) are a widely used model system to study pathological and physiological processes associated withthe cardiovascular system. In this study, next-generation sequencing and high-resolution mass spectrometry to profilethe transcriptome and proteome of primary HUVECs is employed. Analysis of 145 million paired-end reads from next-generation sequencing confirmed expression of 12 186 protein-coding genes (FPKM >= 0.1), 439 novel long non-codingRNAs, and revealed 6089 novel isoforms that were not annotated in GENCODE. Proteomics analysis identifies 6477proteins including confirmation ofN-termini for 1091 proteins, isoforms for 149 proteins, and 1034 phosphosites. Adatabase search to specifically identify other post-translational modifications provide evidence for a number of modifi-cation sites on 117 proteins which include ubiquitylation, lysine acetylation, and mono-, di- and tri-methylation events.Evidence for 11 ‘missing proteins’, which are proteins for which there was insufficient or no protein level evidence, isprovided. Peptides supporting missing protein and novel events are validated by comparison of MS/MS fragmentationpatterns with synthetic peptides. Finally, 245 variant peptides derived from 207 expressed proteins in addition to alternatetranslational start sites for seven proteins and evidence for novel proteoforms for five proteins resulting from alternativesplicing are identified. Overall, it is believed that the integrated approach employed in this study is widely applicable tostudy any primary cell type for deeper molecular characterization.

Late-onset Alzheimer’s disease (AD) is a complex age-related neurodegenerative disorder that likely involves epigenetic factors. To better understand the epigenetic state associated with AD, we surveyed 420,852 DNA methylation (DNAm) sites from neurotypical controls (N=49) and late-onset AD patients (N=24) across four brain regions (hippocampus, entorhinal cortex, dorsolateral prefrontal cortex and cerebellum). We identified 858 sites with robust differential methylation collectively annotated to 772 possible genes (FDR<5%, within 10 kb). These sites were overrepresented in AD genetic risk loci (p=0.00655) and were enriched for changes during normal aging (p<2.2×10−16), and nearby genes were enriched for processes related to cell-adhesion, immunity, and calcium homeostasis (FDR<5%). To functionally validate these associations, we generated and analyzed corresponding transcriptome data to prioritize 130 genes within 10 kb of the differentially methylated sites. These 130 genes were differentially expressed between AD cases and controls and their expression was associated with nearby DNAm (p<0.05). This integrated analysis implicates novel genes in Alzheimer’s disease, such as ANKRD30B. These results highlight DNAm differences in Alzheimer’s disease that have gene expression correlates, further implicating DNAm as an epigenetic mechanism underlying pathological molecular changes associated with AD. Furthermore, our framework illustrates the value of integrating epigenetic and transcriptomic data for understanding complex disease.

Background: Antibody class switch recombination (CSR) to IgG, IgA or IgE is a hallmark of adaptive immunity, allowing antibody function diversification beyond IgM. CSR involves a deletion of the IgM/IgD constant region genes placing a new acceptor Constant (CH) gene, downstream of the VDJH exon. CSR depends on non-coding (CSRnc) transcription of donor Iμ and acceptor IH exons, located 5′ upstream of each CH coding gene. Although our knowledge of the role of CSRnc transcription has advanced greatly, its extension and importance in healthy and diseased humans is scarce. Methods: We analyzed CSRnc transcription in 70,603 publicly available RNA-seq samples, including GTEx, TCGA and the Sequence Read Archive (SRA) using recount2, an online resource consisting of normalized RNA-seq gene and exon counts, as well as coverage BigWig files that can be programmatically accessed through R. CSRnc transcription was validated with a qRT-PCR assay for Iμ, Iγ1 and Iγ3 in humans in response to vaccination. Results: We mapped IH transcription for the human IgH locus, including the less understood IGHD gene. CSRnc transcription was restricted to B cells and is widely distributed in normal adult tissues, but predominant in blood, spleen, MALT-containing tissues, visceral adipose tissue and some so-called ‘immune privileged’ tissues. However, significant Iγ4 expression was found even in non-lymphoid fetal tissues. CSRnc expression in cancer tissues mimicked the expression of their normal counterparts, with notable pattern changes in some common cancer subsets. CSRnc transcription in tumors appears to result from tumor infiltration by B cells, since CSRnc transcription was not detected in corresponding tumor-derived immortal cell lines. Additionally, significantly increased Iδ transcription in ileal mucosa in Crohn’s disease with ulceration was found. Conclusions: CSRnc transcription occurs in multiple anatomical locations beyond classical secondary lymphoid organs, representing a potentially useful marker of effector B cell responses in normal and pathological immune responses. The pattern of IH exon expression may reveal clues of the local immune response (i.e. cytokine milieu) in health and disease. This is a great example of how the public recount2 data can be used to further our understanding of transcription, including regions outside the known transcriptome.

Poster presented at ASHG2018, biodata18 and PSB2019

This paper is on smoking and its relation to gene expression at different life stages. These genes are also related to autism spectrum disorder.

Poster presented at SAGES2018

More than 70,000 short-read RNA-sequencing samples are publicly available through the recount2 project, a curated database of summary coverage data. However, no current methods can estimate transcript-level abundances using the reduced-representation information stored in this database. Here we present a linear model utilizing coverage of junctions and subdivided exons to generate transcript abundance estimates of comparable accuracy to those obtained from methods requiring read-level data. Our approach flexibly models bias, produces standard errors, and is easy to refresh given updated annotation. We illustrate our method on simulated and real data and release transcript abundance estimates for the samples in recount2.

Poster presented at IDIES2017

The recount2 resource is composed of over 70,000 uniformly processed human RNA-seq samples spanning TCGA and SRA, including GTEx. The processed data can be accessed via the recount2 website and the recount Bioconductor package. This workflow explains in detail how to use the recount package and how to integrate it with other Bioconductor packages for several analyses that can be carried out with the recount2 resource. In particular, we describe how the coverage count matrices were computed in recount2 as well as different ways of obtaining public metadata, which can facilitate downstream analyses. Step-by-step directions show how to do a gene-level differential expression analysis, visualize base-level genome coverage data, and perform an analyses at multiple feature levels. This workflow thus provides further information to understand the data in recount2 and a compendium of R code to use the data.

Expression of the gene set of HNMT, HRH1, HRH2 and HRH3 was significantly altered between ASD and matched controls, and this finding was replicated with an independent data set.

derfinder analysis using expressed region-level and single base-level approaches provides a compromise between full transcript reconstruction and feature-level analysis. The package is available from Bioconductor.

We describe Rail-RNA, a cloud-enabled spliced aligner that analyzes many samples at once. Rail-RNA eliminates redundant work across samples, making it more efficient as samples are added. For many samples, Rail-RNA is more accurate than annotation-assisted aligners.

regionReport is an R package for generating detailed interactive reports from region-level genomic analyses as well as feature-level RNA-seq. The report includes quality-control checks, an overview of the results, an interactive table of the genomic regions or features of interest and reproducibility information. regionReport provides specialised reports for exploring DESeq2, edgeR, or derfinder differential expression analyses results. regionReport is also flexible and can easily be expanded with report templates for other analysis pipelines.

Poster presented at ASHG2015 and IDIES2015

Transcriptome analysis of human brain provides fundamental insight into development and disease, but it largely relies on existing annotation. We sequenced transcriptomes of 72 prefrontal cortex samples across six life stages and identified 50,650 differentially expression regions (DERs) associated with developmental and aging, agnostic of annotation. While many DERs annotated to non-exonic sequence (41.1%), most were similarly regulated in cytosolic mRNA extracted from independent samples. The DERs were developmentally conserved across 16 brain regions and in the developing mouse cortex, and were expressed in diverse cell and tissue types. The DERs were further enriched for active chromatin marks and clinical risk for neurodevelopmental disorders such as schizophrenia. Lastly, we demonstrate quantitatively that these DERs associate with a changing neuronal phenotype related to differentiation and maturation. These data show conserved molecular signatures of transcriptional dynamics across brain development, have potential clinical relevance and highlight the incomplete annotation of the human brain transcriptome.

There has been a major shift from microarrays to RNA-sequencing (RNA-seq) for measuring gene expression as the price per measurement between these technologies has become comparable. The advantages of RNA-seq are increased measurement flexibility to detect alternative transcription, allele specific transcription, or transcription outside of known coding regions. The price of this increased flexibility is: (a) an increase in raw data size and (b) more decisions that must be made by the data analyst. Here we provide a selective review and extension of our previous work in attempting to measure variability in results due to different choices about how to summarize and analyze RNA-sequencing data. We discuss a standard model for gene expression measurements that breaks variability down into variation due to technology, biology, and measurement error. Finally, wee show the importance of gene model selection, normalization, and choice for statistical model on the ultimate results of an RNA-sequencing experiment.

Poster presented at ENAR2014

Poster presented at Genomics and Bioinformatics Symposium 2013, Delta Omega Poster Competition

Poster presented at JHU Biostatistics Department Retreat

Many different systems of bacterial interactions have been described. However, relatively few studies have explored how interactions between different microorganisms might influence bacterial development. To explore such interspecies interactions, we focused on Bacillus subtilis, which characteristically develops into matrix-producing cannibals before entering sporulation. We investigated whether organisms from the natural environment of B. subtilis—the soil—were able to alter the development of B.subtilis. To test this possibility, we developed a coculture microcolony screen in which we used fluorescent reporters to identify soil bacteria able to induce matrix production in B. subtilis. Most of the bacteria that influence matrix production in B. subtilis are members of the genus Bacillus, suggesting that such interactions may be predominantly with close relatives. The interactions we observed were mediated via two different mechanisms. One resulted in increased expression of matrix genes via the activation of a sensor histidine kinase, KinD. The second was kinase independent and conceivably functions by altering the relative subpopulations of B. subtilis cell types by preferentially killing noncannibals. These two mechanisms were grouped according to the inducing strain’s relatedness to B. subtilis. Our results suggest that bacteria preferentially alter their development in response to secreted molecules from closely related bacteria and do so using mechanisms that depend on the phylogenetic relatedness of the interacting bacteria.

Poster presented at EMBO Conference Series from Functional Genomics to Systems Biology 2010

Poster presented at BioC2010