Linear models enable powerful differential activity analysis in massively parallel reporter assays

Decoding biology with massively parallel reporter assays and machine learning

Massively parallel reporter assays (MPRAs) are powerful tools for quantifying the impacts of sequence variation on gene expression. Reading out molecular phenotypes with sequencing enables interrogating the impact of sequence variation beyond genome scale. Machine learning models integrate and codify information learned from MPRAs and enable generalization by predicting sequences outside the training data set. Models can provide a quantitative understanding of cis-regulatory codes controlling gene expression, enable variant stratification, and guide the design of synthetic regulatory elements for applications from synthetic biology to mRNA and gene therapy. This review focuses on cis-regulatory MPRAs, particularly those that interrogate cotranscriptional and post-transcriptional processes: alternative splicing, cleavage and polyadenylation, translation, and mRNA decay.

Single Cell Sequencing Provides Clues about the Developmental Genetic Basis of Evolutionary Adaptations in Syngnathid Fishes

Seahorses, pipefishes, and seadragons are fishes from the family Syngnathidae that have evolved extraordinary traits including male pregnancy, elongated snouts, loss of teeth, and dermal bony armor. The developmental genetic and cellular changes that led to the evolution of these traits are largely unknown. Recent syngnathid genome assemblies revealed suggestive gene content differences and provide the opportunity for detailed genetic analyses. We created a single cell RNA sequencing atlas of Gulf pipefish embryos to understand the developmental basis of four traits: derived head shape, toothlessness, dermal armor, and male pregnancy. We completed marker gene analyses, built genetic networks, and examined spatial expression of select genes. We identified osteochondrogenic mesenchymal cells in the elongating face that express regulatory genes bmp4, sfrp1a, and prdm16. We found no evidence for tooth primordia cells, and we observed re-deployment of osteoblast genetic networks in developing dermal armor. Finally, we found that epidermal cells expressed nutrient processing and environmental sensing genes, potentially relevant for the brooding environment. The examined pipefish evolutionary innovations are composed of recognizable cell types, suggesting derived features originate from changes within existing gene networks. Future work addressing syngnathid gene networks across multiple stages and species is essential for understanding how their novelties evolved.

Deciphering the functional impact of Alzheimer's Disease-associated variants in resting and proinflammatory immune cells

Genome-wide association studies have identified loci associated with Alzheimer's Disease (AD), but identifying the exact causal variants and genes at each locus is challenging due to linkage disequilibrium and their largely non-coding nature. To address this, we performed a massively parallel reporter assay of 3,576 AD-associated variants in THP-1 macrophages in both resting and proinflammatory states and identified 47 expression-modulating variants (emVars). To understand the endogenous chromatin context of emVars, we built an activity-by-contact model using epigenomic maps of macrophage inflammation and inferred condition-specific enhancer-promoter pairs. Intersection of emVars with enhancer-promoter pairs and microglia expression quantitative trait loci allowed us to connect 39 emVars to 76 putative AD risk genes enriched for AD-associated molecular signatures. Overall, systematic characterization of AD-associated variants enhances our understanding of the regulatory mechanisms underlying AD pathogenesis.

Liver eQTL meta-analysis illuminates potential molecular mechanisms of cardiometabolic traits

Understanding the molecular mechanisms of complex traits is essential for developing targeted interventions. We analyzed liver expression quantitative-trait locus (eQTL) meta-analysis data on 1,183 participants to identify conditionally distinct signals. We found 9,013 eQTL signals for 6,564 genes; 23% of eGenes had two signals, and 6% had three or more signals. We then integrated the eQTL results with data from 29 cardiometabolic genome-wide association study (GWAS) traits and identified 1,582 GWAS-eQTL colocalizations for 747 eGenes. Non-primary eQTL signals accounted for 17% of all colocalizations. Isolating signals by conditional analysis prior to coloc resulted in 37% more colocalizations than using marginal eQTL and GWAS data, highlighting the importance of signal isolation. Isolating signals also led to stronger evidence of colocalization: among 343 eQTL-GWAS signal pairs in multi-signal regions, analyses that isolated the signals of interest resulted in higher posterior probability of colocalization for 41% of tests. Leveraging allelic heterogeneity, we predicted causal effects of gene expression on liver traits for four genes. To predict functional variants and regulatory elements, we colocalized eQTL with liver chromatin accessibility QTL (caQTL) and found 391 colocalizations, including 73 with non-primary eQTL signals and 60 eQTL signals that colocalized with both a caQTL and a GWAS signal. Finally, we used publicly available massively parallel reporter assays in HepG2 to highlight 14 eQTL signals that include at least one expression-modulating variant. This multi-faceted approach to unraveling the genetic underpinnings of liver-related traits could lead to therapeutic development.

Unsupervised discovery of family specific vocal usage in the Mongolian gerbil

In nature, animal vocalizations can provide crucial information about identity, including kinship and hierarchy. However, lab-based vocal behavior is typically studied during brief interactions between animals with no prior social relationship, and under environmental conditions with limited ethological relevance. Here, we address this gap by establishing long-term acoustic recordings from Mongolian gerbil families, a core social group that uses an array of sonic and ultrasonic vocalizations. Three separate gerbil families were transferred to an enlarged environment and continuous 20-day audio recordings were obtained. Using a variational autoencoder (VAE) to quantify 583,237 vocalizations, we show that gerbils exhibit a more elaborate vocal repertoire than has been previously reported and that vocal repertoire usage differs significantly by family. By performing gaussian mixture model clustering on the VAE latent space, we show that families preferentially use characteristic sets of vocal clusters and that these usage preferences remain stable over weeks. Furthermore, gerbils displayed family-specific transitions between vocal clusters. Since gerbils live naturally as extended families in complex underground burrows that are adjacent to other families, these results suggest the presence of a vocal dialect which could be exploited by animals to represent kinship. These findings position the Mongolian gerbil as a compelling animal model to study the neural basis of vocal communication and demonstrates the potential for using unsupervised machine learning with uninterrupted acoustic recordings to gain insights into naturalistic animal behavior.

Translation elongation as a rate limiting step of protein production

The impact of synonymous codon choice on protein output has important implications for understanding endogenous gene expression and design of synthetic mRNAs. Synonymous codons are decoded at different speeds, but simple models predict that this should not drive protein output. Instead, translation initiation should be the rate limiting step for production of protein per mRNA, with little impact of codon choice. Previously, we used a neural network model to design a series of synonymous fluorescent reporters and showed that their protein output in yeast spanned a seven-fold range corresponding to their predicted translation elongation speed. Here, we show that this effect is not due primarily to the established impact of slow elongation on mRNA stability, but rather, that slow elongation further decreases the number of proteins made per mRNA. We combine simulations and careful experiments on fluorescent reporters to show that translation is limited on non-optimally encoded transcripts. Using a genome-wide CRISPRi screen, we find that impairing translation initiation attenuates the impact of slow elongation, showing a dynamic balance between rate limiting steps of protein production. Our results show that codon choice can directly limit protein production across the full range of endogenous variability in codon usage.

Genetic susceptibility to neurodevelopmental conditions associates with neonatal DNA methylation patterns in the general population: an individual participant data meta-analysis

Background

Autism spectrum disorder (ASD), attention-deficit/hyperactivity disorder (ADHD), and schizophrenia (SCZ) are highly heritable and linked to disruptions in foetal (neuro)development. While epigenetic processes are considered an important underlying pathway between genetic susceptibility and neurodevelopmental conditions, it is unclear (i) whether genetic susceptibility to these conditions is associated with epigenetic patterns, specifically DNA methylation (DNAm), already at birth; (ii) to what extent DNAm patterns are unique or shared across conditions, and (iii) whether these neonatal DNAm patterns can be leveraged to enhance genetic prediction of (neuro)developmental outcomes.

Methods

We conducted epigenome-wide meta-analyses of genetic susceptibility to ASD, ADHD, and schizophrenia, quantified using polygenic scores (PGSs) on cord blood DNAm, using four population-based cohorts (n pooled=5,802), all North European. Heterogeneity statistics were used to estimate overlap in DNAm patterns between PGSs. Subsequently, DNAm-based measures of PGSs were built in a target sample, and used as predictors to test incremental variance explained over PGS in 130 (neuro)developmental outcomes spanning birth to 14 years.

Outcomes

In probe-level analyses, SCZ-PGS associated with neonatal DNAm at 246 loci (p<9×10-8), predominantly in the major histocompatibility complex. Functional characterization of these DNAm loci confirmed strong genetic effects, significant blood-brain concordance and enrichment for immune-related pathways. 8 loci were identified for ASD-PGS (mapping to FDFT1 and MFHAS1), and none for ADHD-PGS. Regional analyses indicated a large number of differentially methylated regions for all PGSs (SCZ-PGS: 157, ASD-PGS: 130, ADHD-PGS: 166). DNAm signals showed little overlap between PGSs. We found suggestive evidence that incorporating DNAm-based measures of genetic susceptibility at birth increases explained variance for several child cognitive and motor outcomes over and above PGS.

Interpretation

Genetic susceptibility for neurodevelopmental conditions, particularly schizophrenia, is detectable in cord blood DNAm at birth in a population-based sample, with largely distinct DNAm patterns between PGSs. These findings support an early-origins perspective on schizophrenia.

Funding

HorizonEurope; European Research Council.

Precise measurement of molecular phenotypes with barcode-based CRISPRi systems

Genome-wide CRISPR-Cas9 screens have untangled regulatory networks and revealed the genetic underpinnings of diverse biological processes. Their success relies on experimental designs that interrogate specific molecular phenotypes and distinguish key regulators from background effects. Here, we realize these goals with a generalizable platform for CRISPR interference with barcoded expression reporter sequencing (CiBER-seq) that dramatically improves the sensitivity and scope of genome-wide screens. We systematically address technical factors that distort phenotypic measurements by normalizing expression reporters against closely-matched control promoters, integrated together into the genome at single copy. To test our ability to capture post-transcriptional and post-translational regulation through sequencing, we screened for genes that affected nonsense-mediated mRNA decay and Doa10-mediated cytosolic protein decay. Our optimized CiBER-seq screens accurately capture the known components of well-studied RNA and protein quality control pathways with minimal background. These results demonstrate the precision and versatility of CiBER-seq for dissecting the genetic networks controlling cellular behaviors.

Defining the Regulatory Logic of Breast Cancer Using Single-Cell Epigenetic and Transcriptome Profiling

Annotation of the cis-regulatory elements that drive transcriptional dysregulation in cancer cells is critical to improving our understanding of tumor biology. Herein, we present a compendium of matched chromatin accessibility (scATAC-seq) and transcriptome (scRNA-seq) profiles at single-cell resolution from human breast tumors and healthy mammary tissues processed immediately following surgical resection. We identify the most likely cell-of-origin for luminal breast tumors and basal breast tumors and then introduce a novel methodology that implements linear mixed-effects models to systematically quantify associations between regions of chromatin accessibility (i.e. regulatory elements) and gene expression in malignant cells versus normal mammary epithelial cells. These data unveil regulatory elements with that switch from silencers of gene expression in normal cells to enhancers of gene expression in cancer cells, leading to the upregulation of clinically relevant oncogenes. To translate the utility of this dataset into tractable models, we generated matched scATAC-seq and scRNA-seq profiles for breast cancer cell lines, revealing, for each subtype, a conserved oncogenic gene expression program between in vitro and in vivo cells. Together, this work highlights the importance of non-coding regulatory mechanisms that underlie oncogenic processes and the ability of single-cell multi-omics to define the regulatory logic of BC cells at single-cell resolution.

Depletion of cap-binding protein eIF4E dysregulates amino acid metabolic gene expression

Protein synthesis is metabolically costly and must be tightly coordinated with changing cellular needs and nutrient availability. The cap-binding protein eIF4E makes the earliest contact between mRNAs and the translation machinery, offering a key regulatory nexus. We acutely depleted this essential protein and found surprisingly modest effects on cell growth and recovery of protein synthesis. Paradoxically, impaired protein biosynthesis upregulated genes involved in the catabolism of aromatic amino acids simultaneously with the induction of the amino acid biosynthetic regulon driven by the integrated stress response factor GCN4. We further identified the translational control of Pho85 cyclin 5 (PCL5), a negative regulator of Gcn4, that provides a consistent protein-to-mRNA ratio under varied translation environments. This regulation depended in part on a uniquely long poly(A) tract in the PCL5 5' UTR and poly(A) binding protein. Collectively, these results highlight how eIF4E connects protein synthesis to metabolic gene regulation, uncovering mechanisms controlling translation during environmental challenges.

Massively parallel characterization of regulatory elements in the developing human cortex

Nucleotide changes in gene regulatory elements are important determinants of neuronal development and diseases. Using massively parallel reporter assays in primary human cells from mid-gestation cortex and cerebral organoids, we interrogated the cis-regulatory activity of 102,767 open chromatin regions, including thousands of sequences with cell type-specific accessibility and variants associated with brain gene regulation. In primary cells, we identified 46,802 active enhancer sequences and 164 variants that alter enhancer activity. Activity was comparable in organoids and primary cells, suggesting that organoids provide an adequate model for the developing cortex. Using deep learning we decoded the sequence basis and upstream regulators of enhancer activity. This work establishes a comprehensive catalog of functional gene regulatory elements and variants in human neuronal development.

Combining DNA and protein alignments to improve genome annotation with LiftOn

As the number and variety of assembled genomes continues to grow, the number of annotated genomes is falling behind, particularly for eukaryotes. DNA-based mapping tools help to address this challenge, but they are only able to transfer annotation between closely-related species. Here we introduce LiftOn, a homology-based software tool that integrates DNA and protein alignments to enhance the accuracy of genome-scale annotation and to allow mapping between relatively distant species. LiftOn's protein-centric algorithm considers both types of alignments, chooses optimal open reading frames, resolves overlapping gene loci, and finds additional gene copies where they exist. LiftOn can reliably transfer annotation between genomes representing members of the same species, as we demonstrate on human, mouse, honey bee, rice, and Arabidopsis thaliana. It can further map annotation effectively across species pairs as far apart as mouse and rat or Drosophila melanogaster and D. erecta.

Isogenic hiPSC models of Turner syndrome development reveal shared roles of inactive X and Y in the human cranial neural crest network

Modeling the developmental etiology of viable human aneuploidy can be challenging in rodents due to syntenic boundaries, or primate-specific biology. In humans, monosomy-X (45,X) causes Turner syndrome (TS), altering craniofacial, skeletal, endocrine, and cardiovascular development, which in contrast remain unaffected in 39,X-mice. To learn how human monosomy-X may impact early embryonic development, we turned to human 45,X and isogenic euploid induced pluripotent stem cells (hiPSCs) from male and female mosaic donors. Because neural crest (NC) derived cell types are hypothesized to underpin craniofacial and cardiovascular changes in TS, we performed a highly-powered differential expression study on hiPSC-derived anterior neural crest cells (NCCs). Across three independent isogenic panels, 45,X NCCs show impaired acquisition of PAX7+SOX10+ markers, and disrupted expression of other NCC-specific genes, relative to their isogenic euploid controls. In particular, 45,X NCCs increase cholesterol biosynthesis genes while reducing transcripts that feature 5' terminal oligopyrimidine (TOP) motifs, including those of ribosomal protein and nuclear-encoded mitochondrial genes. Such metabolic pathways are also over-represented in weighted co-expression gene modules that are preserved in monogenic neurocristopathy. Importantly, these gene modules are also significantly enriched in 28% of all TS-associated terms of the human phenotype ontology. Our analysis identifies specific sex-linked genes that are expressed from two copies in euploid males and females alike and qualify as candidate haploinsufficient drivers of TS phenotypes in NC-derived lineages. This study demonstrates that isogenic hiPSC-derived NCC panels representing monosomy-X can serve as a powerful model of early NC development in TS and inform new hypotheses towards its etiology.

Serine synthesis pathway upregulated by E-cadherin is essential for the proliferation and metastasis of breast cancers

The loss of E-cadherin (E-cad), an epithelial cell adhesion molecule, has been implicated in the epithelial-mesenchymal transition (EMT), promoting invasion and migration of cancer cells and, consequently, metastasis. However, recent studies have demonstrated that E-cad supports the survival and proliferation of metastatic cancer cells, suggesting that our understanding of E-cad in metastasis is far from comprehensive. Here, we report that E-cad upregulates the de novo serine synthesis pathway (SSP) in breast cancer cells. The SSP provides metabolic precursors for biosynthesis and resistance to oxidative stress, critically beneficial for E-cad-positive breast cancer cells to achieve faster tumor growth and more metastases. Inhibition of PHGDH, a rate-limiting enzyme in the SSP, significantly and specifically hampered the proliferation of E-cad-positive breast cancer cells and rendered them vulnerable to oxidative stress, inhibiting their metastatic potential. Our findings reveal that E-cad adhesion molecule significantly reprograms cellular metabolism, promoting tumor growth and metastasis of breast cancers.

Identification of 27 allele-specific regulatory variants in Parkinson's disease using a massively parallel reporter assay

Genome wide association studies (GWAS) have identified a number of genomic loci that are associated with Parkinson's disease (PD) risk. However, the majority of these variants lie in non-coding regions, and thus the mechanisms by which they influence disease development, and/or potential subtypes, remain largely elusive. To address this, we used a massively parallel reporter assay (MPRA) to screen the regulatory function of 5254 variants that have a known or putative connection to PD. We identified 138 loci with enhancer activity, of which 27 exhibited allele-specific regulatory activity in HEK293 cells. The identified regulatory variant(s) typically did not match the original tag variant within the PD associated locus, supporting the need for deeper exploration of these loci. The existence of allele specific transcriptional impacts within HEK293 cells, confirms that at least a subset of the PD associated regions mark functional gene regulatory elements. Future functional studies that confirm the putative targets of the empirically verified regulatory variants will be crucial for gaining a greater understanding of how gene regulatory network(s) modulate PD risk.

Temporal and spatial dynamics of Plasmodium falciparum clonal lineages in Guyana

Plasmodium parasites, the causal agents of malaria, are eukaryotic organisms that obligately undergo sexual recombination within mosquitoes. However, in low transmission settings where most mosquitoes become infected with only a single parasite clone, parasites recombine with themselves, and the clonal lineage is propagated rather than broken up by outcrossing. We investigated whether stochastic/neutral factors drive the persistence and abundance of Plasmodium falciparum clonal lineages in Guyana, a country with relatively low malaria transmission, but the only setting in the Americas in which an important artemisinin resistance mutation (pfk13 C580Y) has been observed. To investigate whether this clonality was potentially associated with the persistence and spatial spread of the mutation, we performed whole genome sequencing on 1,727 Plasmodium falciparum samples collected from infected patients across a five-year period (2016-2021). We characterized the relatedness between each pair of monoclonal infections (n=1,409) through estimation of identity by descent (IBD) and also typed each sample for known or candidate drug resistance mutations. A total of 160 clones (mean IBD ≥ 0.90) were circulating in Guyana during the study period, comprising 13 highly related clusters (mean IBD ≥ 0.40). In the five-year study period, we observed a decrease in frequency of a mutation associated with artemisinin partner drug (piperaquine) resistance (pfcrt C350R) and limited co-occurence of pfcrt C350R with duplications of plasmepsin 2/3, an epistatic interaction associated with piperaquine resistance. We additionally report polymorphisms exhibiting evidence of selection for drug resistance or other phenotypes and reported a novel pfk13 mutation (G718S) as well as 61 nonsynonymous substitutions that increased markedly in frequency. However, P. falciparum clonal dynamics in Guyana appear to be largely driven by stochastic factors, in contrast to other geographic regions. The use of multiple artemisinin combination therapies in Guyana may have contributed to the disappearance of the pfk13 C580Y mutation.

Integrative functional genomic analyses identify genetic variants influencing skin pigmentation in Africans

Skin color is highly variable in Africans, yet little is known about the underlying molecular mechanism. Here we applied massively parallel reporter assays to screen 1,157 candidate variants influencing skin pigmentation in Africans and identified 165 single-nucleotide polymorphisms showing differential regulatory activities between alleles. We combine Hi-C, genome editing and melanin assays to identify regulatory elements for MFSD12, HMG20B, OCA2, MITF, LEF1, TRPS1, BLOC1S6 and CYB561A3 that impact melanin levels in vitro and modulate human skin color. We found that independent mutations in an OCA2 enhancer contribute to the evolution of human skin color diversity and detect signals of local adaptation at enhancers of MITF, LEF1 and TRPS1, which may contribute to the light skin color of Khoesan-speaking populations from Southern Africa. Additionally, we identified CYB561A3 as a novel pigmentation regulator that impacts genes involved in oxidative phosphorylation and melanogenesis. These results provide insights into the mechanisms underlying human skin color diversity and adaptive evolution.

A Generalized Higher-order Correlation Analysis Framework for Multi-Omics Network Inference

Multiple -omics (genomics, proteomics, etc.) profiles are commonly generated to gain insight into a disease or physiological system. Constructing multi-omics networks with respect to the trait(s) of interest provides an opportunity to understand relationships between molecular features but integration is challenging due to multiple data sets with high dimensionality. One approach is to use canonical correlation to integrate one or two omics types and a single trait of interest. However, these types of methods may be limited due to (1) not accounting for higher-order correlations existing among features, (2) computational inefficiency when extending to more than two omics data when using a penalty term-based sparsity method, and (3) lack of flexibility for focusing on specific correlations (e.g., omics-to-phenotype correlation versus omics-to-omics correlations). In this work, we have developed a novel multi-omics network analysis pipeline called Sparse Generalized Tensor Canonical Correlation Analysis Network Inference (SGTCCA-Net) that can effectively overcome these limitations. We also introduce an implementation to improve the summarization of networks for downstream analyses. Simulation and real-data experiments demonstrate the effectiveness of our novel method for inferring omics networks and features of interest.

Charting Single Cell Lineage Dynamics and Mutation Networks via Homing CRISPR

Single cell lineage tracing, essential for unraveling cellular dynamics in disease evolution is critical for developing targeted therapies. CRISPR-Cas9, known for inducing permanent and cumulative mutations, is a cornerstone in lineage tracing. The novel homing guide RNA (hgRNA) technology enhances this by enabling dynamic retargeting and facilitating ongoing genetic modifications. Charting these mutations, especially through successive hgRNA edits, poses a significant challenge. Our solution, LINEMAP, is a computational framework designed to trace and map these mutations with precision. LINEMAP meticulously discerns mutation alleles at single-cell resolution and maps their complex interrelationships through a mutation evolution network. By utilizing a Markov Process model, we can predict mutation transition probabilities, revealing potential mutational routes and pathways. Our reconstruction algorithm, anchored in the Markov model's attributes, reconstructs cellular lineage pathways, shedding light on the cell's evolutionary journey to the minutiae of single-cell division. Our findings reveal an intricate network of mutation evolution paired with a predictive Markov model, advancing our capability to reconstruct single-cell lineage via hgRNA. This has substantial implications for advancing our understanding of biological mechanisms and propelling medical research forward.

Signatures of transmission in within-host M. tuberculosis variation

Background

Because M. tuberculosis evolves slowly, transmission clusters often contain multiple individuals with identical consensus genomes, making it difficult to reconstruct transmission chains. Finding additional sources of shared M. tuberculosis variation could help overcome this problem. Previous studies have reported M. tuberculosis diversity within infected individuals; however, whether within-host variation improves transmission inferences remains unclear.

Methods

To evaluate the transmission information present in within-host M. tuberculosis variation, we re-analyzed publicly available sequence data from three household transmission studies, using household membership as a proxy for transmission linkage between donor-recipient pairs.

Findings

We found moderate levels of minority variation present in M. tuberculosis sequence data from cultured isolates that varied significantly across studies (mean: 6, 7, and 170 minority variants above a 1% minor allele frequency threshold, outside of PE/PPE genes). Isolates from household members shared more minority variants than did isolates from unlinked individuals in the three studies (mean 98 shared minority variants vs. 10; 0.8 vs. 0.2, and 0.7 vs. 0.2, respectively). Shared within-host variation was significantly associated with household membership (OR: 1.51 [1.30,1.71], for one standard deviation increase in shared minority variants). Models that included shared within-host variation improved the accuracy of predicting household membership in all three studies as compared to models without within-host variation (AUC: 0.95 versus 0.92, 0.99 versus 0.95, and 0.93 versus 0.91).

Interpretation

Within-host M. tuberculosis variation persists through culture and could enhance the resolution of transmission inferences. The substantial differences in minority variation recovered across studies highlights the need to optimize approaches to recover and incorporate within-host variation into automated phylogenetic and transmission inference.

Funding

NIAID: 5K01AI173385.

Systematic investigation of allelic regulatory activity of schizophrenia-associated common variants

Genome-wide association studies (GWASs) have successfully identified 145 genomic regions that contribute to schizophrenia risk, but linkage disequilibrium makes it challenging to discern causal variants. We performed a massively parallel reporter assay (MPRA) on 5,173 fine-mapped schizophrenia GWAS variants in primary human neural progenitors and identified 439 variants with allelic regulatory effects (MPRA-positive variants). Transcription factor binding had modest predictive power, while fine-map posterior probability, enhancer overlap, and evolutionary conservation failed to predict MPRA-positive variants. Furthermore, 64% of MPRA-positive variants did not exhibit expressive quantitative trait loci signature, suggesting that MPRA could identify yet unexplored variants with regulatory potentials. To predict the combinatorial effect of MPRA-positive variants on gene regulation, we propose an accessibility-by-contact model that combines MPRA-measured allelic activity with neuronal chromatin architecture.

High-throughput functional dissection of noncoding SNPs with biased allelic enhancer activity for insulin resistance-relevant phenotypes

Most of the single-nucleotide polymorphisms (SNPs) associated with insulin resistance (IR)-relevant phenotypes by genome-wide association studies (GWASs) are located in noncoding regions, complicating their functional interpretation. Here, we utilized an adapted STARR-seq to evaluate the regulatory activities of 5,987 noncoding SNPs associated with IR-relevant phenotypes. We identified 876 SNPs with biased allelic enhancer activity effects (baaSNPs) across 133 loci in three IR-relevant cell lines (HepG2, preadipocyte, and A673), which showed pervasive cell specificity and significant enrichment for cell-specific open chromatin regions or enhancer-indicative markers (H3K4me1, H3K27ac). Further functional characterization suggested several transcription factors (TFs) with preferential allelic binding to baaSNPs. We also incorporated multi-omics data to prioritize 102 candidate regulatory target genes for baaSNPs and revealed prevalent long-range regulatory effects and cell-specific IR-relevant biological functional enrichment on them. Specifically, we experimentally verified the distal regulatory mechanism at IRS1 locus, in which rs952227-A reinforces IRS1 expression by long-range chromatin interaction and preferential binding to the transcription factor HOXC6 to augment the enhancer activity. Finally, based on our STARR-seq screening data, we predicted the enhancer activity of 227,343 noncoding SNPs associated with IR-relevant phenotypes (fasting insulin adjusted for BMI, HDL cholesterol, and triglycerides) from the largest available GWAS summary statistics. We further provided an open resource (http://www.bigc.online/fnSNP-IR) for better understanding genetic regulatory mechanisms of IR-relevant phenotypes.

HATCHet2: clone- and haplotype-specific copy number inference from bulk tumor sequencing data

Multi-region DNA sequencing of primary tumors and metastases from individual patients helps identify somatic aberrations driving cancer development. However, most methods to infer copy-number aberrations (CNAs) analyze individual samples. We introduce HATCHet2 to identify haplotype- and clone-specific CNAs simultaneously from multiple bulk samples. HATCHet2 introduces a novel statistic, the mirrored haplotype B-allele frequency (mhBAF), to identify mirrored-subclonal CNAs having different numbers of copies of parental haplotypes in different tumor clones. HATCHet2 also has high accuracy in identifying focal CNAs and extends the earlier HATCHet method in several directions. We demonstrate HATCHet2's improved accuracy using simulations and a single-cell sequencing dataset. HATCHet2 analysis of 50 prostate cancer samples from 10 patients reveals previously-unreported mirrored-subclonal CNAs affecting cancer genes.

Rare variants in long non-coding RNAs are associated with blood lipid levels in the TOPMed Whole Genome Sequencing Study

Long non-coding RNAs (lncRNAs) are known to perform important regulatory functions. Large-scale whole genome sequencing (WGS) studies and new statistical methods for variant set tests now provide an opportunity to assess the associations between rare variants in lncRNA genes and complex traits across the genome. In this study, we used high-coverage WGS from 66,329 participants of diverse ancestries with blood lipid levels (LDL-C, HDL-C, TC, and TG) in the National Heart, Lung, and Blood Institute (NHLBI) Trans-Omics for Precision Medicine (TOPMed) program to investigate the role of lncRNAs in lipid variability. We aggregated rare variants for 165,375 lncRNA genes based on their genomic locations and conducted rare variant aggregate association tests using the STAAR (variant-Set Test for Association using Annotation infoRmation) framework. We performed STAAR conditional analysis adjusting for common variants in known lipid GWAS loci and rare coding variants in nearby protein coding genes. Our analyses revealed 83 rare lncRNA variant sets significantly associated with blood lipid levels, all of which were located in known lipid GWAS loci (in a ±500 kb window of a Global Lipids Genetics Consortium index variant). Notably, 61 out of 83 signals (73%) were conditionally independent of common regulatory variations and rare protein coding variations at the same loci. We replicated 34 out of 61 (56%) conditionally independent associations using the independent UK Biobank WGS data. Our results expand the genetic architecture of blood lipids to rare variants in lncRNA, implicating new therapeutic opportunities.

mutscan-a flexible R package for efficient end-to-end analysis of multiplexed assays of variant effect data

Multiplexed assays of variant effect (MAVE) experimentally measure the effect of large numbers of sequence variants by selective enrichment of sequences with desirable properties followed by quantification by sequencing. mutscan is an R package for flexible analysis of such experiments, covering the entire workflow from raw reads up to statistical analysis and visualization. The core components are implemented in C++ for efficiency. Various experimental designs are supported, including single or paired reads with optional unique molecular identifiers. To find variants with changed relative abundance, mutscan employs established statistical models provided in the edgeR and limma packages. mutscan is available from https://github.com/fmicompbio/mutscan .

GALBA: Genome Annotation with Miniprot and AUGUSTUS

The Earth Biogenome Project has rapidly increased the number of available eukaryotic genomes, but most released genomes continue to lack annotation of protein-coding genes. In addition, no transcriptome data is available for some genomes. Various gene annotation tools have been developed but each has its limitations. Here, we introduce GALBA, a fully automated pipeline that utilizes miniprot, a rapid protein- to-genome aligner, in combination with AUGUSTUS to predict genes with high accuracy. Accuracy results indicate that GALBA is particularly strong in the annotation of large vertebrate genomes. We also present use cases in insects, vertebrates, and a previously unannotated land plant. GALBA is fully open source and available as a docker image for easy execution with Singularity in high-performance computing environments. Our pipeline addresses the critical need for accurate gene annotation in newly sequenced genomes, and we believe that GALBA will greatly facilitate genome annotation for diverse organisms.

Arsenic is a potent co-mutagen of ultraviolet light

Environmental co-exposures are widespread and are major contributors to carcinogenic mechanisms. Two well-established environmental agents causing skin cancer are ultraviolet radiation (UVR) and arsenic. Arsenic is a known co-carcinogen that enhances UVR's carcinogenicity. However, the mechanisms of arsenic co-carcinogenesis are not well understood. In this study, we utilized primary human keratinocytes and a hairless mouse model to investigate the carcinogenic and mutagenic properties of co-exposure to arsenic and UVR. In vitro and in vivo exposures revealed that, on its own, arsenic is neither mutagenic nor carcinogenic. However, in combination with UVR, arsenic exposure has a synergistic effect leading to an accelerated mouse skin carcinogenesis as well as to more than 2-fold enrichment of UVR mutational burden. Notably, mutational signature ID13, previously found only in UVR-associated human skin cancers, was observed exclusively in mouse skin tumors and cell lines jointly exposed to arsenic and UVR. This signature was not observed in any model system exposed purely to arsenic or purely to UVR, making ID13 the first co-exposure signature to be reported using controlled experimental conditions. Analysis of existing genomics data from basal cell carcinomas and melanomas revealed that only a subset of human skin cancers harbor ID13 and, consistent with our experimental observations, these cancers exhibited an elevated UVR mutagenesis. Our results provide the first report of a unique mutational signature caused by a co-exposure to two environmental carcinogens and the first comprehensive evidence that arsenic is a potent co-mutagen and co-carcinogen of UVR. Importantly, our findings suggest that a large proportion of human skin cancers are not formed purely due to UVR exposure but rather due to a co-exposure of UVR and other co-mutagens such as arsenic.

Leveraging massively parallel reporter assays for evolutionary questions

A long-standing goal of evolutionary biology is to decode how gene regulation contributes to organismal diversity. Doing so is challenging because it is hard to predict function from non-coding sequence and to perform molecular research with non-model taxa. Massively parallel reporter assays (MPRAs) enable the testing of thousands to millions of sequences for regulatory activity simultaneously. Here, we discuss the execution, advantages, and limitations of MPRAs, with a focus on evolutionary questions. We propose solutions for extending MPRAs to rare taxa and those with limited genomic resources, and we underscore MPRA's broad potential for driving genome-scale, functional studies across organisms.

Functional investigation of inherited noncoding genetic variation impacting the pharmacogenomics of childhood acute lymphoblastic leukemia treatment

Although acute lymphoblastic leukemia (ALL) is the most common childhood cancer, there is limited understanding of the contribution of inherited genetic variation on inter-individual differences in chemotherapy response. Defining genetic factors impacting therapy failure can help better predict response and identify drug resistance mechanisms. We therefore mapped inherited noncoding variants associated with chemotherapeutic drug resistance and/or treatment outcome to ALL cis-regulatory elements and investigated their gene regulatory potential and genomic connectivity using massively parallel reporter assays and promoter capture Hi-C, respectively. We identified 53 variants with reproducible allele-specific effects on transcription and high-confidence gene targets. Subsequent functional interrogation of the top variant (rs1247117) determined that it disrupted a PU.1 consensus motif and PU.1 binding affinity. Importantly, deletion of the genomic interval containing rs1247117 sensitized ALL cells to vincristine. Together, these data demonstrate that noncoding regulatory variation associated with diverse pharmacological traits harbor significant effects on allele-specific transcriptional activity and impact sensitivity to chemotherapeutic agents in ALL.

Modeling common and rare genetic risk factors of neuropsychiatric disorders in human induced pluripotent stem cells

Recent genome-wide association studies (GWAS) and whole-exome sequencing of neuropsychiatric disorders, especially schizophrenia, have identified a plethora of common and rare disease risk variants/genes. Translating the mounting human genetic discoveries into novel disease biology and more tailored clinical treatments is tied to our ability to causally connect genetic risk variants to molecular and cellular phenotypes. When combined with the Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR)/CRISPR-associated (Cas) nuclease-mediated genome editing system, human induced pluripotent stem cell (hiPSC)-derived neural cultures (both 2D and 3D organoids) provide a promising tractable cellular model for bridging the gap between genetic findings and disease biology. In this review, we first conceptualize the advances in understanding the disease polygenicity and convergence from the past decade of iPSC modeling of different types of genetic risk factors of neuropsychiatric disorders. We then discuss the major cell types and cellular phenotypes that are most relevant to neuropsychiatric disorders in iPSC modeling. Finally, we critically review the limitations of iPSC modeling of neuropsychiatric disorders and outline the need for implementing and developing novel methods to scale up the number of iPSC lines and disease risk variants in a systematic manner. Sufficiently scaled-up iPSC modeling and a better functional interpretation of genetic risk variants, in combination with cutting-edge CRISPR/Cas9 gene editing and single-cell multi-omics methods, will enable the field to identify the specific and convergent molecular and cellular phenotypes in precision for neuropsychiatric disorders.

Plasmid and Sequencing Library Preparation for CRISPRi Barcoded Expression Reporter Sequencing (CiBER-seq) in Saccharomyces cerevisiae

Genetic networks regulate nearly all biological processes, including cellular differentiation, homeostasis, and immune responses. Determining the precise role of each gene within a regulatory network can explain its overall, integrated function, and pinpoint mechanisms underlying misregulation in disease states. Transcriptional reporter assays are a useful tool for dissecting these genetic networks, because they link a molecular process to a measurable readout, such as the expression of a fluorescent protein. Here, we introduce a new technique that uses expressed RNA barcodes as reporters, to measure transcriptional changes induced by CRISPRi-mediated genetic perturbation across a diverse, genome-wide library of guide RNAs. We describe an exemplary reporter based on the promoter that drives His4 expression in these guidelines, which can be used as a framework to interrogate other expression phenotypes. In this workflow, a library of plasmids is assembled, encoding a CRISPRi guide RNA (gRNA) along with one or more transcriptional reporters that drive expression of guide-specific nucleotide barcode sequences. For example, when interrogating regulation of the budding yeast HIS4 promoter normalized against a control housekeeping promoter that drives Pgk1 expression, this plasmid library contains a gRNA expression cassette, a HIS4 reporter driving expression of one gRNA-specific nucleotide barcode, and a PGK1 reporter driving expression of a second, gRNA-specific barcode. Long-read sequencing is used to determine which gRNA is associated with these nucleotide barcodes. The plasmid library is then transformed into yeast cells, where each cell receives one plasmid, and experiences a genetic perturbation driven by the guide on that plasmid. The expressed RNA barcodes are extracted in bulk and quantified using high-throughput sequencing, thereby measuring the effect of their corresponding gRNA on barcoded reporter expression. In the case of the HIS4 reporter described above, guides disrupting translation elongation will increase expression of the associated HIS4 barcode specifically, without changing expression of the PGK1 control barcode. It is further possible to quantify plasmid abundance by DNA sequencing, as an additional approach to normalize for differences in plasmid abundance within the population of cells. This protocol outlines the steps to prepare barcode reporter CRISPRi plasmid libraries, link guides to barcodes with long-read sequencing, and measure expression changes through barcode RNA and DNA sequencing. This method is ideal for probing transcriptional or post-transcriptional regulation, as it measures the effects of a genetic perturbation by directly quantifying reporter RNA abundance, rather than relying on indirect growth or fluorescence readouts. Graphic abstract.

Functional noncoding SNPs in human endothelial cells fine-map vascular trait associations

Functional consequences of genetic variation in the noncoding human genome are difficult to ascertain despite demonstrated associations to common, complex disease traits. To elucidate properties of functional noncoding SNPs with effects in human endothelial cells (ECs), we utilized our previous molecular quantitative trait locus (molQTL) analysis for transcription factor binding, chromatin accessibility, and H3K27 acetylation to nominate a set of likely functional noncoding SNPs. Together with information from genome-wide association studies (GWASs) for vascular disease traits, we tested the ability of 34,344 variants to perturb enhancer function in ECs using the highly multiplexed STARR-seq assay. Of these, 5711 variants validated, whose enriched attributes included: (1) mutations to TF binding motifs for ETS or AP-1 that are regulators of the EC state; (2) location in accessible and H3K27ac-marked EC chromatin; and (3) molQTL associations whereby alleles associate with differences in chromatin accessibility and TF binding across genetically diverse ECs. Next, using pro-inflammatory IL1B as an activator of cell state, we observed robust evidence (>50%) of context-specific SNP effects, underscoring the prevalence of noncoding gene-by-environment (GxE) effects. Lastly, using these cumulative data, we fine-mapped vascular disease loci and highlighted evidence suggesting mechanisms by which noncoding SNPs at two loci affect risk for pulse pressure/large artery stroke and abdominal aortic aneurysm through respective effects on transcriptional regulation of POU4F1 and LDAH. Together, we highlight the attributes and context dependence of functional noncoding SNPs and provide new mechanisms underlying vascular disease risk.

Single-Cell Epigenomics and Functional Fine-Mapping of Atherosclerosis GWAS Loci

Supplemental Digital Content is available in the text.

Genome-wide association studies have identified hundreds of loci associated with coronary artery disease (CAD). Many of these loci are enriched in cisregulatory elements but not linked to cardiometabolic risk factors nor to candidate causal genes, complicating their functional interpretation.

MPRAdecoder: Processing of the Raw MPRA Data With a priori Unknown Sequences of the Region of Interest and Associated Barcodes

Massively parallel reporter assays (MPRAs) enable high-throughput functional evaluation of numerous DNA regulatory elements and/or their mutant variants. The assays are based on the construction of reporter plasmid libraries containing two variable parts, a region of interest (ROI) and a barcode (BC), located outside and within the transcription unit, respectively. Importantly, each plasmid molecule in a such a highly diverse library is characterized by a unique BC-ROI association. The reporter constructs are delivered to target cells and expression of BCs at the transcript level is assayed by RT-PCR followed by next-generation sequencing (NGS). The obtained values are normalized to the abundance of BCs in the plasmid DNA sample. Altogether, this allows evaluating the regulatory potential of the associated ROI sequences. However, depending on the MPRA library construction design, the BC and ROI sequences as well as their associations can be a priori unknown. In such a case, the BC and ROI sequences, their possible mutant variants, and unambiguous BC-ROI associations have to be identified, whereas all uncertain cases have to be excluded from the analysis. Besides the preparation of additional "mapping" samples for NGS, this also requires specific bioinformatics tools. Here, we present a pipeline for processing raw MPRA data obtained by NGS for reporter construct libraries with a priori unknown sequences of BCs and ROIs. The pipeline robustly identifies unambiguous (so-called genuine) BCs and ROIs associated with them, calculates the normalized expression level for each BC and the averaged values for each ROI, and provides a graphical visualization of the processed data.

A novel family of secreted insect proteins linked to plant gall development

In an elaborate form of inter-species exploitation, many insects hijack plant development to induce novel plant organs called galls that provide the insect with a source of nutrition and a temporary home. Galls result from dramatic reprogramming of plant cell biology driven by insect molecules, but the roles of specific insect molecules in gall development have not yet been determined. Here we study the aphidHormaphis cornu, which makes distinctive “cone” galls on leaves of witch hazelHamamelis virginiana. We found that derived genetic variants in the aphid genedeterminant of gall color(dgc) are associated with strong downregulation ofdgctranscription in aphid salivary glands, upregulation in galls of seven genes involved in anthocyanin synthesis, and deposition of two red anthocyanins in galls. We hypothesize that aphids inject DGC protein into galls, and that this results in differential expression of a small number of plant genes.Dgcis a member of a large, diverse family of novel predicted secreted proteins characterized by a pair of widely spaced cysteine-tyrosine-cysteine (CYC) residues, which we named BICYCLE proteins.Bicyclegenes are most strongly expressed in the salivary glands specifically of galling aphid generations, suggesting that they may regulate many aspects of gall development.Bicyclegenes have experienced unusually frequent diversifying selection, consistent with their potential role controlling gall development in a molecular arms race between aphids and their host plants.

One Sentence Summary

Aphidbicyclegenes, which encode diverse secreted proteins, contribute to plant gall development.

Massively parallel discovery of human-specific substitutions that alter enhancer activity

Genetic changes that altered the function of gene regulatory elements have been implicated in the evolution of human traits such as the expansion of the cerebral cortex. However, identifying the particular changes that modified regulatory activity during human evolution remain challenging. Here we used massively parallel enhancer assays in neural stem cells to quantify the functional impact of >32,000 human-specific substitutions in >4,300 human accelerated regions (HARs) and human gain enhancers (HGEs), which include enhancers with novel activities in humans. We found that >30% of active HARs and HGEs exhibited differential activity between human and chimpanzee. We isolated the effects of human-specific substitutions from background genetic variation to identify the effects of genetic changes most relevant to human evolution. We found that substitutions interacted in both additive and nonadditive ways to modify enhancer function. Substitutions within HARs, which are highly constrained compared to HGEs, showed smaller effects on enhancer activity, suggesting that the impact of human-specific substitutions is buffered in enhancers with constrained ancestral functions. Our findings yield insight into how human-specific genetic changes altered enhancer function and provide a rich set of candidates for studies of regulatory evolution in humans.

Massively Parallel Reporter Assays: Defining Functional Psychiatric Genetic Variants Across Biological Contexts

Neuropsychiatric phenotypes have long been known to be influenced by heritable risk factors, directly confirmed by the past decade of genetic studies that have revealed specific genetic variants enriched in disease cohorts. However, the initial hope that a small set of genes would be responsible for a given disorder proved false. The more complex reality is that a given disorder may be influenced by myriad small-effect noncoding variants and/or by rare but severe coding variants, many de novo. Noncoding genomic sequences-for which molecular functions cannot usually be inferred-harbor a large portion of these variants, creating a substantial barrier to understanding higher-order molecular and biological systems of disease. Fortunately, novel genetic technologies-scalable oligonucleotide synthesis, RNA sequencing, and CRISPR (clustered regularly interspaced short palindromic repeats)-have opened novel avenues to experimentally identify biologically significant variants en masse. Massively parallel reporter assays (MPRAs) are an especially versatile technique resulting from such innovations. MPRAs are powerful molecular genetics tools that can be used to screen thousands of untranscribed or untranslated sequences and their variants for functional effects in a single experiment. This approach, though underutilized in psychiatric genetics, has several useful features for the field. We review methods for assaying putatively functional genetic variants and regions, emphasizing MPRAs and the opportunities they hold for dissection of psychiatric polygenicity. We discuss literature applying functional assays in neurogenetics, highlighting strengths, caveats, and design considerations-especially regarding disease-relevant variables (cell type, neurodevelopment, and sex), and we ultimately propose applications of MPRA to both computational and experimental neurogenetics of polygenic disease risk.

Gamete binning: chromosome-level and haplotype-resolved genome assembly enabled by high-throughput single-cell sequencing of gamete genomes

Generating chromosome-level, haplotype-resolved assemblies of heterozygous genomes remains challenging. To address this, we developed gamete binning, a method based on single-cell sequencing of haploid gametes enabling separation of the whole-genome sequencing reads into haplotype-specific reads sets. After assembling the reads of each haplotype, the contigs are scaffolded to chromosome level using a genetic map derived from the gametes. We assemble the two genomes of a diploid apricot tree based on whole-genome sequencing of 445 individual pollen grains. The two haplotype assemblies (N50: 25.5 and 25.8 Mb) feature a haplotyping precision of greater than 99% and are accurately scaffolded to chromosome-level.

CiBER-seq dissects genetic networks by quantitative CRISPRi profiling of expression phenotypes

To realize the promise of CRISPR-Cas9-based genetics, approaches are needed to quantify a specific, molecular phenotype across genome-wide libraries of genetic perturbations. We addressed this challenge by profiling transcriptional, translational, and posttranslational reporters using CRISPR interference (CRISPRi) with barcoded expression reporter sequencing (CiBER-seq). Our barcoding approach allowed us to connect an entire library of guides to their individual phenotypic consequences using pooled sequencing. CiBER-seq profiling fully recapitulated the integrated stress response (ISR) pathway in yeast. Genetic perturbations causing uncharged transfer RNA (tRNA) accumulation activated ISR reporter transcription. Notably, tRNA insufficiency also activated the reporter, independent of the uncharged tRNA sensor. By uncovering alternate triggers for ISR activation, we illustrate how precise, comprehensive CiBER-seq profiling provides a powerful and broadly applicable tool for dissecting genetic networks.

Systematic identification of cis-regulatory variants that cause gene expression differences in a yeast cross

Sequence variation in regulatory DNA alters gene expression and shapes genetically complex traits. However, the identification of individual, causal regulatory variants is challenging. Here, we used a massively parallel reporter assay to measure the cis-regulatory consequences of 5832 natural DNA variants in the promoters of 2503 genes in the yeast Saccharomyces cerevisiae. We identified 451 causal variants, which underlie genetic loci known to affect gene expression. Several promoters harbored multiple causal variants. In five promoters, pairs of variants showed non-additive, epistatic interactions. Causal variants were enriched at conserved nucleotides, tended to have low derived allele frequency, and were depleted from promoters of essential genes, which is consistent with the action of negative selection. Causal variants were also enriched for alterations in transcription factor binding sites. Models integrating these features provided modest, but statistically significant, ability to predict causal variants. This work revealed a complex molecular basis for cis-acting regulatory variation.

Common genetic risk variants identified in the SPARK cohort support DDHD2 as a candidate risk gene for autism

Autism spectrum disorder (ASD) is a highly heritable neurodevelopmental disorder. Large genetically informative cohorts of individuals with ASD have led to the identification of a limited number of common genome-wide significant (GWS) risk loci to date. However, many more common genetic variants are expected to contribute to ASD risk given the high heritability. Here, we performed a genome-wide association study (GWAS) on 6222 case-pseudocontrol pairs from the Simons Foundation Powering Autism Research for Knowledge (SPARK) dataset to identify additional common genetic risk factors and molecular mechanisms underlying risk for ASD. We identified one novel GWS locus from the SPARK GWAS and four significant loci, including an additional novel locus from meta-analysis with a previous GWAS. We replicated the previous observation of significant enrichment of ASD heritability within regulatory regions of the developing cortex, indicating that disruption of gene regulation during neurodevelopment is critical for ASD risk. We further employed a massively parallel reporter assay (MPRA) and identified a putative causal variant at the novel locus from SPARK GWAS with strong impacts on gene regulation (rs7001340). Expression quantitative trait loci data demonstrated an association between the risk allele and decreased expression of DDHD2 (DDHD domain containing 2) in both adult and prenatal brains. In conclusion, by integrating genetic association data with multi-omic gene regulatory annotations and experimental validation, we fine-mapped a causal risk variant and demonstrated that DDHD2 is a novel gene associated with ASD risk.

Evaluation of Davis et al.: Exploring Sequence of Determinants of Transcriptional Regulation-The Case of c-AMP Response Element

One snapshot of the peer review process for "Dissection of c-AMP Response Element Architecture by Using Genomic and Episomal Massively Parallel Reporter Assays" (Davis et al., 2020).

Statistical considerations for the analysis of massively parallel reporter assays data

Noncoding DNA contains gene regulatory elements that alter gene expression, and the function of these elements can be modified by genetic variation. Massively parallel reporter assays (MPRA) enable high-throughput identification and characterization of functional genetic variants, but the statistical methods to identify allelic effects in MPRA data have not been fully developed. In this study, we demonstrate how the baseline allelic imbalance in MPRA libraries can produce biased results, and we propose a novel, nonparametric, adaptive testing method that is robust to this bias. We compare the performance of this method with other commonly used methods, and we demonstrate that our novel adaptive method controls Type I error in a wide range of scenarios while maintaining excellent power. We have implemented these tests along with routines for simulating MPRA data in the Analysis Toolset for MPRA (@MPRA), an R package for the design and analyses of MPRA experiments. It is publicly available at http://github.com/redaq/atMPRA.

Bayesian modelling of high-throughput sequencing assays with malacoda

NGS studies have uncovered an ever-growing catalog of human variation while leaving an enormous gap between observed variation and experimental characterization of variant function. High-throughput screens powered by NGS have greatly increased the rate of variant functionalization, but the development of comprehensive statistical methods to analyze screen data has lagged. In the massively parallel reporter assay (MPRA), short barcodes are counted by sequencing DNA libraries transfected into cells and the cell's output RNA in order to simultaneously measure the shifts in transcription induced by thousands of genetic variants. These counts present many statistical challenges, including overdispersion, depth dependence, and uncertain DNA concentrations. So far, the statistical methods used have been rudimentary, employing transformations on count level data and disregarding experimental and technical structure while failing to quantify uncertainty in the statistical model. We have developed an extensive framework for the analysis of NGS functionalization screens available as an R package called malacoda (available from github.com/andrewGhazi/malacoda). Our software implements a probabilistic, fully Bayesian model of screen data. The model uses the negative binomial distribution with gamma priors to model sequencing counts while accounting for effects from input library preparation and sequencing depth. The method leverages the high-throughput nature of the assay to estimate the priors empirically. External annotations such as ENCODE data or DeepSea predictions can also be incorporated to obtain more informative priors-a transformative capability for data integration. The package also includes quality control and utility functions, including automated barcode counting and visualization methods. To validate our method, we analyzed several datasets using malacoda and alternative MPRA analysis methods. These data include experiments from the literature, simulated assays, and primary MPRA data. We also used luciferase assays to experimentally validate several hits from our primary data, as well as variants for which the various methods disagree and variants detectable only with the aid of external annotations.

A screen of 1,049 schizophrenia and 30 Alzheimer's-associated variants for regulatory potential

Recent genome-wide association studies (GWAS) identified numerous schizophrenia (SZ) and Alzheimer's disease (AD) associated loci, most outside protein-coding regions and hypothesized to affect gene transcription. We used a massively parallel reporter assay to screen, 1,049 SZ and 30 AD variants in 64 and nine loci, respectively for allele differences in driving reporter gene expression. A library of synthetic oligonucleotides assaying each allele five times was transfected into K562 chronic myelogenous leukemia lymphoblasts and SK-SY5Y human neuroblastoma cells. One hundred forty eight variants showed allelic differences in K562 and 53 in SK-SY5Y cells, on average 2.6 variants per locus. Nine showed significant differences in both lines, a modest overlap reflecting different regulatory landscapes of these lines that also differ significantly in chromatin marks. Eight of nine were in the same direction. We observe no preference for risk alleles to increase or decrease expression. We find a positive correlation between the number of SNPs in linkage disequilibrium and the proportion of functional SNPs supporting combinatorial effects that may lead to haplotype selection. Our results prioritize future functional follow up of disease associated SNPs to determine the driver GWAS variant(s), at each locus and enhance our understanding of gene regulation dynamics.

MPRAnalyze: statistical framework for massively parallel reporter assays

Massively parallel reporter assays (MPRAs) can measure the regulatory function of thousands of DNA sequences in a single experiment. Despite growing popularity, MPRA studies are limited by a lack of a unified framework for analyzing the resulting data. Here we present MPRAnalyze: a statistical framework for analyzing MPRA count data. Our model leverages the unique structure of MPRA data to quantify the function of regulatory sequences, compare sequences' activity across different conditions, and provide necessary flexibility in an evolving field. We demonstrate the accuracy and applicability of MPRAnalyze on simulated and published data and compare it with existing methods.