Characterizing the genetic basis of transcriptome diversity through RNA-sequencing of 922 individuals

Genetic variant predictors of gene expression provide new insight into risk of colorectal cancer

Genome-wide association studies have reported 56 independently associated colorectal cancer (CRC) risk variants, most of which are non-coding and believed to exert their effects by modulating gene expression. The computational method PrediXcan uses cis-regulatory variant predictors to impute expression and perform gene-level association tests in GWAS without directly measured transcriptomes. In this study, we used reference datasets from colon (n = 169) and whole blood (n = 922) transcriptomes to test CRC association with genetically determined expression levels in a genome-wide analysis of 12,186 cases and 14,718 controls. Three novel associations were discovered from colon transverse models at FDR ≤ 0.2 and further evaluated in an independent replication including 32,825 cases and 39,933 controls. After adjusting for multiple comparisons, we found statistically significant associations using colon transcriptome models with TRIM4 (discovery P = 2.2 × 10- 4, replication P = 0.01), and PYGL (discovery P = 2.3 × 10- 4, replication P = 6.7 × 10- 4). Interestingly, both genes encode proteins that influence redox homeostasis and are related to cellular metabolic reprogramming in tumors, implicating a novel CRC pathway linked to cell growth and proliferation. Defining CRC risk regions as one megabase up- and downstream of one of the 56 independent risk variants, we defined 44 non-overlapping CRC-risk regions. Among these risk regions, we identified genes associated with CRC (P < 0.05) in 34/44 CRC-risk regions. Importantly, CRC association was found for two genes in the previously reported 2q25 locus, CXCR1 and CXCR2, which are potential cancer therapeutic targets. These findings provide strong candidate genes to prioritize for subsequent laboratory follow-up of GWAS loci. This study is the first to implement PrediXcan in a large colorectal cancer study and findings highlight the utility of integrating transcriptome data in GWAS for discovery of, and biological insight into, risk loci.

Probabilistic fine-mapping of transcriptome-wide association studies

Transcriptome-wide association studies using predicted expression have identified thousands of genes whose locally regulated expression is associated with complex traits and diseases. In this work, we show that linkage disequilibrium induces significant gene-trait associations at non-causal genes as a function of the expression quantitative trait loci weights used in expression prediction. We introduce a probabilistic framework that models correlation among transcriptome-wide association study signals to assign a probability for every gene in the risk region to explain the observed association signal. Importantly, our approach remains accurate when expression data for causal genes are not available in the causal tissue by leveraging expression prediction from other tissues. Our approach yields credible sets of genes containing the causal gene at a nominal confidence level (for example, 90%) that can be used to prioritize genes for functional assays. We illustrate our approach by using an integrative analysis of lipid traits, where our approach prioritizes genes with strong evidence for causality.

Using genetic drug-target networks to develop new drug hypotheses for major depressive disorder

The major depressive disorder (MDD) working group of the Psychiatric Genomics Consortium (PGC) has published a genome-wide association study (GWAS) for MDD in 130,664 cases, identifying 44 risk variants. We used these results to investigate potential drug targets and repurposing opportunities. We built easily interpretable bipartite drug-target networks integrating interactions between drugs and their targets, genome-wide association statistics, and genetically predicted expression levels in different tissues, using the online tool Drug Targetor ( drugtargetor.com ). We also investigated drug-target relationships that could be impacting MDD. MAGMA was used to perform pathway analyses and S-PrediXcan to investigate the directionality of tissue-specific expression levels in patients vs. controls. Outside the major histocompatibility complex (MHC) region, 153 protein-coding genes are significantly associated with MDD in MAGMA after multiple testing correction; among these, five are predicted to be down or upregulated in brain regions and 24 are known druggable genes. Several drug classes were significantly enriched, including monoamine reuptake inhibitors, sex hormones, antipsychotics, and antihistamines, indicating an effect on MDD and potential repurposing opportunities. These findings not only require validation in model systems and clinical examination, but also show that GWAS may become a rich source of new therapeutic hypotheses for MDD and other psychiatric disorders that need new-and better-treatment options.

GRIK5 Genetically Regulated Expression Associated with Eye and Vascular Phenomes: Discovery through Iteration among Biobanks, Electronic Health Records, and Zebrafish

Although the use of model systems for studying the mechanism of mutations that have a large effect is common, we highlight here the ways that zebrafish-model-system studies of a gene, GRIK5, that contributes to the polygenic liability to develop eye diseases have helped to illuminate a mechanism that implicates vascular biology in eye disease. A gene-expression prediction derived from a reference transcriptome panel applied to BioVU, a large electronic health record (EHR)-linked biobank at Vanderbilt University Medical Center, implicated reduced GRIK5 expression in diverse eye diseases. We tested the function of GRIK5 by depletion of its ortholog in zebrafish, and we observed reduced blood vessel numbers and integrity in the eye and increased vascular permeability. Analyses of EHRs in >2.6 million Vanderbilt subjects revealed significant comorbidity of eye and vascular diseases (relative risks 2-15); this comorbidity was confirmed in 150 million individuals from a large insurance claims dataset. Subsequent studies in >60,000 genotyped BioVU participants confirmed the association of reduced genetically predicted expression of GRIK5 with comorbid vascular and eye diseases. Our studies pioneer an approach that allows a rapid iteration of the discovery of gene-phenotype relationships to the primary genetic mechanism contributing to the pathophysiology of human disease. Our findings also add dimension to the understanding of the biology driven by glutamate receptors such as GRIK5 (also referred to as GLUK5 in protein form) and to mechanisms contributing to human eye diseases.

Evolution of DNA Methylation in Papio Baboons

Changes in gene regulation have long been thought to play an important role in primate evolution. However, although a number of studies have compared genome-wide gene expression patterns across primate species, fewer have investigated the gene regulatory mechanisms that underlie such patterns, or the relative contribution of drift versus selection. Here, we profiled genome-scale DNA methylation levels in blood samples from five of the six extant species of the baboon genus Papio (4-14 individuals per species). This radiation presents the opportunity to investigate DNA methylation divergence at both shallow and deeper timescales (0.380-1.4 My). In contrast to studies in human populations, but similar to studies in great apes, DNA methylation profiles clearly mirror genetic and geographic structure. Divergence in DNA methylation proceeds fastest in unannotated regions of the genome and slowest in regions of the genome that are likely more constrained at the sequence level (e.g., gene exons). Both heuristic approaches and Ornstein-Uhlenbeck models suggest that DNA methylation levels at a small set of sites have been affected by positive selection, and that this class is enriched in functionally relevant contexts, including promoters, enhancers, and CpG islands. Our results thus indicate that the rate and distribution of DNA methylation changes across the genome largely mirror genetic structure. However, at some CpG sites, DNA methylation levels themselves may have been a target of positive selection, pointing to loci that could be important in connecting sequence variation to fitness-related traits.

A Multivariate Genome-Wide Association Study of Wing Shape in Drosophila melanogaster

Due to the complexity of genotype-phenotype relationships, simultaneous analyses of genomic associations with multiple traits will be more powerful and informative than a series of univariate analyses. However, in most cases, studies of genotype-phenotype relationships have been analyzed only one trait at a time. Here, we report the results of a fully integrated multivariate genome-wide association analysis of the shape of the Drosophila melanogaster wing in the Drosophila Genetic Reference Panel. Genotypic effects on wing shape were highly correlated between two different laboratories. We found 2396 significant SNPs using a 5% false discovery rate cutoff in the multivariate analyses, but just four significant SNPs in univariate analyses of scores on the first 20 principal component axes. One quarter of these initially significant SNPs retain their effects in regularized models that take into account population structure and linkage disequilibrium. A key advantage of multivariate analysis is that the direction of the estimated phenotypic effect is much more informative than a univariate one. We exploit this fact to show that the effects of knockdowns of genes implicated in the initial screen were on average more similar than expected under a null model. A subset of SNP effects were replicable in an unrelated panel of inbred lines. Association studies that take a phenomic approach, considering many traits simultaneously, are an important complement to the power of genomics.

Nuclear genetic regulation of the human mitochondrial transcriptome

Mitochondria play important roles in cellular processes and disease, yet little is known about how the transcriptional regime of the mitochondrial genome varies across individuals and tissues. By analyzing >11,000 RNA-sequencing libraries across 36 tissue/cell types, we find considerable variation in mitochondrial-encoded gene expression along the mitochondrial transcriptome, across tissues and between individuals, highlighting the importance of cell-type specific and post-transcriptional processes in shaping mitochondrial-encoded RNA levels. Using whole-genome genetic data we identify 64 nuclear loci associated with expression levels of 14 genes encoded in the mitochondrial genome, including missense variants within genes involved in mitochondrial function (TBRG4, MTPAP and LONP1), implicating genetic mechanisms that act in trans across the two genomes. We replicate ~21% of associations with independent tissue-matched datasets and find genetic variants linked to these nuclear loci that are associated with cardio-metabolic phenotypes and Vitiligo, supporting a potential role for variable mitochondrial-encoded gene expression in complex disease.

CERENKOV2: improved detection of functional noncoding SNPs using data-space geometric features

BACKGROUND

We previously reported on CERENKOV, an approach for identifying regulatory single nucleotide polymorphisms (rSNPs) that is based on 246 annotation features. CERENKOV uses the xgboost classifier and is designed to be used to find causal noncoding SNPs in loci identified by genome-wide association studies (GWAS). We reported that CERENKOV has state-of-the-art performance (by two traditional measures and a novel GWAS-oriented measure, AVGRANK) in a comparison to nine other tools for identifying functional noncoding SNPs, using a comprehensive reference SNP set (OSU17, 15,331 SNPs). Given that SNPs are grouped within loci in the reference SNP set and given the importance of the data-space manifold geometry for machine-learning model selection, we hypothesized that within-locus inter-SNP distances would have class-based distributional biases that could be exploited to improve rSNP recognition accuracy. We thus defined an intralocus SNP "radius" as the average data-space distance from a SNP to the other intralocus neighbors, and explored radius likelihoods for five distance measures.

RESULTS

We expanded the set of reference SNPs to 39,083 (the OSU18 set) and extracted CERENKOV SNP feature data. We computed radius empirical likelihoods and likelihood densities for rSNPs and control SNPs, and found significant likelihood differences between rSNPs and control SNPs. We fit parametric models of likelihood distributions for five different distance measures to obtain ten log-likelihood features that we combined with the 248-dimensional CERENKOV feature matrix. On the OSU18 SNP set, we measured the classification accuracy of CERENKOV with and without the new distance-based features, and found that the addition of distance-based features significantly improves rSNP recognition performance as measured by AUPVR, AUROC, and AVGRANK. Along with feature data for the OSU18 set, the software code for extracting the base feature matrix, estimating ten distance-based likelihood ratio features, and scoring candidate causal SNPs, are released as open-source software CERENKOV2.

CONCLUSIONS

Accounting for the locus-specific geometry of SNPs in data-space significantly improved the accuracy with which noncoding rSNPs can be computationally identified.

Heritability estimation and differential analysis of count data with generalized linear mixed models in genomic sequencing studies

Motivation

Genomic sequencing studies, including RNA sequencing and bisulfite sequencing studies, are becoming increasingly common and increasingly large. Large genomic sequencing studies open doors for accurate molecular trait heritability estimation and powerful differential analysis. Heritability estimation and differential analysis in sequencing studies requires the development of statistical methods that can properly account for the count nature of the sequencing data and that are computationally efficient for large datasets.

Results

Here, we develop such a method, PQLseq (Penalized Quasi-Likelihood for sequencing count data), to enable effective and efficient heritability estimation and differential analysis using the generalized linear mixed model framework. With extensive simulations and comparisons to previous methods, we show that PQLseq is the only method currently available that can produce unbiased heritability estimates for sequencing count data. In addition, we show that PQLseq is well suited for differential analysis in large sequencing studies, providing calibrated type I error control and more power compared to the standard linear mixed model methods. Finally, we apply PQLseq to perform gene expression heritability estimation and differential expression analysis in a large RNA sequencing study in the Hutterites.

Availability and implementation

PQLseq is implemented as an R package with source code freely available at www.xzlab.org/software.html and https://cran.r-project.org/web/packages/PQLseq/index.html.

Supplementary information

Supplementary data are available at Bioinformatics online.

Identification of Alzheimer's Disease-Related Genes Based on Data Integration Method

Alzheimer disease (AD) is the fourth major cause of death in the elderly following cancer, heart disease and cerebrovascular disease. Finding candidate causal genes can help in the design of Gene targeted drugs and effectively reduce the risk of the disease. Complex diseases such as AD are usually caused by multiple genes. The Genome-wide association study (GWAS), has identified the potential genetic variants for most diseases. However, because of linkage disequilibrium (LD), it is difficult to identify the causative mutations that directly cause diseases. In this study, we combined expression quantitative trait locus (eQTL) studies with the GWAS, to comprehensively define the genes that cause Alzheimer disease. The method used was the Summary Mendelian randomization (SMR), which is a novel method to integrate summarized data. Two GWAS studies and five eQTL studies were referenced in this paper. We found several candidate SNPs that have a strong relationship with AD. Most of these SNPs overlap in different data sets, providing relatively strong reliability. We also explain the function of the novel AD-related genes we have discovered.

Genetic effects on promoter usage are highly context-specific and contribute to complex traits

Genetic variants regulating RNA splicing and transcript usage have been implicated in both common and rare diseases. Although transcript usage quantitative trait loci (tuQTLs) have been mapped across multiple cell types and contexts, it is challenging to distinguish between the main molecular mechanisms controlling transcript usage: promoter choice, splicing and 3' end choice. Here, we analysed RNA-seq data from human macrophages exposed to three inflammatory and one metabolic stimulus. In addition to conventional gene-level and transcript-level analyses, we also directly quantified promoter usage, splicing and 3' end usage. We found that promoters, splicing and 3' ends were predominantly controlled by independent genetic variants enriched in distinct genomic features. Promoter usage QTLs were also 50% more likely to be context-specific than other tuQTLs and constituted 25% of the transcript-level colocalisations with complex traits. Thus, promoter usage might be an underappreciated molecular mechanism mediating complex trait associations in a context-specific manner.

GIGSEA: genotype imputed gene set enrichment analysis using GWAS summary level data

Summary

level data of GWAS becomes increasingly important in post-GWAS data mining. Here, we present GIGSEA (Genotype Imputed Gene Set Enrichment Analysis), a novel method that uses GWAS summary statistics and eQTL to infer differential gene expression and interrogate gene set enrichment for the trait-associated SNPs. By incorporating empirical eQTL of the disease relevant tissue, GIGSEA naturally accounts for factors such as gene size, gene boundary, SNP distal regulation and multiple-marker regulation. The weighted linear regression model was used to perform the enrichment test, properly adjusting for imputation accuracy, model incompleteness and redundancy in different gene sets. The significance level of enrichment is assessed by the permutation test, where matrix operation was employed to dramatically improve computation speed. GIGSEA has appropriate type I error rates, and discovers the plausible biological findings on the real data set.

Availability and implementation

GIGSEA is implemented in R, and freely available at www.github.com/zhushijia/GIGSEA.

Supplementary information

Supplementary data are available at Bioinformatics online.

Gene expression variation and parental allele inheritance in a Xiphophorus interspecies hybridization model

Understanding the genetic mechanisms underlying segregation of phenotypic variation through successive generations is important for understanding physiological changes and disease risk. Tracing the etiology of variation in gene expression enables identification of genetic interactions, and may uncover molecular mechanisms leading to the phenotypic expression of a trait, especially when utilizing model organisms that have well-defined genetic lineages. There are a plethora of studies that describe relationships between gene expression and genotype, however, the idea that global variations in gene expression are also controlled by genotype remains novel. Despite the identification of loci that control gene expression variation, the global understanding of how genome constitution affects trait variability is unknown. To study this question, we utilized Xiphophorus fish of different, but tractable genetic backgrounds (inbred, F₁ interspecies hybrids, and backcross hybrid progeny), and measured each individual’s gene expression concurrent with the degrees of inter-individual expression variation. We found, (a) F₁ interspecies hybrids exhibited less variability than inbred animals, indicting gene expression variation is not affected by the fraction of heterozygous loci within an individual genome, and (b), that mixing genotypes in backcross populations led to higher levels of gene expression variability, supporting the idea that expression variability is caused by heterogeneity of genotypes of cis or trans loci. In conclusion, heterogeneity of genotype, introduced by inheritance of different alleles, accounts for the largest effects on global phenotypical variability.

Phenotypical variability is a multi-factorial phenomenon. Although it has been shown that inheriting certain gene is associated with lower phenotypical variability, how genome complexity affect phenotypical variability is still unclear. To study this question, we used inbred Xiphophorus fish, backcross interspecies hybrids, and F₁ interspecies hybrids between select Xiphophorus species to model genetic composition with minimum, medium, and maximum heterozygosity respectively, and measured their global gene expression variability. We found gene expression variation is not affected by the percentage of heterozygous loci in individual genome, but instead related to heterogeneity of genotype at local or remote loci.

An integrative approach for building personalized gene regulatory networks for precision medicine

Only a small fraction of patients respond to the drug prescribed to treat their disease, which means that most are at risk of unnecessary exposure to side effects through ineffective drugs. This inter-individual variation in drug response is driven by differences in gene interactions caused by each patient's genetic background, environmental exposures, and the proportions of specific cell types involved in disease. These gene interactions can now be captured by building gene regulatory networks, by taking advantage of RNA velocity (the time derivative of the gene expression state), the ability to study hundreds of thousands of cells simultaneously, and the falling price of single-cell sequencing. Here, we propose an integrative approach that leverages these recent advances in single-cell data with the sensitivity of bulk data to enable the reconstruction of personalized, cell-type- and context-specific gene regulatory networks. We expect this approach will allow the prioritization of key driver genes for specific diseases and will provide knowledge that opens new avenues towards improved personalized healthcare.

Evidence for Weak Selective Constraint on Human Gene Expression

Variation in human complex traits is connected to variation in gene expression, and selection on complex traits can be reflected in selection on gene expression. Here, Glassberg and Gao et al. analyze polymorphic....

Gene expression variation is a major contributor to phenotypic variation in human complex traits. Selection on complex traits may therefore be reflected in constraint on gene expression. Here, we explore the effects of stabilizing selection on cis-regulatory genetic variation in humans. We analyze patterns of expression variation at copy number variants and find evidence for selection against large increases in gene expression. Using allele-specific expression (ASE) data, we further show evidence of selection against smaller-effect variants. We estimate that, across all genes, singletons in a sample of 122 individuals have ∼2.2× greater effects on expression variation than the average variant across allele frequencies. Despite their increased effect size relative to common variants, we estimate that singletons in the sample studied explain, on average, only 5% of the heritability of gene expression from cis-regulatory variants. Finally, we show that genes depleted for loss-of-function variants are also depleted for cis-eQTLs and have low levels of allelic imbalance, confirming tighter constraint on the expression levels of these genes. We conclude that constraint on gene expression is present, but has relatively weak effects on most cis-regulatory variants, thus permitting high levels of gene-regulatory genetic variation.

False positives in trans-eQTL and co-expression analyses arising from RNA-sequencing alignment errors

Sequence similarity among distinct genomic regions can lead to errors in alignment of short reads from next-generation sequencing. While this is well known, the downstream consequences of misalignment have not been fully characterized.  We assessed the potential for incorrect alignment of RNA-sequencing reads to cause false positives in both gene expression quantitative trait locus (eQTL) and co-expression analyses. Trans-eQTLs identified from human RNA-sequencing studies appeared to be particularly affected by this phenomenon, even when only uniquely aligned reads are considered. Over 75% of trans-eQTLs using a standard pipeline occurred between regions of sequence similarity and therefore could be due to alignment errors. Further, associations due to mapping errors are likely to misleadingly replicate between studies. To help address this problem, we quantified the potential for "cross-mapping'' to occur between every pair of annotated genes in the human genome. Such cross-mapping data can be used to filter or flag potential false positives in both trans-eQTL and co-expression analyses. Such filtering substantially alters the detection of significant associations and can have an impact on the assessment of false discovery rate, functional enrichment, and replication for RNA-sequencing association studies.

False positives in trans-eQTL and co-expression analyses arising from RNA-sequencing alignment errors

Sequence similarity among distinct genomic regions can lead to errors in alignment of short reads from next-generation sequencing. While this is well known, the downstream consequences of misalignment have not been fully characterized.  We assessed the potential for incorrect alignment of RNA-sequencing reads to cause false positives in both gene expression quantitative trait locus (eQTL) and co-expression analyses. Trans-eQTLs identified from human RNA-sequencing studies appeared to be particularly affected by this phenomenon, even when only uniquely aligned reads are considered. Over 75\% of trans-eQTLs using a standard pipeline occurred between regions of sequence similarity and therefore could be due to alignment errors. Further, associations due to mapping errors are likely to misleadingly replicate between studies. To help address this problem, we quantified the potential for "cross-mapping'' to occur between every pair of annotated genes in the human genome. Such cross-mapping data can be used to filter or flag potential false positives in both trans-eQTL and co-expression analyses. Such filtering substantially alters the detection of significant associations and can have an impact on the assessment of false discovery rate, functional enrichment, and replication for RNA-sequencing association studies.

High-resolution genetic mapping of putative causal interactions between regions of open chromatin

Physical interaction of regulatory elements in three-dimensional space poses a challenge for studies of disease because non-coding risk variants may be great distances from the genes they regulate. Experimental methods to capture these interactions, such as chromosome conformation capture, usually cannot assign causal direction of effect between regulatory elements, an important component of fine-mapping studies. We developed a Bayesian hierarchical approach that uses two-stage least squares and applied it to an ATAC-seq (assay for transposase-accessible chromatin using sequencing) data set from 100 individuals, to identify over 15,000 high-confidence causal interactions. Most (60%) interactions occurred over <20 kb, where chromosome conformation capture-based methods perform poorly. For a fraction of loci, we identified a single variant that alters accessibility across multiple regions, and experimentally validated the BLK locus, which is associated with multiple autoimmune diseases, using CRISPR genome editing. Our study highlights how association genetics of chromatin state is a powerful approach for identifying interactions between regulatory elements.

Bid maintains mitochondrial cristae structure and function and protects against cardiac disease in an integrative genomics study

Bcl-2 family proteins reorganize mitochondrial membranes during apoptosis, to form pores and rearrange cristae. In vitro and in vivo analysis integrated with human genetics reveals a novel homeostatic mitochondrial function for Bcl-2 family protein Bid. Loss of full-length Bid results in apoptosis-independent, irregular cristae with decreased respiration. Bid-/- mice display stress-induced myocardial dysfunction and damage. A gene-based approach applied to a biobank, validated in two independent GWAS studies, reveals that decreased genetically determined BID expression associates with myocardial infarction (MI) susceptibility. Patients in the bottom 5% of the expression distribution exhibit >4 fold increased MI risk. Carrier status with nonsynonymous variation in Bid’s membrane binding domain, Bid^M148T, associates with MI predisposition. Furthermore, Bid but not Bid^M148T associates with Mcl-1^Matrix, previously implicated in cristae stability; decreased MCL-1 expression associates with MI. Our results identify a role for Bid in homeostatic mitochondrial cristae reorganization, that we link to human cardiac disease.

Cells contain specialized structures called mitochondria, which help to convert fuel into energy. These tiny energy factories have a unique double membrane, with a smooth outer and a folded inner lining. The folds, called cristae, provide a scaffold for the molecular machinery that produces chemical energy that the cell can use. The cristae are dynamic, and can change shape, condensing to increase energy output.

Mitochondria also play a role in cell death. In certain situations, cristae can widen and release the proteins held within their folds. This can trigger a program of self-destruction in the cell. A family of proteins called Bcl-2 control such a ‘programmed cell death’ through the release of mitochondrial proteins. Some family members, including a protein called Bid, can reorganize cristae to regulate this cell-death program. When cells die, Bid proteins that had been split move to the mitochondria. But, even when cells are healthy, Bid molecules that are intact are always there, suggesting that this form of the protein may have another purpose.

To investigate this further, Salisbury-Ruf, Bertram et al. used mice with Bid, and mice that lacked the protein. Without Bid, cells – including heart cells – struggled to work properly and used less oxygen than their normal counterparts. A closer look using electron microscopy revealed abnormalities in the cristae. However, adding ‘intact’ Bid proteins back in to the deficient cells restored them to normal.

Moreover, without Bid, the mice hearts were less able to respond to an increased demand for energy. This decreased their performance and caused the formation of scars in the heart muscle called fibrosis, similar to a pattern observed in human patients following a heart attack.

DNA data from an electronic health record database revealed a link between low levels of Bid genes and heart attack in humans, which was confirmed in further studies. In addition, a specific mutation in the Bid gene was found to affect its ability to regulate the formation of proper cristae.

Combining evidence from mice with human genetics revealed new information about heart diseases. Mitochondrial health may be affected by a combination of specific variations in genes and changes in the Bid protein, which could affect heart attack risk. Understanding more about this association could help to identify and potentially reduce certain risk factors for heart attack.

FIRE: functional inference of genetic variants that regulate gene expression

Motivation

Interpreting genetic variation in noncoding regions of the genome is an important challenge for personal genome analysis. One mechanism by which noncoding single nucleotide variants (SNVs) influence downstream phenotypes is through the regulation of gene expression. Methods to predict whether or not individual SNVs are likely to regulate gene expression would aid interpretation of variants of unknown significance identified in whole-genome sequencing studies.

Results

We developed FIRE (Functional Inference of Regulators of Expression), a tool to score both noncoding and coding SNVs based on their potential to regulate the expression levels of nearby genes. FIRE consists of 23 random forests trained to recognize SNVs in cis-expression quantitative trait loci (cis-eQTLs) using a set of 92 genomic annotations as predictive features. FIRE scores discriminate cis-eQTL SNVs from non-eQTL SNVs in the training set with a cross-validated area under the receiver operating characteristic curve (AUC) of 0.807, and discriminate cis-eQTL SNVs shared across six populations of different ancestry from non-eQTL SNVs with an AUC of 0.939. FIRE scores are also predictive of cis-eQTL SNVs across a variety of tissue types.

Availability and implementation

FIRE scores for genome-wide SNVs in hg19/GRCh37 are available for download at https://sites.google.com/site/fireregulatoryvariation/.

Contact

nilah@stanford.edu.

Supplementary information

Supplementary data are available at Bioinformatics online.

Genetic architecture of gene expression traits across diverse populations

For many complex traits, gene regulation is likely to play a crucial mechanistic role. How the genetic architectures of complex traits vary between populations and subsequent effects on genetic prediction are not well understood, in part due to the historical paucity of GWAS in populations of non-European ancestry. We used data from the MESA (Multi-Ethnic Study of Atherosclerosis) cohort to characterize the genetic architecture of gene expression within and between diverse populations. Genotype and monocyte gene expression were available in individuals with African American (AFA, n = 233), Hispanic (HIS, n = 352), and European (CAU, n = 578) ancestry. We performed expression quantitative trait loci (eQTL) mapping in each population and show genetic correlation of gene expression depends on shared ancestry proportions. Using elastic net modeling with cross validation to optimize genotypic predictors of gene expression in each population, we show the genetic architecture of gene expression for most predictable genes is sparse. We found the best predicted gene in each population, TACSTD2 in AFA and CHURC1 in CAU and HIS, had similar prediction performance across populations with R2 > 0.8 in each population. However, we identified a subset of genes that are well-predicted in one population, but poorly predicted in another. We show these differences in predictive performance are due to allele frequency differences between populations. Using genotype weights trained in MESA to predict gene expression in independent populations showed that a training set with ancestry similar to the test set is better at predicting gene expression in test populations, demonstrating an urgent need for diverse population sampling in genomics. Our predictive models and performance statistics in diverse cohorts are made publicly available for use in transcriptome mapping methods at https://github.com/WheelerLab/DivPop.

Genetics of trans-regulatory variation in gene expression

Heritable variation in gene expression forms a crucial bridge between genomic variation and the biology of many traits. However, most expression quantitative trait loci (eQTLs) remain unidentified. We mapped eQTLs by transcriptome sequencing in 1012 yeast segregants. The resulting eQTLs accounted for over 70% of the heritability of mRNA levels, allowing comprehensive dissection of regulatory variation. Most genes had multiple eQTLs. Most expression variation arose from trans-acting eQTLs distant from their target genes. Nearly all trans-eQTLs clustered at 102 hotspot locations, some of which influenced the expression of thousands of genes. Fine-mapped hotspot regions were enriched for transcription factor genes. While most genes had a local eQTL, most of these had no detectable effects on the expression of other genes in trans. Hundreds of non-additive genetic interactions accounted for small fractions of expression variation. These results reveal the complexity of genetic influences on transcriptome variation in unprecedented depth and detail.

Genetic determinants of co-accessible chromatin regions in activated T cells across humans

Over 90% of genetic variants associated with complex human traits map to non-coding regions, but little is understood about how they modulate gene regulation in health and disease. One possible mechanism is that genetic variants affect the activity of one or more cis-regulatory elements leading to gene expression variation in specific cell types. To identify such cases, we analyzed ATAC-seq and RNA-seq profiles from stimulated primary CD4⁺ T cells in up to 105 healthy donors. We found that regions of accessible chromatin (ATAC-peaks) are co-accessible at kilobase and megabase resolution, consistent with the three-dimensional chromatin organization measured by in situ Hi-C in T cells. Fifteen percent of genetic variants located within ATAC-peaks affected the accessibility of the corresponding peak (local-ATAC-QTLs). Local-ATAC-QTLs have the largest effects on co-accessible peaks, are associated with gene expression and are enriched for autoimmune disease variants. Our results provide insights into how natural genetic variants modulate cis-regulatory elements, in isolation or in concert, to influence gene expression.

Using an atlas of gene regulation across 44 human tissues to inform complex disease- and trait-associated variation

We apply integrative approaches to expression quantitative loci (eQTLs) from 44 tissues from the Genotype-Tissue Expression project and genome-wide association study data. About 60% of known trait-associated loci are in linkage disequilibrium with a cis-eQTL, over half of which were not found in previous large-scale whole blood studies. Applying polygenic analyses to metabolic, cardiovascular, anthropometric, autoimmune, and neurodegenerative traits, we find that eQTLs are significantly enriched for trait associations in relevant pathogenic tissues and explain a substantial proportion of the heritability (40-80%). For most traits, tissue-shared eQTLs underlie a greater proportion of trait associations, although tissue-specific eQTLs have a greater contribution to some traits, such as blood pressure. By integrating information from biological pathways with eQTL target genes and applying a gene-based approach, we validate previously implicated causal genes and pathways, and propose new variant and gene associations for several complex traits, which we replicate in the UK BioBank and BioVU.

Host genetics and microbiome associations through the lens of genome wide association studies

Recent studies indicate that the gut microbiome is partially heritable, motivating the need to investigate microbiome-host genome associations via microbial genome-wide association studies (mGWAS). Existing mGWAS demonstrate that microbiome-host genotype associations are typically weak and are spread across multiple variants, similar to associations often observed in genome-wide association studies (GWAS) of complex traits. Here we reconsider mGWAS by viewing them through the lens of GWAS, and demonstrate that there are striking similarities between the challenges and pitfalls faced by the two study designs. We further advocate the mGWAS community to adopt three key lessons learned over the history of GWAS: firstly, adopting uniform data and reporting formats to facilitate replication and meta-analysis efforts; secondly, enforcing stringent statistical criteria to reduce the number of false positive findings; and thirdly, considering the microbiome and the host genome as distinct entities, rather than studying different taxa and single nucleotide polymorphism (SNPs) separately. Finally, we anticipate that mGWAS sample sizes will have to increase by orders of magnitude to reproducibly associate the host genome with the gut microbiome.

Principled multi-omic analysis reveals gene regulatory mechanisms of phenotype variation

Recent studies have analyzed large-scale data sets of gene expression to identify genes associated with interindividual variation in phenotypes ranging from cancer subtypes to drug sensitivity, promising new avenues of research in personalized medicine. However, gene expression data alone is limited in its ability to reveal cis-regulatory mechanisms underlying phenotypic differences. In this study, we develop a new probabilistic model, called pGENMi, that integrates multi-omic data to investigate the transcriptional regulatory mechanisms underlying interindividual variation of a specific phenotype—that of cell line response to cytotoxic treatment. In particular, pGENMi simultaneously analyzes genotype, DNA methylation, gene expression, and transcription factor (TF)-DNA binding data, along with phenotypic measurements, to identify TFs regulating the phenotype. It does so by combining statistical information about expression quantitative trait loci (eQTLs) and expression-correlated methylation marks (eQTMs) located within TF binding sites, as well as observed correlations between gene expression and phenotype variation. Application of pGENMi to data from a panel of lymphoblastoid cell lines treated with 24 drugs, in conjunction with ENCODE TF ChIP data, yielded a number of known as well as novel (TF, Drug) associations. Experimental validations by TF knockdown confirmed 41% of the predicted and tested associations, compared to a 12% confirmation rate of tested nonassociations (controls). An extensive literature survey also corroborated 62% of the predicted associations above a stringent threshold. Moreover, associations predicted only when combining eQTL and eQTM data showed higher precision compared to an eQTL-only or eQTM-only analysis using pGENMi, further demonstrating the value of multi-omic integrative analysis.

From gene networks to drugs: systems pharmacology approaches for AUD

The alcohol research field has amassed an impressive number of gene expression datasets spanning key brain areas for addiction, species (humans as well as multiple animal models), and stages in the addiction cycle (binge/intoxication, withdrawal/negative effect, and preoccupation/anticipation). These data have improved our understanding of the molecular adaptations that eventually lead to dysregulation of brain function and the chronic, relapsing disorder of addiction. Identification of new medications to treat alcohol use disorder (AUD) will likely benefit from the integration of genetic, genomic, and behavioral information included in these important datasets. Systems pharmacology considers drug effects as the outcome of the complex network of interactions a drug has rather than a single drug-molecule interaction. Computational strategies based on this principle that integrate gene expression signatures of pharmaceuticals and disease states have shown promise for identifying treatments that ameliorate disease symptoms (called in silico gene mapping or connectivity mapping). In this review, we suggest that gene expression profiling for in silico mapping is critical to improve drug repurposing and discovery for AUD and other psychiatric illnesses. We highlight studies that successfully apply gene mapping computational approaches to identify or repurpose pharmaceutical treatments for psychiatric illnesses. Furthermore, we address important challenges that must be overcome to maximize the potential of these strategies to translate to the clinic and improve healthcare outcomes.

Omics research project on prospective cohort studies from the Tohoku Medical Megabank Project

Population-based prospective cohort studies are indispensable for modern medical research as they provide important knowledge on the influences of many kinds of genetic and environmental factors on the cause of disease. Although traditional cohort studies are mainly conducted using questionnaires and physical examinations, modern cohort studies incorporate omics and genomic approaches to obtain comprehensive physical information, including genetic information. Here, we report the design and midterm results of multi-omics analysis on population-based prospective cohort studies from the Tohoku Medical Megabank (TMM) Project. We have incorporated genomic and metabolomic studies in the TMM cohort study as both metabolome and genome analyses are suitable for high-throughput analysis of large-scale cohort samples. Moreover, an association study between the metabolome and genome show that metabolites are an important intermediate phenotype connecting genetic and lifestyle factors to physical and pathologic phenotypes. We apply our metabolome and genome analyses to large-scale cohort samples in the following studies.

The dorsoventral patterning of Musca domestica embryos: insights into BMP/Dpp evolution from the base of the lower cyclorraphan flies

Background

In the last few years, accumulated information has indicated that the evolution of an extra-embryonic membrane in dipterans was accompanied by changes in the gene regulatory network controlled by the BMP/Dpp pathway, which is responsible for dorsal patterning in these insects. However, only comparative analysis of gene expression levels between distant species with two extra-embryonic membranes, like A. gambiae or C. albipunctata, and D. melanogaster, has been conducted. Analysis of gene expression in ancestral species, which evolved closer to the amnioserosa origin, could provide new insights into the evolution of dorsoventral patterning in dipterans.

Results

Here we describe the spatial expression of several key and downstream elements of the Dpp pathway and show the compared patterns of expression between Musca and Drosophila embryos, both dipterans with amnioserosa. Most of the analyzed gene showed a high degree of expression conservation, however, we found several differences in the gene expression pattern of M. domestica orthologs for sog and tolloid. Bioinformatics analysis of the promoter of both genes indicated that the variations could be related to the gain of several binding sites for the transcriptional factor Dorsal in the Md.tld promoter and Snail in the Md.sog enhancer. These altered expressions could explain the unclear formation of the pMad gradient in the M. domestica embryo, compared to the formation of the gradient in D. melanogaster.

Conclusion

Gene expression changes during the dorsal–ventral patterning in insects contribute to the differentiation of extra-embryonic tissues as a consequence of changes in the gene regulatory network controlled by BMP/Dpp. In this work, in early M. domestica embryos, we identified the expression pattern of several genes members involved in the dorsoventral specification of the embryo. We believe that these data can contribute to understanding the evolution of the BMP/Dpp pathway, the regulation of BMP ligands, and the formation of a Dpp gradient in higher cyclorraphan flies.

Electronic supplementary material

The online version of this article (10.1186/s13227-018-0102-5) contains supplementary material, which is available to authorized users.

A molecular network of the aging human brain provides insights into the pathology and cognitive decline of Alzheimer's disease

There is a need for new therapeutic targets with which to prevent Alzheimer’s disease (AD), a major contributor to aging-related cognitive decline. Here, we report the construction and validation of a molecular network of the aging human frontal cortex. Using RNA sequence data from 478 individuals, we first build a molecular network using modules of coexpressed genes and then relate these modules to AD and its neuropathologic and cognitive endophenotypes. We confirm these associations in two independent AD datasets as well as in epigenomic data. We also illustrate the use of the network in prioritizing amyloid-associated genes for in vitro validation in human neurons and astrocytes. These analyses based on unique cohorts enable us to resolve the role of distinct cortical modules that have a direct effect on the accumulation of AD pathology from those that have a direct effect on cognitive decline, exemplifying a network approach to complex diseases.

Systems biology analysis of RNA sequencing data from the aging human cortex identifies a molecular network which prioritizes groups of genes that influence cognitive decline or neuropathology in Alzheimer’s disease.

Putative bovine topological association domains and CTCF binding motifs can reduce the search space for causative regulatory variants of complex traits

BACKGROUND

Topological association domains (TADs) are chromosomal domains characterised by frequent internal DNA-DNA interactions. The transcription factor CTCF binds to conserved DNA sequence patterns called CTCF binding motifs to either prohibit or facilitate chromosomal interactions. TADs and CTCF binding motifs control gene expression, but they are not yet well defined in the bovine genome. In this paper, we sought to improve the annotation of bovine TADs and CTCF binding motifs, and assess whether the new annotation can reduce the search space for cis-regulatory variants.

RESULTS

We used genomic synteny to map TADs and CTCF binding motifs from humans, mice, dogs and macaques to the bovine genome. We found that our mapped TADs exhibited the same hallmark properties of those sourced from experimental data, such as housekeeping genes, transfer RNA genes, CTCF binding motifs, short interspersed elements, H3K4me3 and H3K27ac. We showed that runs of genes with the same pattern of allele-specific expression (ASE) (either favouring paternal or maternal allele) were often located in the same TAD or between the same conserved CTCF binding motifs. Analyses of variance showed that when averaged across all bovine tissues tested, TADs explained 14% of ASE variation (standard deviation, SD: 0.056), while CTCF explained 27% (SD: 0.078). Furthermore, we showed that the quantitative trait loci (QTLs) associated with gene expression variation (eQTLs) or ASE variation (aseQTLs), which were identified from mRNA transcripts from 141 lactating cows' white blood and milk cells, were highly enriched at putative bovine CTCF binding motifs. The linearly-furthermost, and most-significant aseQTL and eQTL for each genic target were located within the same TAD as the gene more often than expected (Chi-Squared test P-value < 0.001).

CONCLUSIONS

Our results suggest that genomic synteny can be used to functionally annotate conserved transcriptional components, and provides a tool to reduce the search space for causative regulatory variants in the bovine genome.

Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics

Scalable, integrative methods to understand mechanisms that link genetic variants with phenotypes are needed. Here we derive a mathematical expression to compute PrediXcan (a gene mapping approach) results using summary data (S-PrediXcan) and show its accuracy and general robustness to misspecified reference sets. We apply this framework to 44 GTEx tissues and 100+ phenotypes from GWAS and meta-analysis studies, creating a growing public catalog of associations that seeks to capture the effects of gene expression variation on human phenotypes. Replication in an independent cohort is shown. Most of the associations are tissue specific, suggesting context specificity of the trait etiology. Colocalized significant associations in unexpected tissues underscore the need for an agnostic scanning of multiple contexts to improve our ability to detect causal regulatory mechanisms. Monogenic disease genes are enriched among significant associations for related traits, suggesting that smaller alterations of these genes may cause a spectrum of milder phenotypes.

An Osteoporosis Risk SNP at 1p36.12 Acts as an Allele-Specific Enhancer to Modulate LINC00339 Expression via Long-Range Loop Formation

Genome-wide association studies (GWASs) have reproducibly associated variants within intergenic regions of 1p36.12 locus with osteoporosis, but the functional roles underlying these noncoding variants are unknown. Through an integrative functional genomic and epigenomic analyses, we prioritized rs6426749 as a potential causal SNP for osteoporosis at 1p36.12. Dual-luciferase assay and CRISPR/Cas9 experiments demonstrate that rs6426749 acts as a distal allele-specific enhancer regulating expression of a lncRNA (LINC00339) (∼360 kb) via long-range chromatin loop formation and that this loop is mediated by CTCF occupied near rs6426749 and LINC00339 promoter region. Specifically, rs6426749-G allele can bind transcription factor TFAP2A, which efficiently elevates the enhancer activity and increases LINC00339 expression. Downregulation of LINC00339 significantly increases the expression of CDC42 in osteoblast cells, which is a pivotal regulator involved in bone metabolism. Our study provides mechanistic insight into how a noncoding SNP affects osteoporosis by long-range interaction, a finding that could indicate promising therapeutic targets for osteoporosis.

Using Droplet Digital PCR to Analyze Allele-Specific RNA Expression

Genome-wide association studies have discovered thousands of common alleles that associate with human phenotypes and disease. Many of these variants are in non-protein-coding (regulatory) regions and are believed to affect phenotypes by modifying gene expression. In any organism with a diploid genome, such as humans, measuring the expression of each allele of a gene provides a well-controlled way to identify allelic influences on that gene's expression. Here, we describe a protocol for precisely measuring the allele-specific expression of individual genes. This method targets the nucleotide differences between the two alleles of a gene within an individual and measures the "allelic skew," the extent to which one allele is expressed more than the other. We cover the design of effective assays, the optimization of reactions, and the interpretation of the resulting data.

A Mixed-Effects Model for Powerful Association Tests in Integrative Functional Genomics

Genome-wide association studies (GWASs) have successfully identified thousands of genetic variants for many complex diseases; however, these variants explain only a small fraction of the heritability. Recently, genetic association studies that leverage external transcriptome data have received much attention and shown promise for discovering novel variants. One such approach, PrediXcan, is to use predicted gene expression through genetic regulation. However, there are limitations in this approach. The predicted gene expression may be biased, resulting from regularized regression applied to moderately sample-sized reference studies. Further, some variants can individually influence disease risk through alternative functional mechanisms besides expression. Thus, testing only the association of predicted gene expression as proposed in PrediXcan will potentially lose power. To tackle these challenges, we consider a unified mixed effects model that formulates the association of intermediate phenotypes such as imputed gene expression through fixed effects, while allowing residual effects of individual variants to be random. We consider a set-based score testing framework, MiST (mixed effects score test), and propose two data-driven combination approaches to jointly test for the fixed and random effects. We establish the asymptotic distributions, which enable rapid calculation of p values for genome-wide analyses, and provide p values for fixed and random effects separately to enhance interpretability over GWASs. Extensive simulations demonstrate that our approaches are more powerful than existing ones. We apply our approach to a large-scale GWAS of colorectal cancer and identify two genes, POU5F1B and ATF1, which would have otherwise been missed by PrediXcan, after adjusting for all known loci.

Bayesian nonparametric discovery of isoforms and individual specific quantification

Most human protein-coding genes can be transcribed into multiple distinct mRNA isoforms. These alternative splicing patterns encourage molecular diversity, and dysregulation of isoform expression plays an important role in disease etiology. However, isoforms are difficult to characterize from short-read RNA-seq data because they share identical subsequences and occur in different frequencies across tissues and samples. Here, we develop biisq, a Bayesian nonparametric model for isoform discovery and individual specific quantification from short-read RNA-seq data. biisq does not require isoform reference sequences but instead estimates an isoform catalog shared across samples. We use stochastic variational inference for efficient posterior estimates and demonstrate superior precision and recall for simulations compared to state-of-the-art isoform reconstruction methods. biisq shows the most gains for low abundance isoforms, with 36% more isoforms correctly inferred at low coverage versus a multi-sample method and 170% more versus single-sample methods. We estimate isoforms in the GEUVADIS RNA-seq data and validate inferred isoforms by associating genetic variants with isoform ratios.

Alternative splicing leads to transcript isoform diversity. Here, Aguiar et al. develop biisq, a Bayesian nonparametric approach to discover and quantify isoforms from RNA-seq data.

A Common Type 2 Diabetes Risk Variant Potentiates Activity of an Evolutionarily Conserved Islet Stretch Enhancer and Increases C2CD4A and C2CD4B Expression

Genome-wide association studies (GWASs) and functional genomics approaches implicate enhancer disruption in islet dysfunction and type 2 diabetes (T2D) risk. We applied genetic fine-mapping and functional (epi)genomic approaches to a T2D- and proinsulin-associated 15q22.2 locus to identify a most likely causal variant, determine its direction of effect, and elucidate plausible target genes. Fine-mapping and conditional analyses of proinsulin levels of 8,635 non-diabetic individuals from the METSIM study support a single association signal represented by a cluster of 16 strongly associated (p < 10^-17) variants in high linkage disequilibrium (r² > 0.8) with the GWAS index SNP rs7172432. These variants reside in an evolutionarily and functionally conserved islet and β cell stretch or super enhancer; the most strongly associated variant (rs7163757, p = 3 × 10^-19) overlaps a conserved islet open chromatin site. DNA sequence containing the rs7163757 risk allele displayed 2-fold higher enhancer activity than the non-risk allele in reporter assays (p < 0.01) and was differentially bound by β cell nuclear extract proteins. Transcription factor NFAT specifically potentiated risk-allele enhancer activity and altered patterns of nuclear protein binding to the risk allele in vitro, suggesting that it could be a factor mediating risk-allele effects. Finally, the rs7163757 proinsulin-raising and T2D risk allele (C) was associated with increased expression of C2CD4B, and possibly C2CD4A, both of which were induced by inflammatory cytokines, in human islets. Together, these data suggest that rs7163757 contributes to genetic risk of islet dysfunction and T2D by increasing NFAT-mediated islet enhancer activity and modulating C2CD4B, and possibly C2CD4A, expression in (patho)physiologic states.

Understanding Neurodevelopmental Disorders: The Promise of Regulatory Variation in the 3'UTRome

Neurodevelopmental disorders have a strong genetic component, but despite widespread efforts, the specific genetic factors underlying these disorders remain undefined for a large proportion of affected individuals. Given the accessibility of exome sequencing, this problem has thus far been addressed from a protein-centric standpoint; however, protein-coding regions only make up ∼1% to 2% of the human genome. With the advent of whole genome sequencing we are in the midst of a paradigm shift as it is now possible to interrogate the entire sequence of the human genome (coding and noncoding) to fill in the missing heritability of complex disorders. These new technologies bring new challenges, as the number of noncoding variants identified per individual can be overwhelming, making it prudent to focus on noncoding regions of known function, for which the effects of variation can be predicted and directly tested to assess pathogenicity. The 3'UTRome is a region of the noncoding genome that perfectly fulfills these criteria and is of high interest when searching for pathogenic variation related to complex neurodevelopmental disorders. Herein, we review the regulatory roles of the 3'UTRome as binding sites for microRNAs or RNA binding proteins, or during alternative polyadenylation. We detail existing evidence that these regions contribute to neurodevelopmental disorders and outline strategies for identification and validation of novel putatively pathogenic variation in these regions. This evidence suggests that studying the 3'UTRome will lead to the identification of new risk factors, new candidate disease genes, and a better understanding of the molecular mechanisms contributing to neurodevelopmental disorders.

Gene-by-environment interactions in urban populations modulate risk phenotypes

Uncovering the interaction between genomes and the environment is a principal challenge of modern genomics and preventive medicine. While theoretical models are well defined, little is known of the G × E interactions in humans. We used an integrative approach to comprehensively assess the interactions between 1.6 million data points, encompassing a range of environmental exposures, health, and gene expression levels, coupled with whole-genome genetic variation. From ∼1000 individuals of a founder population in Quebec, we reveal a substantial impact of the environment on the transcriptome and clinical endophenotypes, overpowering that of genetic ancestry. Air pollution impacts gene expression and pathways affecting cardio-metabolic and respiratory traits, when controlling for genetic ancestry. Finally, we capture four expression quantitative trait loci that interact with the environment (air pollution). Our findings demonstrate how the local environment directly affects disease risk phenotypes and that genetic variation, including less common variants, can modulate individual’s response to environmental challenges.

Individuals with different genotypes may respond differently to environmental variation. Here, Favé et al. find substantial impacts of different environment exposures on the transcriptome and clinical endophenotypes when controlling for genetic ancestry by analyzing data from ∼1000 individuals from a founder population in Quebec.

Characterising cis-regulatory variation in the transcriptome of histologically normal and tumour-derived pancreatic tissues

OBJECTIVE

To elucidate the genetic architecture of gene expression in pancreatic tissues.

DESIGN

We performed expression quantitative trait locus (eQTL) analysis in histologically normal pancreatic tissue samples (n=95) using RNA sequencing and the corresponding 1000 genomes imputed germline genotypes. Data from pancreatic tumour-derived tissue samples (n=115) from The Cancer Genome Atlas were included for comparison.

RESULTS

We identified 38 615 cis-eQTLs (in 484 genes) in histologically normal tissues and 39 713 cis-eQTL (in 237 genes) in tumour-derived tissues (false discovery rate <0.1), with the strongest effects seen near transcriptional start sites. Approximately 23% and 42% of genes with significant cis-eQTLs appeared to be specific for tumour-derived and normal-derived tissues, respectively. Significant enrichment of cis-eQTL variants was noted in non-coding regulatory regions, in particular for pancreatic tissues (1.53-fold to 3.12-fold, p≤0.0001), indicating tissue-specific functional relevance. A common pancreatic cancer risk locus on 9q34.2 (rs687289) was associated with ABO expression in histologically normal (p=5.8×10-8) and tumour-derived (p=8.3×10-5) tissues. The high linkage disequilibrium between this variant and the O blood group generating deletion variant in ABO (exon 6) suggested that nonsense-mediated decay (NMD) of the 'O' mRNA might explain this finding. However, knockdown of crucial NMD regulators did not influence decay of the ABO 'O' mRNA, indicating that a gene regulatory element influenced by pancreatic cancer risk alleles may underlie the eQTL.

CONCLUSIONS

We have identified cis-eQTLs representing potential functional regulatory variants in the pancreas and generated a rich data set for further studies on gene expression and its regulation in pancreatic tissues.

Molecular-genetic characterization of common, noncoding UBASH3A variants associated with type 1 diabetes

Genome-wide association and fine-mapping studies have identified over 40 susceptibility regions for type 1 diabetes (T1D), a common autoimmune disease; however, most of the disease-associated variants are noncoding, and it remains a challenge to understand their biological contributions to T1D pathogenesis. One identified T1D risk locus is located at chromosome 21q22.3 where the most likely candidate gene is UBASH3A, a negative regulator of NF-κB signaling. Various noncoding variants in UBASH3A have been shown to be associated with T1D or other autoimmune diseases. Here we investigated four such SNPs-rs11203202, rs80054410, rs11203203, and rs1893592. We discovered a novel role for rs1893592 in T1D and showed that its minor allele protects against T1D. Our haplotype analysis identified three T1D-associated UBASH3A haplotypes, and revealed that risk for T1D is affected by additive effects of these four UBASH3A variants. In human primary CD4⁺ T cells, upon T-cell receptor stimulation, the minor allele of rs1893592 was associated with both a significant reduction in the overall mRNA levels of UBASH3A, and an increase in the proportion of a normally occurring, but low-abundant, UBASH3A transcript that retains intron-9 sequences and cannot produce full-length UBASH3A protein. This reduction in UBASH3A, as a consequence of the minor allele at rs1893592, resulted in increased secretion of IL-2, a key cytokine that is required for T-cell activation and function but is deficient in some T1D subjects. Our study provides new mechanistic insights into how rs1893592 affects T1D and autoimmunity, and how interactions between multiple T1D-associated, noncoding variants influence the disease risk.

The Expanding Landscape of Alternative Splicing Variation in Human Populations

Alternative splicing is a tightly regulated biological process by which the number of gene products for any given gene can be greatly expanded. Genomic variants in splicing regulatory sequences can disrupt splicing and cause disease. Recent developments in sequencing technologies and computational biology have allowed researchers to investigate alternative splicing at an unprecedented scale and resolution. Population-scale transcriptome studies have revealed many naturally occurring genetic variants that modulate alternative splicing and consequently influence phenotypic variability and disease susceptibility in human populations. Innovations in experimental and computational tools such as massively parallel reporter assays and deep learning have enabled the rapid screening of genomic variants for their causal impacts on splicing. In this review, we describe technological advances that have greatly increased the speed and scale at which discoveries are made about the genetic variation of alternative splicing. We summarize major findings from population transcriptomic studies of alternative splicing and discuss the implications of these findings for human genetics and medicine.

Evaluation of PrediXcan for prioritizing GWAS associations and predicting gene expression

Genome-wide association studies (GWAS) have been successful in facilitating the understanding of genetic architecture behind human diseases, but this approach faces many challenges. To identify disease-related loci with modest to weak effect size, GWAS requires very large sample sizes, which can be computational burdensome. In addition, the interpretation of discovered associations remains difficult. PrediXcan was developed to help address these issues. With built in SNP-expression models, PrediXcan is able to predict the expression of genes that are regulated by putative expression quantitative trait loci (eQTLs), and these predicted expression levels can then be used to perform gene-based association studies. This approach reduces the multiple testing burden from millions of variants down to several thousand genes. But most importantly, the identified associations can reveal the genes that are under regulation of eQTLs and consequently involved in disease pathogenesis. In this study, two of the most practical functions of PrediXcan were tested: 1) predicting gene expression, and 2) prioritizing GWAS results. We tested the prediction accuracy of PrediXcan by comparing the predicted and observed gene expression levels, and also looked into some potential influential factors and a filter criterion with the aim of improving PrediXcan performance. As for GWAS prioritization, predicted gene expression levels were used to obtain gene-trait associations, and background regions of significant associations were examined to decrease the likelihood of false positives. Our results showed that 1) PrediXcan predicted gene expression levels accurately for some but not all genes; 2) including more putative eQTLs into prediction did not improve the prediction accuracy; and 3) integrating predicted gene expression levels from the two PrediXcan whole blood models did not eliminate false positives. Still, PrediXcan was able to prioritize GWAS associations that were below the genome-wide significance threshold in GWAS, while retaining GWAS significant results. This study suggests several ways to consider PrediXcan's performance that will be of value to eQTL and complex human disease research.

How powerful are summary-based methods for identifying expression-trait associations under different genetic architectures?

Transcriptome-wide association studies (TWAS) have recently been employed as an approach that can draw upon the advantages of genome-wide association studies (GWAS) and gene expression studies to identify genes associated with complex traits. Unlike standard GWAS, summary level data suffices for TWAS and offers improved statistical power. Two popular TWAS methods include either (a) imputing the cis genetic component of gene expression from smaller sized studies (using multi-SNP prediction or MP) into much larger effective sample sizes afforded by GWAS - TWAS-MP or (b) using summary-based Mendelian randomization - TWAS-SMR. Although these methods have been effective at detecting functional variants, it remains unclear how extensive variability in the genetic architecture of complex traits and diseases impacts TWAS results. Our goal was to investigate the different scenarios under which these methods yielded enough power to detect significant expression-trait associations. In this study, we conducted extensive simulations based on 6000 randomly chosen, unrelated Caucasian males from Geisinger's MyCode population to compare the power to detect cis expression-trait associations (within 500 kb of a gene) using the above-described approaches. To test TWAS across varying genetic backgrounds we simulated gene expression and phenotype using different quantitative trait loci per gene and cis-expression /trait heritability under genetic models that differentiate the effect of causality from that of pleiotropy. For each gene, on a training set ranging from 100 to 1000 individuals, we either (a) estimated regression coefficients with gene expression as the response using five different methods: LASSO, elastic net, Bayesian LASSO, Bayesian spike-slab, and Bayesian ridge regression or (b) performed eQTL analysis. We then sampled with replacement 50,000, 150,000, and 300,000 individuals respectively from the testing set of the remaining 5000 individuals and conducted GWAS on each set. Subsequently, we integrated the GWAS summary statistics derived from the testing set with the weights (or eQTLs) derived from the training set to identify expression-trait associations using (a) TWAS-MP (b) TWAS-SMR (c) eQTL-based GWAS, or (d) standalone GWAS. Finally, we examined the power to detect functionally relevant genes using the different approaches under the considered simulation scenarios. In general, we observed great similarities among TWAS-MP methods although the Bayesian methods resulted in improved power in comparison to LASSO and elastic net as the trait architecture grew more complex while training sample sizes and expression heritability remained small. Finally, we observed high power under causality but very low to moderate power under pleiotropy.

Principles and methods of in-silico prioritization of non-coding regulatory variants

Over a decade of genome-wide association, studies have made great strides toward the detection of genes and genetic mechanisms underlying complex traits. However, the majority of associated loci reside in non-coding regions that are functionally uncharacterized in general. Now, the availability of large-scale tissue and cell type-specific transcriptome and epigenome data enables us to elucidate how non-coding genetic variants can affect gene expressions and are associated with phenotypic changes. Here, we provide an overview of this emerging field in human genomics, summarizing available data resources and state-of-the-art analytic methods to facilitate in-silico prioritization of non-coding regulatory mutations. We also highlight the limitations of current approaches and discuss the direction of much-needed future research.

A genomic perspective on HLA evolution

Several decades of research have convincingly shown that classical human leukocyte antigen (HLA) loci bear signatures of natural selection. Despite this conclusion, many questions remain regarding the type of selective regime acting on these loci, the time frame at which selection acts, and the functional connections between genetic variability and natural selection. In this review, we argue that genomic datasets, in particular those generated by next-generation sequencing (NGS) at the population scale, are transforming our understanding of HLA evolution. We show that genomewide data can be used to perform robust and powerful tests for selection, capable of identifying both positive and balancing selection at HLA genes. Importantly, these tests have shown that natural selection can be identified at both recent and ancient timescales. We discuss how findings from genomewide association studies impact the evolutionary study of HLA genes, and how genomic data can be used to survey adaptive change involving interaction at multiple loci. We discuss the methodological developments which are necessary to correctly interpret genomic analyses involving the HLA region. These developments include adapting the NGS analysis framework so as to deal with the highly polymorphic HLA data, as well as developing tools and theory to search for signatures of selection, quantify differentiation, and measure admixture within the HLA region. Finally, we show that high throughput analysis of molecular phenotypes for HLA genes-namely transcription levels-is now a feasible approach and can add another dimension to the study of genetic variation.

PRECISION MEDICINE: FROM DIPLOTYPES TO DISPARITIES TOWARDS IMPROVED HEALTH AND THERAPIES

Precision medicine research efforts both in basic science discovery and clinical implementation are well underway and promise to provide individualized preventions and treatments, improving overall health care delivery. To achieve these goals, advances in data capture and analysis are needed spanning different types of 'omic and clinical data. The efforts to enhance precise treatments for all may accentuate healthcare disparities unless specific challenges are identified and addressed. This session of the 2018 Pacific Symposium on Biocomputing presents the latest developments in this transdisciplinary research space of genomics, medicine, and population health.

Resolving TYK2 locus genotype-to-phenotype differences in autoimmunity

Thousands of genetic variants have been identified, which contribute to the development of complex diseases, but determining how to elucidate their biological consequences for translation into clinical benefit is challenging. Conflicting evidence regarding the functional impact of genetic variants in the tyrosine kinase 2 (TYK2) gene, which is differentially associated with common autoimmune diseases, currently obscures the potential of TYK2 as a therapeutic target. We aimed to resolve this conflict by performing genetic meta-analysis across disorders; subsequent molecular, cellular, in vivo, and structural functional follow-up; and epidemiological studies. Our data revealed a protective homozygous effect that defined a signaling optimum between autoimmunity and immunodeficiency and identified TYK2 as a potential drug target for certain common autoimmune disorders.

Epigenetic and Transcriptional Variability Shape Phenotypic Plasticity

Epigenetic and transcriptional variability contribute to the vast diversity of cellular and organismal phenotypes and are key in human health and disease. In this review, we describe different types, sources, and determinants of epigenetic and transcriptional variability, enabling cells and organisms to adapt and evolve to a changing environment. We highlight the latest research and hypotheses on how chromatin structure and the epigenome influence gene expression variability. Further, we provide an overview of challenges in the analysis of biological variability. An improved understanding of the molecular mechanisms underlying epigenetic and transcriptional variability, at both the intra- and inter-individual level, provides great opportunity for disease prevention, better therapeutic approaches, and personalized medicine.