Screening the human exome: a comparison of whole genome and whole transcriptome sequencing

The study of common human diseases is rapidly moving away from an exclusive focus on common variants using genome-wide association studies and toward sequencing approaches that represent most variants, including those that are rare in the general population.

Although rapidly falling, the per base costs of next generation sequencing platforms still preclude the generation of large sample sizes of entirely sequenced genomes at high coverage. In addition to this economic constraint, it is widely appreciated that the very large number of variants identified in such studies will make it difficult to use association evidence alone to identify causal sites. For these reasons, there has been considerable interest in focusing attention on coding variants as a first step at complete representation of human variation. Part of the motivation for this approach stems from the experience with Mendelian diseases, in which 59% of the causal variants are either missense or nonsense mutations [1]. Although there has been considerable speculation on the topic, there are in fact no solid data showing that the picture is any different for common diseases, which may also be influenced by variants that are in or near protein coding sequence [1].

The most comprehensive approach for focusing on exons alone is clearly exome capture, where regions matching a defined set of coding exons are pulled from the genomic DNA (gDNA) using microarrays and then sequenced. However, this approach requires an initial and costly hybridization step. The cost of exome sequencing has contributed to the interest in sequencing the transcriptome (RNA-Seq) as an alternative, and possibly easier and less expensive strategy [2]. While this approach will clearly miss poorly expressed genes in whatever tissue is being studied, it does have the advantage of generating additional information, such as gene expression level and splicing patterns.

Although exome capture was demonstrated to identify approximately 95% of genomic single nucleotide variants (SNVs) in curated and non-paralogous exons [3], it is not currently known to what extent SNVs identified by RNA-Seq capture the full set of exonic SNVs identified by genomic sequencing. If the ability to capture SNVs by RNA-Seq is highly dependent on expression level, then this method would be useful only when performed in the appropriate tissue type. If, on the other hand, RNA-Seq at high coverage allows SNVs to be captured even in genes that are not highly expressed, then both methods could be useful for opening up sequencing studies to larger datasets in more diverse scientific studies.

Here, we have sequenced the entire genome and transcriptome of a single individual to high coverage. By comparing the SNVs identified in the transcriptome at different levels of coverage to those identified in the gDNA, we are able to directly evaluate how well RNA-Seq captures genomic variants.

If one simply considers all coding SNVs, our study suggests that about 40% can be identified by RNA-Seq using PBMCs as the RNA source. If we focus, however, on only PBMC-expressed genes, we find that 81% of coding variants are identified. This suggests that RNA-Seq may be a workable alternative for identifying exonic variants when performed in the appropriate tissue for the trait of interest.

One limiting factor in variant identification by RNA-Seq is the ability to uniquely align a given read. Although we were able to align 78% of our reads to a location in the transcriptome, only 69% aligned to exactly one location and could be kept in the analysis. Because exons are some of the most conserved regions of the genome, without the help of intervening and variable intron sequences it is much harder to align a read to the correct gene, and especially to have it align to only one location. This limits one's ability to identify variants in certain genomic locations, lowering the sensitivity of this method. Furthermore, if a read is uniquely aligned to the wrong gene, such as a paralog, then this can result in false positive SNVs being called in the cDNA. Restricting SNV calls to genes without paralogs did increase specificity from 0.47 to 0.54 and sensitivity from 0.39 to 0.42.

Another limiting factor is coverage. Because genes are expressed at different levels, in the random sampling of transcripts that are sequenced there will be gross imbalances. Some transcripts will have more than 1,000-fold coverage while other transcripts, although also expressed to some level in that tissue, will have coverage that is too low for variants to be accurately called. Given the diminishing returns of additional sequencing in terms of variants called (Figure 3), it simply will not be possible to use RNA-Seq to capture all exonic variants in genes with low expression levels. RNA-Seq is most useful for identifying variants in a tissue type that highly expresses genes related to the trait under study. One interesting possibility to improve the proportion of SNVs that can be called would be to use more than one tissue type as a source for the RNA. For example, data available on expression of core exons in muscle shows that adding cDNA from this tissue to PBMC cDNA would increase the number of adequately expressed transcripts by 68% [6, 7]. There are no expression data for a more easily accessible tissue such as skin for this exon array publicly available; however, one can extrapolate that adding cDNA from almost any tissue would be similarly beneficial to analysis.

A large number of false positive SNVs were identified in this dataset: even when restricted to PBMC-expressed genes without paralogs the specificity was only 0.72. Many of these false positives, however, were due to the very high coverage produced in this study. At low levels of coverage, reads that could produce false positive calls are not yet abundant enough to pass the quality control filters used, and the specificity remains high. However, as more reads are added to the dataset, the number of incorrect alignments and sequencing mistakes increases, pushing more and more of these false positive calls over the quality control filters. Thus, the more coverage that is added in the search for true positives, the more false positives that will appear in the data. We found that the specificity for the dataset as a whole could be increased from 0.47 to 0.67 (and from 0.72 to 0.9 for PBMC-expressed genes without paralogs), at a substantial cost to sensitivity, by restricting the permissible read depth range for SNV calls. The specificity was low at low read depths, as is expected when a call is supported by less data; interestingly, the minimum read depth required for high specificity increased as coverage increased, supporting the view that high coverage introduces more noise and comes with a requirement for stricter quality control. An additional support to this view was the finding of low specificity at very high read depths. Both a minimum and a maximum read depth cutoff are advisable for increasing the specificity.

Some of the false positives found in this dataset are certainly due to the incorrect alignment of reads, as even when a gene has no annotated paralog it can have sections with sequence similarities to other genes that permit errors during alignment. Another alignment problem that is unique to RNA-Seq is incorrect alignment at the very ends of a read due to splicing. Because reads that span exons require a certain number (four in our study) of bases to fall on both sides of the exon-exon boundary for proper alignment, reads that cross an exon boundary at the very edge of the read will have between one and three bases aligned to the wrong location. Our study removed most of this type of false positive SNV during the quality control process. Another scenario that would produce a seemingly false positive SNV would be sufficient coverage in the cDNA to call the variant but insufficient coverage in the gDNA sequencing to do the same: however, in our study such bases had a median coverage of 27× in the gDNA, which does not support this explanation. Also, previous studies have shown that variants can be present in the RNA but not the gDNA due to RNA editing [2, 8].

There are also likely to be disagreements between gDNA and cDNA sequencing stemming from expression differences. For example, some SNVs found in the genomic sequence may be missed in the cDNA due to expression balances, when an individual is heterozygous for a given SNV yet the reference allele is much more highly expressed. Also, an SNV may be called heterozygous in the gDNA but homozygous in the cDNA due to expression imbalances. In our dataset, 10% of the 19,054 SNVs that overlapped by location between gDNA and cDNA were heterozygous in gDNA and yet homozygous in cDNA. Differences in zygosity between the two methods can also result from insufficient coverage in one dataset or the other: for example, 1% of these 19,054 SNVs were homozygous in gDNA and yet heterozygous in cDNA, and the median coverage for these SNVs in gDNA was only 12×, compared to 28× for SNVs that matched for zygosity. It is likely that many of the discrepancies where the SNV was heterozygous in gDNA and homozygous in cDNA also resulted from low coverage, as the median coverage for this group was only 8× in the cDNA. It should also be noted that only 13 of the 19,054 SNVs that overlapped by location had completely mismatched alleles, such as the cDNA being homozygous for a G and the gDNA homozygous for C when the reference allele was A.

Our study found that although the percent overlap of our false positive SNVs with dbSNP entries was far less than that of true positive (94%), false negative (89%), or all exonic gDNA (90%) SNVs, it was still substantially greater than zero (23%). Furthermore, we found that the percent of false positives corresponding to dbSNP entries increased as the PBMC expression level increased (Supplemental figure S4 in Additional file 1), implying that the 'false positives' seen at high expression levels may actually be true positives that were simply not seen in the gDNA due to issues with coverage or alignment or because of RNA editing. Additionally, we found that only 13% of the false positive SNVs found in dbSNP had been genotyped in HapMap, compared with 69% of the true positives found in dbSNP. This suggests that many of the dbSNP entries matching false positives are less well-curated, with less supporting evidence. A brief inspection of a subset of the SNVs showed that false positives were more likely to have cDNA evidence (as opposed to gDNA evidence) supporting their dbSNP entry than were true positives; this could be a reflection of either misalignment of the cDNA reads for dbSNP entries, or RNA edits that would not be seen in the gDNA. Finally, we found that the percent overlap of false positives with dbSNP decreased as the number of lanes of sequence data increased (Supplemental figure S3 in Additional file 1). This finding corresponds with the fact that the overall specificity decreased as coverage increased (Figure 1); these phenomena are likely caused by the increasing ability of the random noise inherently present in the data to overcome quality control cutoffs as more and more reads were added, as described above.

Two previous studies have also looked at using RNA-Seq to identify SNVs. Chepelev et al. [2] used a dataset of 27 million uniquely aligned 30-bp RNA-Seq reads; while they detected 50% of known exons with at least 1× coverage, and identified approximately 11,000 SNVs, they did not compare these data to gDNA sequence and thus could not calculate sensitivity or specificity of their methods. Shah et al. [8] used a dataset of 183 million RNA-Seq reads of 38.9 bp each, 55 million of which aligned to exons or exon junctions. They compared these data to 2.5 billion aligned gDNA reads of 48.2 bp to better understand the evolution of mutation in a lobular breast tumor. Although they showed that the number of SNVs called through RNA-Seq increased as the number of reads increased, they did not discuss the sensitivity or specificity of SNV calls when compared to whole genome sequencing, nor did they analyze how the number of SNVs changed as known expression level changed.

Another technology that has been used to sequence coding variants is exome capture. Because this method sequences reads that are from gDNA but enriched for the portions of interest, there are no complications with aligning splice junctions or being limited by expression level. A recent study showed that when restricted to the non-paralogous exons of 16,496 curated protein coding genes, exome capture utilizing 41 million 76-bp reads (one quarter the number of bases we aligned) captured SNVs with a sensitivity of 95% and a specificity of 90% [3]. Although this performance outshines the RNA-Seq data presented here, RNA-Seq does have some advantages beyond the ability to inexpensively genotype exonic SNVs. It can also be used to identify expression differences between individuals or even between alleles within an individual, which can lead to discovery of a nearby causal variant even if it is not exonic. RNA-Seq may also provide insight into novel exons, splice junctions or splice forms in the tissue or cell type being studied that might not be recognized as protein coding in genomic sequencing, or captured with targeted exomic sequencing.

While it is most useful to perform RNA-Seq in a tissue relevant to the trait under study, other tissues can also be of some use. For example, data from Heinzen et al. [6] revealed an r² of 0.23 between transcript expression in PBMCs and in brain, and that 63% of the transcripts highly expressed in brain were also expressed in blood to a level that allowed for consistent SNV detection by RNA-Seq (defined as expression level of at least 4% of the most highly expressed transcript in both).

Here we show that RNA-Seq captured 81% of the exonic variants from genes that were well expressed in the source tissue. Although its usefulness is limited to these genes, the cheaper cost involved, as well as the extra information gained about expression and splice variants, may make this method a workable alternative to genomic sequencing or exome capture for groups that have access to the right types of tissue.