Technical Comments

Comment on “Widespread RNA and DNA Sequence Differences in the Human Transcriptome”

Science  16 Mar 2012:
Vol. 335, Issue 6074, pp. 1302
DOI: 10.1126/science.1210484

Abstract

Li et al. (Research Articles, 1 July 2011, p. 53; published online 19 May 2011) reported more than 10,000 mismatches between messenger RNA and DNA sequences from the same individuals, which they attributed to previously unrecognized mechanisms of gene regulation. We found that at least 88% of these sequence mismatches can likely be explained by technical artifacts such as errors in mapping sequencing reads to a reference genome, sequencing errors, and genetic variation.

Li et al. (1) sequenced cDNA from lymphoblastoid cell lines derived from 27 individuals whose genomes have been sequenced at low coverage (2) and identified 10,210 sites of mismatches between an individual’s mRNA and DNA sequences [RNA-DNA differences (RDDs)]. RDD sites included all possible combinations of sequence mismatches, and the authors validated a subset of these mismatches by additional assays. These observations were interpreted as evidence for novel mechanisms of gene regulation, analogous perhaps to A→I RNA editing (3).

An alternative explanation is that some RDD sites are technical artifacts due to errors in mapping sequencing reads to a reference genome or systematic sequencing errors. To evaluate this possibility, we examined the sequence alignments used to call RDD sites [see supporting online material (SOM)]. Visualizing these alignments revealed a number of anomalies. For example, at the RDD site presented in Fig. 1A, all mismatches to the genome occur at the last base of reads aligned to the negative DNA strand. No such anomalies are seen in alignments around a positive control site (Fig. 1B). The biases in the first example are consistent with several known issues that cause spurious differences between Illumina sequencing reads and a reference genome; these include read-mapping errors between paralogous genomic regions and around insertions and deletions (2, 4), as well as position and strand biases in the error rate of Illumina sequencing (57).

Fig. 1

Identifying false-positive RDD calls. (A) RNA-seq read alignments around an RDD call from Li et al. (1). Plotted are the positions of read alignments to the genome surrounding the RDD site at chromosome 11, position 105,473,792. The solid lines show sequencing reads aligning to the (+) strand of the genome, and dotted lines are alignments to the (–) strand of the genome. At the center of the plot is the base corresponding to the RDD site; the reference base is in black, the nonreference base is in red, and both are labeled with respect to the (+) DNA strand. Alignments have been organized such that the mismatches to the genome are at the bottom of the figure. For plotting, we randomly sampled 20 alignments that match the genome at the RDD site; all 11 alignments that mismatch the genome are shown. (B) Read alignments around a positive control RDD site. Plotted are the positions of read alignments to the genome surrounding the known A→I editing site in AZIN1 (12) (on the forward strand, this site appears as T→C). The format is the same as in (A). For plotting, we randomly sampled 15 alignments that match the genome at the RDD site and 15 alignments that do not match the genome at the site. (C) Position biases in alignments around RDD sites. For each RDD site with at least five reads mismatching the genome, we calculated the fraction of reads with the mismatch (or the match) at each position in the alignment of the RNA-seq read to the genome (on the + DNA strand). Plotted is the average of this fraction across all sites, separately for the alignments that match and mismatch the genome. (D) Histogram of P values for the position bias test. For each RDD site with at least five reads mismatching the genome, we calculated a P value for the position bias test (SOM). Plotted is the histogram of these P values. If these sites were not consistently biased, the distribution of P values would be uniform; this is indicated with the dashed gray line.

We asked whether the patterns seen in Fig. 1A are typical among RDD sites. Indeed, mismatches to the genome at RDD sites are dramatically enriched at the ends of RNA sequencing reads; this contrasts with reads that match the genome at these sites (Fig. 1C). This pattern is evidence that many of the RDD sites are false positives due to mapping or sequencing errors.

To quantify what fraction of RDD sites may be false positives, we used metrics developed for calling single-nucleotide polymorphisms (SNPs) from Illumina sequencing data. In this context, it is known that a search for mismatches between aligned reads and a genome will result in large numbers of false-positive SNPs, many of which can be filtered out based on various criteria (2, 4, 8, 9). We used two criteria based on comparing, at each RDD site, the alignments of RNA sequencing reads that match the genome with the alignments of reads that mismatch the genome—a test for position bias and a test for strand bias (SOM). These tests provide quantitative measures for the intuition that there should be no systematic differences in strand or start position between alignments of reads covering the two alternative genotypes at a site and are similar to tests implemented in SNP-calling packages (4, 9).

The histogram of P values for the position bias test for the 7812 RDD sites with at least five reads supporting both bases reveals a clear skew toward low P values, indicating pervasive technical artifacts (Fig. 1D). At a P-value threshold of 0.01, 87% of these RDD sites fail either the strand bias test or the position bias test (at a P-value threshold of 0.05, the corresponding number is 93%). To test the specificity of these filters, we compared the reported RDD sites to a database of known A→I RNA editing sites (10). There are 23 sites in common between the two data sets; of these, 21 (91%) pass both of the filters. This indicates that we are largely only removing false positives.

Genetic variation is another source of false positives; an additional 1% of the putative RDD sites appear instead to be known genetic variants in these individuals (SOM). In total, we estimate that at least 88% (at a P-value threshold of 0.01) to 94% (at a P-value threshold of 0.05) of the RDD sites are likely false positives. This is probably an underestimate of the true false-positive rate, because some false-positive sites will pass the bias tests by chance and there are additional, unannotated SNPs in the genome.

Given the above results, we reexamined the validation experiments done by Li et al.. (1). These experiments are of two types. First, at 11 sites, the authors confirmed that the RDD event was absent from genomic DNA but present in cDNA by Sanger sequencing. At 6 of these 11 sites, the event is of the type A→G, and 4 of these 6 are present in a database of known A→I RNA editing sites (10); these are likely true positives. Of the remaining five sites, three fall in a single gene (HLA-DQB2) that is copy number variable in these individuals (11), and one (in the gene DPP7) overlaps a known SNP (at which the reported RDD type matches the known alleles) (2). We suggest that the authors have detected genetic variation rather than RNA-DNA differences at these sites. In sum, these experiments identify two previously unknown sites of A→I RNA editing and provide evidence for a single G→A event.

The second validation experiment involved identifying peptides corresponding to RDD events. In their table 3, Li et al.. (1) provide 17 examples where both the “DNA form” (the unaltered version) and the “RNA form” (the modified version) of peptides were detected by mass spectrometry. All but one of these sites fail the bias tests described above. We propose that the “RNA forms” of these peptides are in most cases normal forms produced by paralogous genes. Indeed, examination of the “RNA forms” revealed that seven match both the reported protein and additional proteins equally well, and four of the remaining ten match other proteins (in addition to the reported protein) with a single additional mismatch (Table 1 and SOM). It cannot be ruled out that the “RNA forms” of these proteins are instead normal forms caused by genetic variation in their paralogs. An additional possibility is that some “RNA forms” result from sequencing errors in the peptides.

Table 1

Characteristics of RDD sites reported in peptides. We reevaluated the peptides presented in table 3 of Li et al. (1). Repeated from that table are the gene names, positions and types of RDD sites, and “RNA forms” of protein sequences. We additionally show the numbers of aligned reads that mismatch the genome at each site, and the P values from the tests for position bias and strand bias at each site. P values in bold are less than 0.01. We used BLAST to search the human genome for matches to the peptides; given are the names of additional genes [apart from the one reported by Li et al. (1)] that match the peptide equally well (because these are the “RNA forms” of the peptides, the best matches have a single mismatching amino acid) and the number of genes with one additional mismatch (for a total of two mismatches) to the peptide. Mismatches are defined as either a substitution or an insertion/deletion of a single amino acid (13).

View this table:

In summary, we estimate that a minimum of 88 to 94% of the RDD sites identified by Li et al.. (1) are false positives due to mapping errors, sequencing errors, and genetic variation. It is possible that the remainder of RDD sites contain examples of novel mechanisms of gene regulation.

Supporting Online Material

www.sciencemag.org/cgi/content/full/335/6074/1302-d/DC1

Materials and Methods

Figs. S1 and S2

References

References and Notes

  1. Single-letter abbreviations for the amino acid residues are as follows: A, Ala; C, Cys; D, Asp; E, Glu; F, Phe; G, Gly; H, His; I, Ile; K, Lys; L, Leu; M, Met; N, Asn; P, Pro; Q, Gln; R, Arg; S, Ser; T, Thr; V, Val; W, Trp; and Y, Tyr.
  2. Acknowledgments: J. K. Pritchard is on the scientific advisory board and consults for the genotyping company 23AndMe. This work was supported by NIH grants MH084703 to J. K. Pritchard and GM077959 to Y.G.
View Abstract

Subjects

Navigate This Article