Blog

Figure 1. The RNA-seq data from the 1000 Genomes Project was processed through two parallel pipelines to generate imputed genotypes and transcript quantification for the eQTL analysis. STAR was used to align the sequencing reads to the reference genome followed by genotype imputation through the Gencove pipeline, producing the final imputed VCF files. The Salmon quasi-mapping and transcript quantification pipeline was used to process the transcripts and generate TPM values. The imputed genotypes and transcript quantifications were then integrated for downstream eQTL analysis. 

We processed RNA-seq data from 100 randomly selected individuals in the 1000 Genomes Project, consisting of 80 European and 20 African samples. These samples were processed through our standard leave-one-out imputation pipeline, which is typically used for low-pass whole genome sequencing data. Due to the unique characteristics of RNA-seq data such as non-uniform coverage and variability in transcript abundance, we changed our default alignment software from BWA to STAR ³, a specialized RNA-seq aligner, to handle RNA-seq spliced reads more effectively. Our imputation results shows an impressive overall concordance rate of > 99% and non-reference concordance (NRC) of > 82% sites passing filter between the imputed genotypes from RNA-seq data and truth genotypes obtained from whole genome sequencing of the same subjects, indicating that our standard imputation pipeline performs reliably even with RNA-seq data. Despite the uneven coverage in RNA-seq data, we found that the relationship between the effective coverage of RNA-seq samples and the imputation performance was only slightly worse than standard whole-genome low pass WGS data ⁴, suggesting that the high coverage in transcribed regions provided sufficient information for accurate genotype imputation.

Figure 2. Box plots showing the overall and non-reference concordance of imputed genotypes compared to true WGS genotypes. The panel on the left represents overall concordance. The second panel on the right represents non-reference concordance (NRC). Each data point corresponds to one of the samples. Blue box plot shows all sites. Yellow box plot shows PASS sites. 

Figure 3. Volcano plot showing the strength and direction of gene expression changes associated with the top eQTLs. Fold change shown on the x-axis and -log10 p-values on the y-axis. Fold change represents the gene expression difference between homozygous reference and homozygous alternate individuals. Upregulated genes identified by eQTLs are shown in blue, while downregulated genes are shown in red. The top 10 significant signals are highlighted with gene names, including the top 5 upregulated and top 5 downregulated genes. The analysis was performed using eQTLs imputed using the Gencove pipeline. 

We further validated these findings by comparing the gene expression changes influenced by these eQTLs with data from the GTEx database ⁹, confirming high concordance in the direction of effect between the two datasets, further reinforcing the reliability of our imputed genotypes in quantifying biologically meaningful eQTL associations. Below, we highlight two of these examples. 

For example, rs6674909, associated with the gene MTCH2, shows distinct expression levels depending on the genotype, with the alternate allele significantly increasing MTCH2 expression. MTCH2 plays a role in mitochondrial function and apoptosis, and its dysregulation has been linked to metabolic disorders such as obesity.

Similarly, rs12936231, an eQTL associated with the gene ORMDL3, shows an opposite effect, where individuals with the homozygous reference genotype have significantly higher expression of this gene compared to those with the homozygous alternate genotype. Changes in the expression levels of ORMDL3 have been implicated in various inflammatory and immune responses such as asthma and ulcerative colitis. These results demonstrate how these genetic variants, acting as eQTLs, influence gene expression and reinforce both the biological significance of these findings and the accuracy of our imputed genotypes for carrying out these analyses.

Figure 4. The relationship between eQTLs and gene expression for the MTCH2 and ORMDL3 genes. The left two panels display the expression of MTCH2 associated with rs66749409, with the top panel showing the boxplot from our experiment and the bottom panel presenting data from the GTEx database. The right two panels show the expression of ORMDL3 associated with rs12936231. This comparison demonstrates the high consistency of eQTL signals between our imputed genotypes and independent GTEx data. 

A cost effective solution to obtain accurate genotype calls from RNA-seq data.

Despite the uneven genome coverage of RNA-seq data, we found that it can be effectively leveraged by our imputation pipeline for accurate genotype imputation. With approximately half of the genome being intronic, the low but still detectable RNA-seq coverage in these regions provides additional value information for genotype imputation. This off-target sequencing, combined with the high coverage in transcribed regions, provides a comprehensive foundation for inferring haplotype blocks for genotype imputation. This approach provides a cost-effective and efficient alternative to additional sequencing or genotyping, enabling investigators to broaden the potential of RNA-seq data for more integrative downstream analyses that leverage both transcriptomic and genotypic data. Our experiment serves as a clear example of the feasibility of this approach by successfully conducting eQTL analysis on this dataset.

In addition to integrating the STAR RNA-seq aligner into our pipeline to ensure optimal alignment accuracy, we also tested the performance of RNA-seq imputation with our standard pipeline which utilizes BWA for alignment. Unlike STAR, BWA is not splice-aware, which reduces its ability to handle spliced reads effectively. This resulted in a somewhat reduced NRC of 0.78, a 5% drop in concordance, which highlights the importance of using the right tools for RNA-seq data to achieve optimal imputation accuracy.

Many organizations and research groups already possess vast RNA-seq datasets from numerous samples. By processing these existing RNA-seq datasets through our imputation pipeline, researchers can easily complement their transcriptomic data with high-quality genome-wide genotypes without additional sequencing effort and expand the scope for downstream analyses of RNA-seq data as demonstrated in this report.

As we continue to add new features to our pipeline, we expect further improvements in both accuracy and efficiency, making RNA-seq genotype imputation a highly feasible and valuable tool for large-scale research initiatives as well as clinical and drug development applications.