Transcriptome genetics using second generation sequencing in a Caucasian population

Abstract

Gene expression is an important phenotype that informs about genetic and environmental effects on cellular state. Many studies have previously identified genetic variants for gene expression phenotypes using custom and commercially available microarrays. Second generation sequencing technologies are now providing unprecedented access to the fine structure of the transcriptome. We have sequenced the mRNA fraction of the transcriptome in 60 extended HapMap individuals of European descent and have combined these data with genetic variants from the HapMap3 project. We have quantified exon abundance based on read depth and have also developed methods to quantify whole transcript abundance. We have found that approximately 10 million reads of sequencing can provide access to the same dynamic range as arrays with better quantification of alternative and highly abundant transcripts. Correlation with SNPs (small nucleotide polymorphisms) leads to a larger discovery of eQTLs (expression quantitative trait loci) than with arrays. We also detect a substantial number of variants that influence the structure of mature transcripts indicating variants responsible for alternative splicing. Finally, measures of allele-specific expression allowed the identification of rare eQTLs and allelic differences in transcript structure. This analysis shows that high throughput sequencing technologies reveal new properties of genetic effects on the transcriptome and allow the exploration of genetic effects in cellular processes.

Figure 1. Array association p-values for RNA-Seq significant eQTLs

The p-value distribution using array data for RNA-Seq eQTLs significant at the 0.01 and 0.001 permutation threshold from both the exon and transcript quantification data is plotted. In each of the plots, the significant tail of the p-value distribution is substantially enriched indicating that eQTLs discovered through transcriptome sequencing are also significant in arrays. For each plot this excess is quantified using the q-value statistic 1-pi0 to estimate the proportion of true positives. Enrichment in the p-value distribution is higher for eQTLs discovered via transcript quantification and for eQTLs that are more significant at a more stringent permutation threshold.

Figure 2. Exon eQTLs by exon relative location

We investigated the proportion of discovered exon eQTLs relative to exon location in multi-exonic genes (normalized by number of exons tested within class). For 0.001 eQTLs we see that the proportion of discoveries is increased relative to middle exons in both first, second and last exons. We also see that we make proportionally more discoveries than any other class in the last exon of the gene.

Figure 3. Haplotype homozygosity for shared ASE haplotypes versus shared and unshared ASE haplotypes

(A) We assessed the degree to which shared ASE indicated a rare regulatory haplotype. We selected heterozygotes which were present in 6 or more individuals and assessed haplotype homozygosity between haplotypes which shared a significant ASE effect (p<0.05) (labelled as ASE in plot) versus those which shared and did not share an ASE effect (labelled as Control in plot) for all exons for which we did not have evidence for an eQTL (not significant at 0.05 permutation threshold). We see that when comparing the among significant (ASE) extent of haplotype homozygosity to that of significant vs. non-significant haplotypes (Control) where 2 or 3 individuals share the ASE significant signal there is greater haplotype homozygosity for the haplotypes which share the ASE signal indicating that these are on a more recent and rarer haplotype. This signal decreases when the ASE signal is shared in 4 or more individuals. Here the derived allele was selected as the one with the longest haplotype homozygosity without reference to the ASE signal. (B) For each heterozygote we plotted the extent of haplotype homozygosity for significant ASE haplotypes (X-axis) versus significant by not significant ASE haplotypes (Y-axis). We observed that the length of homozygosity is greater in the significant haplotypes compared to each other than when compared to non-significant corresponding haplotypes. Here the derived allele was selected as the one with the longest haplotype homozygosity without reference to the ASE signal.

Figure 4. Allelic alternative splicing effects

The two panels show examples of alternative use of exons centered on two significant ASE SNPs. In panel (A) there is a very significant signal of alternative use of exons (KS, P<10⁻¹⁴) between alleles where larger abundance coincides with larger diversity in transcript structure. In panel (B) there is significant ASE signal but not significant diversity in transcript structure between alleles.