ASHG 2014 Poster - eQTL in Diploid Genomes

I presented this poster at the 2014 American Society of Human Genetics Meeting in San Diego.

Transcript

1. www.jax.org Allele Specific Expression and eQTL in Diploid Genomes Steven C. Munger1, Narayanan Raghupathy1, Kwangbom Choi1, Daniel M. Gatti1, Petr Simecek1, and Gary A. Churchill1 1The Jackson Laboratory, Bar Harbor, Maine 04609 USA INTRODUCTION The Diversity Outbred (DO) heterogeneous mouse stock is derived from eight inbred founder strains that together capture ~90% of the known genetic variation in the mouse (40+M SNPs, 6+M indels). The DO has been maintained for multiple generations by randomized outcrossing. As a result each chromosome is a genetically unique, balanced mixture of eight founder haplotypes, and each animal contains hundreds of recombination events. DO mice are heterozygous with respect to founder origin at 87.5% of loci and there are 36 possible genotypes (8 homozygous + 28 heterozygous). Over half of all 100bp RNA-seq reads will have at least one SNP segregating in the DO population. Genotyping arrays can distinguish among the DO founder genotypes (Welsh and McMillan 2012) and a hidden Markov model (HMM) has been developed to reconstruct the 36-state haplotype mosaic of each DO genome (Svenson et al. 2012). We have developed Seqnature, a tool that can construct individualized diploid genomes with founder strain-specific SNPs, insertions, and deletions (indels) from haplotype data imputed for each DO sample. The tool also builds a diploid gene annotation file with gene coordinates adjusted to account for indels. With these files, we construct a diploid set of isoform sequences (transcriptome) for each sample. RNA-seq reads are aligned with Bowtie (Langmead et al. 2009) to this transcriptome, and gene and isoform abundance is estimated with an EM algorithm (RSEM, Li and Dewey 2010). ABSTRACT The emergence of high throughput sequencing (HTS) technologies has coincided with the development of advanced genetic reference populations including the mouse Diversity Outbred (DO) heterogeneous stock. The application of HTS to genetically diverse mapping populations has the potential to provide nucleotide resolution of causal variants underlying phenotypic differences. However the increase in information content comes at the cost of increased analytical complexity. We have developed novel methods and software to exploit the high genetic diversity and heterogeneous diploid genomes of the DO to yield new layers of information and inform fine mapping of phenotypic and expression QTL (eQTL). Importantly, these methods can be readily applied to any human expression dataset where genotyping data is available. RNA-seq alignment to individualized diploid genomes yields direct, accurate estimates of allele specific expression (ASE), and improves eQTL identification and resolution. To illustrate this, we profiled the liver transcriptomes of 453 DO mice by RNA-seq and estimated gene, isoform, and allele expression. We correlated gene expression differences to genetic variation and identified 9,000 local expression QTL (eQTL) and 900 distant eQTL, showing that most of the variation in transcript abundance derives from segregating local genetic variation. Allelic expression differences confirmed that cis-acting mechanisms underlie most local eQTL, and DO allele estimates correlate well with gene expression in livers from the eight DO founder strains. Most eQTL appear biallelic suggestive of a single causal variant, however complex 3- and 4- allele patterns are observed. Cis- eQTL with allelic expression patterns that deviate from the known strain ancestry are most amenable to fine mapping, and in some cases we predict a single causal variant. Distant eQTL with large effects are rare in the adult liver transcriptome, however we have developed conditioning methods to amplify trans effects and test candidate regulatory genes in the interval. ACCURATE ESTIMATION OF ALLELE SPECIFIC EXPRESSION We simulated 100bp single end Illumina reads with the Flux Simulator (Griebel et al. 2012) from a genetically heterogeneous Diversity Outbred (DO) genome. We aligned reads with Bowtie (v0.12.8) to each of the NCBIM37 reference transcriptome and individualized transcriptome. Read alignment files were input into RSEM (v1.2.1) to estimate gene level abundance. Importantly, because we align to a diploid version of the transcriptome, our pipeline generates direct estimates of allelic abundance. Gene abundance estimates are improved by read alignment to individualized transcriptomes (Munger et al. 2014). Moreover, estimates of allele-level gene abundance are well-correlated to the simulated ground truth values (r = 0.82), with the estimated frequency of the lower-expressed allele differing by on average by less than 7% (median = 4%) from the ground truth value (A). Increasing the sequencing depth to 30 million reads and stringency to require a minimum of 10 unique allele alignments produces a corresponding increase in the accuracy of allelic gene-level abundance estimates (B). PREDICTION OF TRANS REGULATORS WITH CONDITIONAL SCANS Trans eQTL are indicative of a genetic interaction between the controlled transcript and an upstream regulatory gene. We have taken a conditional mapping strategy to predict the most likely candidate gene(s) underlying the trans eQTL, with the assumption that the expression level of the causal gene is responsible for the variance in the controlled transcript. For example, Lrtm1 has both a cis-eQTL and trans-eQTL (A). We regress out the local genotype to amplify the trans signal (blue line in B), and then scan Lrtm1 conditioning on the expression of each of the 331 genes within +/- 5Mb of the trans peak. Conditioning on the expression of Igsf23 almost completely abolishes the trans association for Lrtm1 (red line in B), and the founder effect plots suggest that Igsf23 is a repressor of Lrtm1 (C, D). ALLELIC IMBALANCE IS PERVASIVE Our analytical pipeline outputs direct estimates of allelic abundance that are tagged with the founder origin of the gene allele. We observe significant strain biases in allelic expression for most expressed genes in the liver (n=13,068/ 16,985 genes; FDR 1%) by comparing allele abundance estimates across 453 DO samples (n= 906 allele estimates). These biases indicate extensive local genetic variation in the DO population that act in cis to regulate allele-level gene abundance, and their 8-strain patterns provide high predictive power to identify the specific variant(s) underlying the effect. MOST LIVER eQTL ARE LOCAL Next we profiled liver RNA-seq in 453 DO samples, aligned all samples to their individualized diploid transcriptome, estimated gene and allele abundnance using our EMASE algorithm (Ragupathy et al., in preparation), and mapped loci that were significantly linked to variation in gene expression (eQTL). We find that eQTL mapping with gene abundance estimates derived from individualized alignments identifies nearly 3,000 novel local eQTL that are masked by the common reference alignment (Munger et al. 2014). The vast majority of significant eQTL (10,243/ 11,932 total, FDR 5%) map within 5Mb of the gene they control, as evidenced by the thick diagonal row of dots in the transcriptome map below. Position of the gene is plotted on the y-axis, and position of the significant eQTL is plotted on the x-axis. CONCLUSIONS •  Applying a common reference alignment strategy to genetically diverse RNA-seq reads can cause spurious eQTL associations (Munger et al. 2014). One erroneous alignment can produce two (or more) spurious eQTLs. Importantly, alignment to a common reference may obscure real biological signals (e.g. local eQTLs) underlying GWAS hits. •  Most eQTL in the adult liver are local and cis-acting. Because the DO is a highly recombinant multiparental outbred population with full genome sequences for the eight founder strains, it is particularly powerful for high resolution mapping of the causal genetic variant(s) underlying variance in quantitative endophenotypes such as transcript and protein abundance. •  This work was supported by The Jackson Laboratory, and National Institutes of Health (NIH) grants P50GM076468 to G.A.C, DK091207, DK66369, and DK58037 to A.D.A., and F32HD074299 to S.C.M. START: Diversity Outbred Mouse Seqnature HMM Abundance Estimates •  Gene •  Isoform •  Allele Merge Build Isoforms Variant Calls (Sanger vcf – Keane et al. 2012) Pseudophased 8-State Haplotypes R L Individualized Diploid Genome & Gene Set R L 36-State Genotypes SNP Genotyping Data (MUGA) RNA-seq Short Reads (fastq) T1L ATCGCGTTTACCGGA_B6! T1R ATCGCATTTACCGGA_WSB! ! T2L CGGAATCGCCTTAACGGT_CAST! T2R CG-AATCGCCTTGACGGT_NOD! ! T3L TTACCGG---CT_129S1! T3R TCACCGGTAGCT_PWK! ! Individualized Diploid Transcriptome EM Algorithm (RSEM – Li et al. 2010) Read Alignment (Bowtie – Trapnell et al. 2009) Individualized Alignment and EM Quantitation Pipeline Collaborative Cross Breeding Funnel Diversity Outbred Stock … Random Outbreeding A B C D E F G H PREDICTION OF CAUSAL VARIANTS UNDERLYING CIS-eQTL We have demonstrated that most eQTL in the adult liver transcriptome act in cis to affect transcript abundance of the linked founder allele. These founder allele patterns can be highly informative for fine-mapping down the locus. For example, Slc22a2 has a strong cis-eQTL after alignment to the individualized transcriptomes (A). The eQTL appears biallelic (B,C), exhibiting one of the 127 possible partitions of the eight founder strain alleles. Compared to an F2 intercross with three possible genotypes (e.g. AA, AB, BB), this complex strain pattern is only observed for 22 out of ~50k SNPs in the flanking region (D). steven.munger@jax.org Twitter: @stevemunger Seqnature can be downloaded at https://github.com/jaxcs/Seqnature Complex 3- and 4-allele patterns are observed for ~10% of cis-eQTL E F Common allele patterns – those that reflect the known strain ancestry (ie. wild-derived strain as outgroup) – resist fine-mapping by this approach (E). However, for a surprising number of genes, the eQTL allele pattern is rare and can be fine-mapped to a few candidates or even a single SNP/indel (F). 129, NOD, PWK Low 129, NOD, PWK High Lrtm1 Igsf23 Lrtm1 Conditional scan on Igsf23 expression Lrtm1 Main Scan Trans eQTL Cis eQTL 331 genes within +/- 5Mb A B C D