Intra-Exon Motif Correlations as a Proxy Measure for Mean Per-Tile Sequence Quality Data in RNA-Seq

Abstract

Given the wide variability in the quality of next-generation sequencing data submitted to public repositories, it is essential to identify methods that can perform quality control on these data sets when additional quality control data, such as mean tile data, are missing from public repositories. In this study, we present evidence that correlating counts of reads corresponding to pairs of motifs separated over specific distances on individual exons can be used as a proxy mean tile data in the data sets we analyzed and hence could be used when mean tile data are not available. As test data sets we use the Homo sapiens in vitro transcribed (IVT) data set, and a Drosophila melanogaster data set comprising wild and mutant types. We find that a FastQC analysis of the available parts of these data sets demonstrates that the per-tile sequencing quality is good for all the data sets apart from the mutant-type data where the mutant-r3 data are worse than the mutant-r2 data. Correspondingly, intra-exon motif correlations are reasonably large for all data sets except this latter case where the mutant-r2 correlations are low and the mutant-r3 correlations close to zero. We propose that these extremely low correlations are indicative of bias of technical origin, such as flowcell errors. In addition to this, the intra-exon motif correlations as a function of both guanosine-cytosine (GC) content parameters are somewhat higher and less dependent on the GC content parameters in the IVT-Plasmids messenger RNA (mRNA) selection free RNA-Seq sample (control) than in the other RNA-Seq samples that did undergo mRNA selection: both ribosomal depletion (IVT-Only) and PolyA selection (IVT-PolyA, wild type, and mutant).

1. INTRODUCTION

Next-generation sequencing (NGS) methods have revolutionized nucleic acid sequencing largely as a result of the employment of fluorescence-based nucleotide chemistry to generate a light signal on nucleotide incorporation (Baumeister et al, 1987; Connell et al, 1985; Soper, 1985), miniaturization, and massively parallel sequencing reactions (Mardis, 2006). Although these have, to a degree, simplified the core sequencing process allowing reactions to be performed in clusters to generate enough signal and in parallel to increase throughput, NGS technologies share the same complex preparatory procedures (Ji and Shendure, 2008).

These are typically fragmentation of fragments to the size appropriate for the target sequencing platform, amplification (e.g., PCR), and ligation of synthetic sequencing adapters for the sequencing platform. Such high-throughput sequencing technologies generate millions to billions of reads in a matter of days and generate large data sets (Metzker, 2010)—it has enabled a number of large-scale sequencing projects. These include the 100,000 genomes project in the United Kingdom (Siva, 2015; Gallagher, 2014), the National Institutes of Health Precision Medicine 1 million genomes project in the United States (Reuters, 2015), and a 1 million genome project by the Beijing Genome Institute in China (Cyranoski, 2016) to name just a few.

RNA-Seq, on which this work focuses, is a high-throughput NGS technique for estimating the contraction of messenger RNA (mRNA) transcripts in a transcriptome. It provides wider coverage of the transcriptome than microarrays as its methods involve the direct sequencing of transcripts of RNA found in the sample (Gerstein et al, 2009; Kukurba and Montgomery, 2015). RNA-Seq can be used to study various types of RNA present: total RNA, mRNA, pre-mRNA, and non-coding RNA (ncRNA), such as microRNA and long ncRNA enabling it to be used to study alternative splicing events (Dickerson et al, 2014; Park et al, 2013). Furthermore, RNA-Seq achieves this at a higher resolution (Kukurba and Montgomery, 2015) than other technologies. To quantify gene expression, a mapping process is combined with gene boundary information so as to count the number of transcripts that map to a given gene or exon region (Garber et al, 2011; Gerstein et al, 2009; McCue et al, 2008).

RNA-Seq has transformed our view of the extent and complexity of the transcriptome through deep sequencing (Gerstein et al, 2009), and also as a result of the increased precision, the technique offers over other methods. While recent developments in the RNA-Seq workflow, from sample preparation to sequencing, have furthered our understanding of the transcriptome, they have also required substantial effort for data analysis and computation, and given the complexity of RNA-Seq workflow necessitates study of the bias that can be introduced in the preparatory steps (Kapranov et al, 2011; Kukurba and Montgomery, 2015). Characterization of bias in RNA-Seq is especially incumbent, given that the method sequences and measures the transcriptome indirectly using reverse-transcribed complementary DNA (cDNA) (Bowers et al, 2009).

Bias introduced in the preparatory steps can have a profound effect on the raw data and typically manifest themselves as sequence-specific or positional biases, while bias introduced by the sequencing process itself is often systematic in nature (Meacham et al, 2011).

The main obstacle to obtaining accurate estimates of transcript expression from RNA-Seq data is nonuniformity in the distribution of mapped reads to the reference genome, which reduces the certainty that the measured counts of mapped reads reflect the true expression of the transcript within the cell's transcriptome. These bias have numerous sources such as wet-lab sample preparatory techniques, the sequencing process itself (Alnasir and Shanahan, 2015; Dohm et al, 2008), and the potential for errors in post-sequencing data processing. They perturb the uniformity of the distribution of mapped reads to a reference genome (Lahens et al, 2014), and such bias manifests itself as sequence-specific or positional (Donaghey et al, 2011). Also, positional biases can occur due to random hexamer priming in sample preparation (Hansen et al, 2010).

Large amounts of these raw RNA-Seq read data are deposited in public repositories, such as the Sequence Read Archive (SRA) (Leinonen et al, 2011) and Gene Expression Omnibus (GEO) (Edgar, 2002). Furthermore, the SRA, for example, does not require a quality check on submission (Nakazato et al, 2013) and has shown poor annotation of sequencing protocol steps—both at the top-level study and at individual experiment record level (Alnasir and Shanahan, 2015). Hence, it is critical that methods are developed to characterize and quantify bias in these data sets. Such methods can augment the analysis of guanosine-cytosine (GC) metadata in the data sets or can serve as an alternative measure when this metadata is not present (Andrews, 2010).

FastQC (Andrews, 2010) is often used to perform a number of quality control checks on the raw reads in high-throughput sequencing data. The per-tile sequence quality measure reported by FastQC is particularly important because it represents the deviation of each nucleotide from the average quality for each flowcell tile. The measure is used by FastQC to produce a heat map plot of this deviation—cold colors indicate base quality scores that are at, or above, the average quality scores for other tiles on the flowcell, and hotter colors indicate base quality scores that are below the average quality scores for other tiles on the flowcell.

Hence, a good heat map plot should be blue all over. Per-tile sequence quality can only be computed from Illumina sequencing data in which the original sequence identifiers are retained, that is, it is only present in raw Illumina reads. However, this information is not retained after processing Illumina reads into, for example, Binary Alignment Map (BAM) files after sequence alignment. This difficulty is compounded when using sequencing data sets from public repositories, such as the SRA, where not all submissions include the raw reads.

In this study, we propose that a method we have previously devised, which applies distributed-computing to quantify the sequence-specific deviations in the uniformity of mapped reads (Alnasir and Shanahan, 2017), can be used as a proxy measure when mean tile data are not available for short-read RNA-Seq data. Our method uses counts of reads overlapping motifs and works at the deep, read level. This approach is based on the assumption that 4-mers in short reads from one region on an exon will be correlated with 4-mers in short reads from another region of the same exon.

To provide the capacity to process the amounts of data typical in transcriptomic data sets (Campbell et al, 2015), our analysis employs parallel distributed computing algorithms and infrastructure using the Apache Spark platform—it is named Hercules and is available on GitHub (Alnasir, 2018) and at Zenodo (Alnasir and Shanahan, 2022; a message passing interface (MPI) version has now also been made available). We demonstrate this using a controlled in vitro transcribed (IVT) data set created by Lahens et al (2014) for the purpose of characterizing bias introduced in RNA-Seq library preparation, as well as a Drosophila melanogaster data set.

The first data set (IVT) was produced utilizing IVT in Escherichia coli to clone a pool of ∼1000 preselected human plasmids from the Mammalian Gene Collection (MGC) (Derge et al, 2009). Because the sequences and expression levels of these plasmids are known, and they do not undergo splicing, this allowed them to generate a highly controlled set of samples, and therefore a controlled data set, in which the source of biological variation in the samples is minimized. The samples in this data set were then subjected to different RNA-Seq preparatory protocols, specifically varying the step in which mRNA is selected.

This enables the study and quantification of the effect of these steps on coverage levels of the MGC transcripts when they were aligned to the human reference genome (hg19/grch37). They found that the mRNA selection methods employed in RNA-Seq protocols, poly-A and ribosomal depletion, both resulted in significant fold changes in the coverage of the IVT MGC plasmids when compared with sequencing the IVT MGC plasmids directly (without mRNA selection). Importantly, as the bias introduced in this data set is well characterized and attributed, we apply our analysis method to a selection of relevant samples.

In previous work, we applied our analysis method to replicates of two samples from a D. melanogaster data set produced from typical biological specimens using conventional RNA-Seq protocols (Alnasir and Shanahan, 2017). These are two small, but whole transcriptomes—those of the fruit fly species D. melanogaster wild type and mutant-r2, comprising ∼12.9 and 15.0 M reads, respectively. We will re-examine our analysis of these data with respect to sequencing tile means. While the IVT data set offers a set of samples that allow us to study the effect of different RNA selection methods, the D. melanogaster data sets have the same RNA selection method applied to all samples and vary only in the glass eye mutation—that is, the technical variation is fixed, and the biological variation should be minimal. Furthermore, the Drosophila species and its reference genome are extremely well studied and annotated, and the data have excellent provenance.

2. MATERIALS AND METHODS

2.1. Quantifying sequence-specific deviation in the distribution of mapped reads across exons

As explained in detail in Alnasir and Shanahan (2017), the uniformity of read distribution across an exon (Fig. 1) can be quantified by computing Pearson's correlations of the counts for the given motif pair in all exons within the data set by aggregating the motif pair counts at a given distance apart (motif-spacing) regardless of position within the exon (Fig. 2). We used motif-spacings of 10, 50, 100, and 200 base pairs (bp). An ideal data set would have perfect correlations for motif pairs (for instance, +1.0 for the Pearson correlation coefficient) for any given motif-pair and motif-spacing.

FIG. 1.

Typical distribution of RNA-Seq reads mapped to an exon (Alnasir and Shanahan, 2018).

FIG. 2.

Quantification of read coverage using short pairs of sequence motifs (4-mers) within the reads shown in light shading. Motif-pairs show variable correlations that are used to quantify sequence-specific deviations in the distribution of mapped reads across exons (Alnasir and Shanahan, 2017).

To thoroughly examine the affect of sequence-specific motifs on the uniformity of read distribution, we analyzed the correlation for all 4-mer motifs ranging from AAAA to GGGG (i.e., 4⁴ combinations) in the RNA-Seq reads of an aligned Sequence Alignment Map (SAM) file. We verified this method by running our analysis on an artificially created transcriptome with in-built uniform distribution of reads (see Section A, Supplementary Table S1, and Supplementary Fig. S1 in the Supplementary Data).

The effect of extremes of GC content in RNA-Seq data has been discussed in numerous studies (Chen et al, 2013; Harrison et al, 2010), and we therefore also investigate the effect the mean GC content of reads within the exon, and the GC content of the 4-mer motif itself, has on the distribution of reads across the exon. To partition reads by mean GC content (of the exon), we define binned GC content ranges: 30%–40%, 40%–50%, 50%–60%, and 60%–70%. Since exons of such extreme mean GC values are not observed, that is, <30% or >70%, at least not at sufficient levels, we are not examining these. Given we are working with 4-mers, the motif GC is fixed to the range 0%, 25%, 50%, 75%, and 100%.

2.2. Homo sapiens IVT RNA-Seq data set

We have analyzed three Homo sapiens samples, which were prepared by IVT RNA-Seq. Lahens et al (2014) in their study used RNA that has been IVT from cDNA clones in E. coli. Their rationale for the development of this data set was that the “nucleotide sequence at every base was known, the splicing pattern established, and the expression the level coverage is uniform across the transcript.” This means that any bias occurring in the coverage of reads in these three samples must be as a result of technical rather than biological origin. In addition, these samples are known to demonstrate intra-exon coverage bias (we have used their H. sapiens data; Lahens et al, 2014) and hence is ideal to study from this perspective.

We explored known-bias in this RNA-Seq, by analyzing intra-exon motif pair correlation within the reads, we performed analysis on three samples from the H. sapiens data set, which applied different library preparation protocols to each sample during RNA-Seq. These samples from Lahens et al (2014) are known to demonstrate intra-exon coverage bias (we have used their H. sapiens data). Although the raw data deposited in GEO for this data set have not been aligned, the library preparation protocol for each sample and the alignment and post-processing strategy applied have been clearly documented. We, therefore, aligned the reads of the samples to the reference hg19 genome according to the documented parameters to generate the SAM files for analysis by our method.

The RNA in these samples has been transcribed from cDNA clones in E. coli DH5α cells. The data set comprises a pool of 1062 RNAs from a full-length human cDNA library sequenced using RNA-Seq. The first sample, IVT-Only, had its IVT RNA subjected to ribosomal RNA depletion before sequencing, while the second sample, IVT-PolyAsel, had polyadenylated selection applied instead of ribosomal depletion—these are two different, routinely used protocols for selecting specifically mature (mRNA) from RNA samples. The third sample, IVT-Plasmids, is our control as it was produced by direct sequencing of the Human IVT plasmids without RNA selection (i.e., neither ribosomal depletion nor PolyA selection was applied). These data are deposited at the GEO database with the ID GSE50445 (Black et al, 2012).

The distribution of protein-coding exon lengths in H. sapiens is shown in Supplementary Figure S2 of the Supplementary Data. We note that, although the median exon length in H. sapiens is 121 bp (shorter than Drosophila) and ∼80% of the exons are <200 bp (Sakharkar et al, 2004), the remaining 20% of exons will contribute to motif pair correlations at 200 bp apart.

2.3. D. melanogaster RNA-Seq data set

To investigate intra-exon motif pair correlation within the reads from a “typical data set,” we have used two Drosophila (species D. melanogaster) transcriptomics data sets, which differ by mutation gl[60j] in the eye–antennal disk. These are the full transcriptomes of the wild-type and mutant glass eye mutations, acquired from Stein Aerts Laboratory of Computational Biology at the University of Leuven, Belgium. These data are deposited at the GEO database with the ID GSE39781 (Aerts, 2012). The D. melanogaster data sets featured in a research publication by Naval-Sánchez et al (2013).

Supplementary Figure S3 shows the distribution of exon lengths in the Drosophila genome—the median exon length is 298 bp. This is important because it shows that most of the exons are longer than 200 bp and therefore can have data to compute correlations.

The information regarding the source of the biological samples, the sample preparatory protocols, and post-sequencing processing that were applied to samples in these two data sets is documented in Table 1.

Table 1.
The Sample and Library Preparation Protocols, Together with the Data Processing Steps, Applied to the RNA-Seq Data Sets Used in This Analysis

Homo sapiens (IVT-Plasmids, IVT-Only, and IVT-PolyAsel) Drosophila melanogaster (wild type and mutant)

Replicates: 1×.sam file for each IVT sample (no replicates) Replicates: 4×.sam files (two replicates for each species)

Sample preparation Sample preparation

Glycerol stocks containing individual cDNAs (cloned into pCMV-Sport 6 plasmid) from the MGC (Derge et al, 2009) were produced. Plasmid DNA was extracted from these glycerol stocks and plated at 50 ng per well in 384-well plates. The contents of three 384-well plates (total of 1062 human transcripts) were collected. The plasmid library was then amplified by transferring 10 ng into Escherichia coli DH5α cells (cat. no. 18258012; Invitrogen, Life Technologies, Carlsbad, CA). The heat shock method was used to transform E. coli (see Lahens et al, 2014, for more details). The plasmids were then purified using Qiagen (Hilden, Germany) Maxiprep Kit (cat. no. 12163) according to the manufacturer's protocol. Samples were sequenced on the Illumina HiSeq 2000 platform For the wild-type D. melanogaster, fly stocks (Canton-S and strain RAL-208) were obtained from the inbred collection of T. Mackay. For the mutant-r2 D. melanogaster, fly stocks (stock 507) were obtained from the Bloomington Stock Center. Eye–antennal and wing imaginal disks were dissected. RNA was extracted, yielding 3 mg of total RNA per sample. The samples were processed into libraries according to the Illumina TruSeq protocol with appropriate indices, pooled. Sequencing of the transcriptomes was performed on the Illumina HiSeq 2000 platform

RNA-Seq RNA-Seq

After sequencing, raw reads from the samples were aligned to the human genome (GRCh37/hg19) using the RNA-Seq Unified Mapper (version 2.0.4) with default parameters. Only reads that mapped to a single location were used (selected from the RUM_Unique aligned reads file) After sequencing, Fastx-clipper (Gordon and Hannon, 2010) was used to discard reads containing residuals of adapter sequences that were discarded (FastX clipper version 0.0.13 with option-M15). Quality control was applied to the raw sequence reads and performed using the FastQC software (version 0.9; Andrews, 2010), checking for PHRED quality >20 and different primer contaminations. The reads were then aligned using TopHat version 2.0 (Pachter et al, 2009) with default parameters to the Flybase D. melanogaster genome version r5.45 (released March 2012; Antonazzo et al, 2017)

Homo sapiens (IVT-Plasmids, IVT-Only, and IVT-PolyAsel)	Drosophila melanogaster (wild type and mutant)
Replicates: 1×.sam file for each IVT sample (no replicates)	Replicates: 4×.sam files (two replicates for each species)
Sample preparation	Sample preparation
Glycerol stocks containing individual cDNAs (cloned into pCMV-Sport 6 plasmid) from the MGC (Derge et al, 2009) were produced. Plasmid DNA was extracted from these glycerol stocks and plated at 50 ng per well in 384-well plates. The contents of three 384-well plates (total of 1062 human transcripts) were collected. The plasmid library was then amplified by transferring 10 ng into Escherichia coli DH5α cells (cat. no. 18258012; Invitrogen, Life Technologies, Carlsbad, CA). The heat shock method was used to transform E. coli (see Lahens et al, 2014, for more details). The plasmids were then purified using Qiagen (Hilden, Germany) Maxiprep Kit (cat. no. 12163) according to the manufacturer's protocol. Samples were sequenced on the Illumina HiSeq 2000 platform	For the wild-type D. melanogaster, fly stocks (Canton-S and strain RAL-208) were obtained from the inbred collection of T. Mackay. For the mutant-r2 D. melanogaster, fly stocks (stock 507) were obtained from the Bloomington Stock Center. Eye–antennal and wing imaginal disks were dissected. RNA was extracted, yielding 3 mg of total RNA per sample. The samples were processed into libraries according to the Illumina TruSeq protocol with appropriate indices, pooled. Sequencing of the transcriptomes was performed on the Illumina HiSeq 2000 platform
RNA-Seq	RNA-Seq
After sequencing, raw reads from the samples were aligned to the human genome (GRCh37/hg19) using the RNA-Seq Unified Mapper (version 2.0.4) with default parameters. Only reads that mapped to a single location were used (selected from the RUM_Unique aligned reads file)	After sequencing, Fastx-clipper (Gordon and Hannon, 2010) was used to discard reads containing residuals of adapter sequences that were discarded (FastX clipper version 0.0.13 with option-M15). Quality control was applied to the raw sequence reads and performed using the FastQC software (version 0.9; Andrews, 2010), checking for PHRED quality >20 and different primer contaminations. The reads were then aligned using TopHat version 2.0 (Pachter et al, 2009) with default parameters to the Flybase D. melanogaster genome version r5.45 (released March 2012; Antonazzo et al, 2017)

Left: H. sapiens IVT RNA-Seq, Right: D. melanogaster.

cDNAs, complementary DNAs; IVT, in vitro transcribed; MGC, Mammalian Gene Collection.

3. RESULTS

3.1. Analysis of IVT RNA in H. sapiens

The first IVT sample we analyzed was the IVT-Plasmids sample, as this was produced from sequencing the Human IVT plasmids directly without applying ribosomal depletion or PolyA selection methods and therefore represents our control. The IVT-Plasmids sample, by virtue of not having RNA selection protocol steps applied, also reduces the technical sources of variation in read distribution. Table 2 shows a number of 4-mer motif-pairs with the highest or lowest correlations. We construct box plots of the correlations across the IVT-Plasmids sample, and we partitioned the results as a function of GC content of the motif and GC content of the exon and produced box and whisker plots (Fig. 3). We observe that the median correlation has a weak dependence on the motif GC content and is in the range of 0.6–0.8, with the exception of 100% motif GC content. The median correlation as a function of the GC content of the exons likewise lies in the range of 0.6–0.8.

FIG. 3.

Homo sapiens IVT-Plasmid replicates Pearson's outliers box and whisker plot of the IVT-Plasmid sample (control) correlations as a function of motif GC and mean exon GC content, for all spacings. GC, guanosine-cytosine; IVT, in vitro transcribed.

Table 2.

Homo sapiens Pearson's Correlation Coefficient Outliers (Top 10 and Lowest 10) for Different Intra-Exon 4-Mer Motif Sequence Pairs at 10, 50, 100, and 200 bp Spacings

IVT-Plasmids H. sapiens
R (10 bp)	R (50 bp)	R (100 bp)	R (200 bp)
Lowest 10 Pearson's correlation outliers and their motifs
ATCG = −0.0543 (23)	ACGA = −0.2715 (15)	TCGA = −0.2600 (16)	ACGA = −0.3050 (10)
CGTC = −0.0209 (44)	GACG = −0.0101 (39)	CGAA = −0.1953 (17)	TCGC = −0.1834 (13)
ACCG = −0.0129 (49)	AACG = 0.0000 (9)	CGTT = −0.1510 (10)	ACCG = −0.1608 (18)
CGTT = −0.0052 (21)	TTCG = 0.0000 (9)	GTCG = −0.0786 (19)	GTAC = −0.1465 (21)
CGTA = 0.0000 (7)	TCGA = 0.0000 (8)	CGGT = −0.0304 (45)	TCGA = −0.1274 (11)
CGAC = 0.0479 (43)	CGTT = 0.0000 (8)	CGTA = 0.0000 (5)	CGAG = −0.1055 (41)
ACGC = 0.1639 (40)	TCGT = 0.0000 (8)	ACGT = 0.0000 (9)	TGCG = −0.0935 (40)
GCGA = 0.1665 (62)	TACG = 0.0000 (7)	TACG = 0.0000 (4)	GTCT = −0.0547 (54)
ACCC = 0.1899 (260)	CGAT = 0.0000 (9)	CCGT = 0.0312 (37)	TAGG = −0.0541 (24)
GGGG = 0.1991 (649)	TCGC = 0.0884 (28)	CGAT = 0.0407 (18)	GCGA = −0.0482 (18)
Highest 10 Pearson's correlation outliers and their motifs
TAGA = 0.9985 (48)	TAGG = 0.9943 (19)	ATCG = 0.9822 (16)	TCAA = 0.8609 (57)
GTCG = 0.9984 (19)	ACTA = 0.9703 (21)	GTAT = 0.9581 (33)	ATCG = 0.8324 (10)
TTAG = 0.9983 (39)	CTAG = 0.9692 (19)	CGTC = 0.9429 (31)	AGGC = 0.7683 (137)
ACGT = 0.9978 (19)	CCTA = 0.9672 (26)	GTTA = 0.9233 (23)	ATCC = 0.7553 (49)
ATAG = 0.9976 (37)	CATT = 0.9644 (103)	AGTC = 0.9189 (59)	GGAT = 0.7430 (40)
ATTC = 0.9974 (85)	ACTT = 0.9594 (84)	GTCA = 0.8984 (67)	TGAC = 0.6765 (52)
GCAC = 0.9971 (110)	GATA = 0.9590 (27)	TCGT = 0.8971 (12)	TACC = 0.6553 (22)
TAGG = 0.9970 (31)	CTAA = 0.9572 (31)	GCAA = 0.8938 (72)	TGGT = 0.6435 (84)
TCGT = 0.9968 (18)	TATC = 0.9495 (27)	CATC = 0.8817 (117)	GGTG = 0.6410 (103)
ACAT = 0.9968 (114)	TTAC = 0.9389 (33)	ATTC = 0.8620 (78)	ACGG = 0.6184 (21)

Actual motif-pair counts are in parentheses.

To compare the effect of applying different RNA selection protocol methods to the IVT samples, specifically ribosomal depletion versus PolyA selection, we compared the IVT-Only and IVT-PolyA samples, respectively. Tables 3 and 4 present that the highest outliers for both these IVT-Seq samples show a number of 4-mer motif pairs that have very high correlations. High correlation outliers are observed across all spacings. We produced box and whisker plots of correlation as a function of GC content of the motif and GC content of the exon (Fig. 4). The trend of the data for IVT-Only and IVT-PolyA, as a function of motif GC content and mean GC content, is similar to that of the IVT-Plasmids data, although the median correlation is somewhat less (0.4–0.6).

FIG. 4.

Correlation (Pearson's) as a function of 4-mer motif and exon GC content, across all spacings, for Homo sapiens for the two IVT samples that had different RNA selection protocols applied: (Top) IVT-Only sample, which underwent ribosomal depletion, and (Bottom) IVT-PolyA, which underwent PolyA selection.

Table 3.

Homo sapiens In Vitro Transcribed-Seq Pearson's Correlation Coefficient Outliers (Top 10 and Lowest 10) for Different Intra-Exon 4-Mer Motif Sequence Pairs at 10, 50, 100, and 200 bp Spacings: In Vitro Transcribed-Only Library Preparation

IVT-Only H. sapiens
R (10 bp)	R (50 bp)	R (100 bp)	R (200 bp)
Lowest 10 Pearson's correlation outliers and their motifs
CGTT = −0.1726 (18)	CGAC = −0.2017 (22)	CGAC = −0.1421 (28)	CTAC = −0.2360 (30)
CGTA = −0.1470 (13)	GCGT = −0.1374 (25)	GCGA = −0.1255 (33)	ACGC = −0.2332 (14)
ACGC = −0.1337 (35)	CGAT = −0.1171 (11)	ACCG = −0.1224 (25)	TCGG = −0.2111 (21)
CCGT = −0.0979 (37)	CCTA = −0.1010 (34)	GACG = −0.0923 (29)	GTAT = −0.1806 (13)
GCGT = −0.0746 (44)	TAAC = −0.0980 (42)	CGAA = −0.0903 (18)	CCGA = −0.1790 (17)
GAGT = −0.0628 (101)	CGAA = −0.0884 (18)	GTTA = −0.0891 (32)	CGGT = −0.1735 (16)
CGTC = −0.0298 (44)	TAGG = −0.0698 (32)	CGTT = −0.0832 (11)	ACGG = −0.1720 (16)
CGGA = −0.0217 (75)	AACC = −0.0657 (81)	CGGT = −0.0666 (34)	ACTA = −0.1571 (23)
GCGA = −0.0156 (52)	CGCA = −0.0563 (40)	GTTC = −0.0627 (48)	GACG = −0.1540 (26)
TACG = 0.0000 (6)	GCTA = −0.0562 (26)	GTAG = −0.0583 (40)	TGCG = −0.1532 (34)
Highest 10 Pearson's correlation outliers and their motifs
TACA = 0.9989 (98)	TCGA = 0.9805 (15)	ATCC = 0.9609 (41)	TAGT = 0.9913 (20)
CTCA = 0.9983 (168)	GACT = 0.9713 (56)	AGTA = 0.9403 (31)	ATCC = 0.9860 (35)
AGTA = 0.9978 (50)	TGCA = 0.9674 (114)	TATA = 0.9261 (50)	GCAT = 0.9734 (40)
TACT = 0.9977 (47)	AGGT = 0.9667 (70)	CGTG = 0.9146 (39)	CTCG = 0.9417 (19)
TGAA = 0.9975 (229)	ATGC = 0.9626 (44)	TCAC = 0.8733 (77)	AGTC = 0.8869 (34)
GTCA = 0.9970 (81)	CCTT = 0.9551 (132)	CACC = 0.8694 (164)	TCAC = 0.8492 (47)
GCAT = 0.9969 (92)	TTGC = 0.9544 (48)	CACT = 0.8642 (109)	TGTG = 0.8017 (111)
AAGA = 0.9969 (257)	ATGG = 0.9543 (93)	CCCT = 0.8530 (275)	ATTC = 0.7894 (52)
AACT = 0.9965 (104)	CTCG = 0.9507 (30)	TGAT = 0.8525 (64)	CTGT = 0.7682 (114)
GTAT = 0.9963 (24)	CCGT = 0.9480 (17)	GTAC = 0.8448 (18)	TACC = 0.7581 (22)

Sample sizes are given in parentheses.

Table 4.

Homo sapiens In Vitro Transcribed-Seq Pearson Correlation Coefficient Outliers (Top 10 and Lowest 10) for Different Intra-Exon 4-Mer Motif Sequence Pairs at 10, 50, 100, and 200 bp Spacings: In Vitro Transcribed with PolyA Library Selection

IVT-PolyAsel H. sapiens
R (10 bp)	R (50 bp)	R (100 bp)	R (200 bp)
Lowest 10 Pearson's correlation outliers and their motifs
CGTA = −0.1953 (11)	ATCG = −0.2190 (16)	CGAA = −0.2046 (13)	TGCG = −0.2025 (42)
CGTT = −0.1012 (20)	CGTT = −0.2000 (10)	TACC = −0.1573 (28)	ATCG = −0.1998 (12)
GCGA = −0.0976 (67)	AACG = −0.1797 (14)	TCGT = −0.1174 (25)	CGTT = −0.1923 (13)
ATCG = −0.0694 (31)	TCTA = −0.1429 (27)	GCGA = −0.1062 (34)	CGAC = −0.1810 (17)
AACG = −0.0510 (20)	ACTA = −0.1184 (23)	CGCA = −0.0956 (84)	ATAG = −0.1664 (21)
GCGT = −0.0388 (46)	CGCA = −0.1003 (59)	ATAG = −0.0786 (49)	CGTG = −0.1445 (49)
CTTA = −0.0373 (44)	GTTA = −0.0876 (31)	GTTG = −0.0708 (115)	GTTA = −0.1404 (40)
TAGC = −0.0171 (44)	CGAC = −0.0747 (40)	TCGC = −0.0676 (46)	AACG = −0.1316 (13)
GGTA = −0.0107 (66)	TACC = −0.0682 (13)	GACG = −0.0637 (41)	CGGT = −0.1314 (32)
GAGT = −0.0079 (110)	GCTA = −0.0583 (35)	TCGG = −0.0576 (63)	GTAG = −0.1307 (41)
Highest 10 Pearson's correlation outliers and their motifs
TACT = 0.9998 (45)	ACTT = 0.9887 (68)	AGTA = 0.9703 (36)	ATCC = 0.9890 (41)
TACC = 0.9997 (39)	AGAG = 0.9886 (129)	CGAT = 0.9035 (18)	ACGA = 0.9751 (11)
CTAC = 0.9992 (60)	CTAA = 0.9726 (25)	TATA = 0.8982 (88)	TAGT = 0.9358 (33)
TAAC = 0.9991 (27)	CACG = 0.9685 (37)	CGTG = 0.8641 (65)	GCAT = 0.8991 (45)
AGTA = 0.9986 (57)	TGGA = 0.9652 (196)	ACAC = 0.8495 (75)	CGAG = 0.7917 (30)
AGAC = 0.9979 (81)	AAGC = 0.9606 (85)	TCAA = 0.8321 (93)	TGTG = 0.7810 (119)
TGAC = 0.9976 (84)	TTGA = 0.9571 (82)	CATA = 0.8264 (39)	GACG = 0.7799 (29)
CCTA = 0.9976 (47)	AACT = 0.9559 (87)	ACCT = 0.8171 (93)	TGAA = 0.7644 (123)
TGAA = 0.9958 (234)	ATAA = 0.9541 (77)	CACT = 0.8108 (114)	AATC = 0.7620 (35)
ATCC = 0.9956 (71)	ATAG = 0.9531 (42)	CCAA = 0.7891 (124)	CATC = 0.7204 (72)

Sample sizes are given in parentheses.

3.2. Analysis of wild-type and mutant-r2 D. melanogaster

We analyzed the whole transcriptomes of two Drosophila fruit fly (species D. melanogaster), which only differ by mutation gl[60j] in the eye–antennal disk (Naval-Sánchez et al, 2013). Overall, outliers for correlations for both wild type and mutant-r2 are given in Tables 5 –8. We observe much lower correlations for the mutant-r2 than in wild type—the difference is especially marked when comparing the motifs with the highest 10 correlations in both samples for all separations.

Table 5.
Drosophila melanogaster Pearson's Correlation Coefficient Outliers (Top 10 and Lowest 10) for Different Intra-Exon 4-Mer Motif Sequence Pairs at 10, 50, 100, and 200 bp Spacings: Wild Type

Wild-type D. melanogaster

R (10 bp) R (50 bp) R (100 bp) R (200 bp)

Lowest 10 Pearson's correlation outliers and their motifs

TAGG = −0.0151 (539) CCTA = −0.0074 (403) GGGG = 0.0024 (1991) GGGG = −0.0187 (1525)

CATA = 0.0108 (1150) GGTT = −0.0051 (1243) CCTA = 0.0066 (627) CCCC = −0.0043 (1721)

CCGA = 0.0250 (1343) CTAG = −0.0001 (331) TACC = 0.0203 (951) ACCC = 0.0030 (1346)

TAAC = 0.0268 (790) CCGG = 0.0225 (2795) GGTT = 0.0259 (2008) AACC = 0.0063 (1559)

CCCC = 0.0461 (2356) CATA = 0.0248 (636) AACC = 0.0460 (2028) CCGG = 0.0182 (3351)

CGTA = 0.0624 (620) GATA = 0.0479 (664) CGCG = 0.0535 (2417) GGGT = 0.0292 (1310)

CGGG = 0.0627 (2534) TAAC = 0.0485 (464) GGTA = 0.0577 (868) GGTA = 0.0350 (680)

GGTA = 0.0639 (779) AGGG = 0.0522 (1313) GGGA = 0.0651 (2585) GTCT = 0.0367 (691)

GGGA = 0.0710 (2197) CGCG = 0.0635 (1551) TCCC = 0.0728 (2637) GGTT = 0.0384 (1674)

AGGG = 0.0832 (1702) AGCC = 0.0744 (1243) TGGG = 0.0800 (3963) CCCA = 0.0392 (3348)

Highest 10 Pearson's correlation outliers and their motifs

TAGA = 0.9860 (602) TAGA = 0.8947 (445) AGTG = 0.9289 (1402) GAGT = 0.8516 (975)

ATAG = 0.9672 (894) GTAC = 0.8623 (529) TATG = 0.8375 (978) ACGG = 0.7927 (837)

TAGC = 0.9325 (729) TTAG = 0.8317 (551) CTAA = 0.8110 (753) ACTA = 0.7202 (596)

AGTG = 0.9002 (1368) AGTG = 0.8158 (920) GCTA = 0.8054 (750) TACG = 0.7173 (517)

AACG = 0.8898 (1330) CTAT = 0.8059 (533) CCAT = 0.7493 (2080) TCGT = 0.7171 (1541)

ATCC = 0.8336 (2306) ACTA = 0.8042 (515) GCTG = 0.7303 (9467) CTAC = 0.7036 (749)

CTGA = 0.8192 (1577) AGCA = 0.7986 (4465) GACT = 0.7296 (1149) ACAG = 0.6966 (1360)

CACT = 0.8175 (1706) GTGT = 0.7916 (1063) TAGT = 0.7249 (802) ACCA = 0.6960 (2168)

CGAG = 0.8028 (2433) CAGT = 0.7776 (1182) GAGT = 0.7149 (1153) CGTA = 0.6860 (560)

AATT = 0.8018 (3447) AGTT = 0.7770 (1506) CGAA = 0.7145 (2171) TGTG = 0.6842 (1867)

Wild-type D. melanogaster
Lowest 10 Pearson's correlation outliers and their motifs
TAGG = −0.0151 (539)	CCTA = −0.0074 (403)	GGGG = 0.0024 (1991)	GGGG = −0.0187 (1525)
CATA = 0.0108 (1150)	GGTT = −0.0051 (1243)	CCTA = 0.0066 (627)	CCCC = −0.0043 (1721)
CCGA = 0.0250 (1343)	CTAG = −0.0001 (331)	TACC = 0.0203 (951)	ACCC = 0.0030 (1346)
TAAC = 0.0268 (790)	CCGG = 0.0225 (2795)	GGTT = 0.0259 (2008)	AACC = 0.0063 (1559)
CCCC = 0.0461 (2356)	CATA = 0.0248 (636)	AACC = 0.0460 (2028)	CCGG = 0.0182 (3351)
CGTA = 0.0624 (620)	GATA = 0.0479 (664)	CGCG = 0.0535 (2417)	GGGT = 0.0292 (1310)
CGGG = 0.0627 (2534)	TAAC = 0.0485 (464)	GGTA = 0.0577 (868)	GGTA = 0.0350 (680)
GGTA = 0.0639 (779)	AGGG = 0.0522 (1313)	GGGA = 0.0651 (2585)	GTCT = 0.0367 (691)
GGGA = 0.0710 (2197)	CGCG = 0.0635 (1551)	TCCC = 0.0728 (2637)	GGTT = 0.0384 (1674)
AGGG = 0.0832 (1702)	AGCC = 0.0744 (1243)	TGGG = 0.0800 (3963)	CCCA = 0.0392 (3348)
Highest 10 Pearson's correlation outliers and their motifs
TAGA = 0.9860 (602)	TAGA = 0.8947 (445)	AGTG = 0.9289 (1402)	GAGT = 0.8516 (975)
ATAG = 0.9672 (894)	GTAC = 0.8623 (529)	TATG = 0.8375 (978)	ACGG = 0.7927 (837)
TAGC = 0.9325 (729)	TTAG = 0.8317 (551)	CTAA = 0.8110 (753)	ACTA = 0.7202 (596)
AGTG = 0.9002 (1368)	AGTG = 0.8158 (920)	GCTA = 0.8054 (750)	TACG = 0.7173 (517)
AACG = 0.8898 (1330)	CTAT = 0.8059 (533)	CCAT = 0.7493 (2080)	TCGT = 0.7171 (1541)
ATCC = 0.8336 (2306)	ACTA = 0.8042 (515)	GCTG = 0.7303 (9467)	CTAC = 0.7036 (749)
CTGA = 0.8192 (1577)	AGCA = 0.7986 (4465)	GACT = 0.7296 (1149)	ACAG = 0.6966 (1360)
CACT = 0.8175 (1706)	GTGT = 0.7916 (1063)	TAGT = 0.7249 (802)	ACCA = 0.6960 (2168)
CGAG = 0.8028 (2433)	CAGT = 0.7776 (1182)	GAGT = 0.7149 (1153)	CGTA = 0.6860 (560)
AATT = 0.8018 (3447)	AGTT = 0.7770 (1506)	CGAA = 0.7145 (2171)	TGTG = 0.6842 (1867)

Sample sizes are given in parentheses.

Table 6.

Drosophila melanogaster Pearson's Correlation Coefficient Outliers (Top 10 and Lowest 10) for Different Intra-Exon 4-Mer Motif Sequence Pairs at 10, 50, 100, and 200 bp Spacings: Wild-Type-r2

Wild-type-r2 D. melanogaster
R (10 bp)	R (50 bp)	R (100 bp)	R (200 bp)
Lowest 10 Pearson's correlation outliers and their motifs
TAGG = −0.0154 (534)	CCTA = −0.0078 (392)	GGGG = 0.0023 (1977)	GGGG = −0.0188 (1507)
CATA = 0.0100 (1141)	CTAG = 0.0005 (320)	CCTA = 0.0065 (611)	CCCC = −0.0045 (1699)
CCGA = 0.0248 (1337)	GGTT = 0.0173 (757)	CGCG = 0.0533 (2398)	ACCC = 0.0030 (1332)
TAAC = 0.0267 (781)	CATA = 0.0204 (630)	GGGA = 0.0648 (2538)	CCGG = 0.0185 (3330)
CCCC = 0.0460 (2346)	CCGG = 0.0224 (2776)	AACC = 0.0713 (1206)	GGGT = 0.0288 (1286)
CGTA = 0.0622 (613)	GATA = 0.0479 (662)	TCCC = 0.0736 (2599)	GTCT = 0.0370 (684)
CGGG = 0.0625 (2519)	TAAC = 0.0484 (456)	TGGG = 0.0797 (3909)	CCCA = 0.0391 (3305)
GGGA = 0.0713 (2185)	AGGG = 0.0520 (1302)	CGGG = 0.0904 (3318)	AACC = 0.0425 (957)
AGGG = 0.0831 (1690)	CGCG = 0.0629 (1524)	TCTC = 0.0930 (1669)	AGGG = 0.0479 (1509)
CGTT = 0.0844 (1242)	AGCC = 0.0747 (1237)	CATA = 0.1012 (1018)	ATAG = 0.0550 (532)
Highest 10 Pearson's correlation outliers and their motifs
TAGA = 0.9860 (600)	TAGA = 0.9033 (440)	TATG = 0.8378 (964)	GAGT = 0.8517 (961)
ATAG = 0.9672 (892)	TACC = 0.8850 (336)	GCTA = 0.8293 (743)	ACGG = 0.7927 (836)
TAGC = 0.9325 (727)	GTAC = 0.8623 (523)	CTAA = 0.8113 (738)	ACTA = 0.7205 (588)
AACG = 0.8898 (1322)	GGTA = 0.8370 (357)	TAGT = 0.7988 (787)	TCGT = 0.7189 (1534)
AGTG = 0.8618 (1575)	TTAG = 0.8317 (544)	CCAT = 0.7495 (2063)	TACG = 0.7173 (514)
ATCC = 0.8336 (2296)	AGTG = 0.8219 (1072)	GCTG = 0.7304 (9423)	CTAC = 0.7047 (740)
CTGA = 0.8193 (1555)	CTAT = 0.8069 (518)	AGTG = 0.7304 (1610)	ACCA = 0.6961 (2155)
CACT = 0.8175 (1694)	ACTA = 0.8049 (509)	GACT = 0.7297 (1132)	CGTA = 0.6874 (552)
CGAG = 0.8028 (2426)	AGCA = 0.7988 (4445)	CGAA = 0.7159 (2143)	TGTG = 0.6844 (1849)
AATT = 0.8019 (3416)	TCTA = 0.7966 (342)	GAGT = 0.7149 (1141)	GAAA = 0.6757 (2726)

Sample sizes are given in parentheses.

Table 7.

Drosophila melanogaster Pearson's Correlation Coefficient Outliers (Top 10 and Lowest 10) for Different Intra-Exon 4-Mer Motif Sequence Pairs at 10, 50, 100, and 200 bp Spacings: Mutant-r2

Mutant-r2 D. melanogaster
R (10 bp)	R (50 bp)	R (100 bp)	R (200 bp)
Lowest 10 Pearson's correlation outliers and their motifs
CCTA = −0.0163 (932)	GACT = −0.0208 (1403)	CCTA = −0.0116 (1083)	TACT = −0.0230 (1301)
TCTA = −0.0137 (1346)	GCTA = −0.0146 (1090)	TCCC = −0.0072 (4951)	GTAA = −0.0215 (1321)
TAGG = −0.0112 (944)	TAAC = −0.0124 (1099)	TCTA = −0.0071 (1505)	GGTA = −0.0203 (1135)
CCGA = −0.0083 (2821)	CATA = −0.0120 (1278)	ACTA = −0.0039 (1535)	AACC = −0.0166 (2747)
TAAC = −0.0040 (1623)	CTAC = −0.0114 (1345)	GGGG = 0.0025 (4280)	TTAC = −0.0150 (1397)
GTAG = −0.0039 (1736)	GTCT = −0.0108 (1466)	GGGA = 0.0091 (4967)	GTGT = −0.0113 (2178)
CATA = 0.0017 (2064)	ATAC = −0.0108 (1180)	GATT = 0.0103 (4102)	GTTA = −0.0107 (1368)
CTAG = 0.0018 (805)	GAGT = −0.0081 (1634)	GTTT = 0.0108 (6502)	TCAC = −0.0091 (2027)
GACC = 0.0025 (2731)	TAAG = −0.0080 (1201)	ACGT = 0.0115 (1450)	TAAC = −0.0081 (1267)
GATA = 0.0032 (1784)	CCTA = −0.0073 (701)	CATA = 0.0115 (2078)	GTGA = −0.0073 (1933)
Highest 10 Pearson's correlation outliers and their motifs
TAGT = 0.7594 (1485)	TACC = 0.5456 (1016)	CAAG = 0.3938 (3971)	CTTG = 0.5247 (3120)
CACT = 0.6408 (1754)	ATGT = 0.5207 (1390)	CTTG = 0.3705 (4073)	GCAG = 0.4106 (6614)
GATG = 0.5231 (3671)	CTTG = 0.4582 (2911)	CTGC = 0.3625 (8133)	TGAT = 0.3769 (1844)
TGCA = 0.4906 (3256)	TTAG = 0.4365 (1190)	TGAT = 0.3560 (2331)	CAAG = 0.3720 (3118)
CAAG = 0.4242 (3575)	CAAG = 0.3459 (2740)	CATC = 0.3467 (3747)	ACAT = 0.3496 (1661)
GCAG = 0.3975 (8213)	AGTT = 0.3399 (1555)	GATG = 0.3396 (3686)	CGAA = 0.3235 (3096)
AAGG = 0.3903 (4231)	TATG = 0.3257 (1270)	TATC = 0.3316 (2020)	GATC = 0.3184 (2394)
TCGA = 0.3841 (2228)	TTGT = 0.3235 (2925)	ATCA = 0.3291 (2324)	CATT = 0.3037 (2055)
ACTC = 0.3804 (2157)	AACT = 0.3221 (1459)	ATCG = 0.3053 (2304)	TGTG = 0.3024 (1977)
AATG = 0.3743 (2211)	GATG = 0.3199 (2601)	CTCG = 0.2964 (2798)	ACGG = 0.3001 (2036)

Sample sizes are given in parentheses.

Table 8.

Drosophila melanogaster Pearson's Correlation Coefficient Outliers (Top 10 and Lowest 10) for Different Intra-Exon 4-Mer Motif Sequence Pairs at 10, 50, 100, and 200 bp Spacings: Mutant-r3

Mutant-r3 type D. melanogaster
R (10 bp)	R (50 bp)	R (100 bp)	R (200 bp)
Lowest 10 Pearson's correlation outliers and their motifs
GTAG = −0.0110 (4451)	GACT = −0.0155 (3709)	CCTA = −0.0079 (3680)	TACT = −0.0167 (3730)
TCTA = −0.0109 (4014)	ATAC = −0.0121 (3061)	ACTA = −0.0067 (3991)	TTAC = −0.0142 (3633)
CCTA = −0.0099 (3058)	CATA = −0.0108 (3062)	ACGT = −0.0028 (4259)	GGTA = −0.0133 (3456)
TAGG = −0.0090 (3009)	TAAC = −0.0100 (2946)	TCCC = −0.0013 (14,111)	AGTG = −0.0114 (5061)
CCGA = −0.0085 (7526)	TAAG = −0.0086 (3286)	TCTA = −0.0010 (5028)	TCAC = −0.0113 (5532)
TAAC = −0.0067 (3901)	TCAT = −0.0073 (5180)	TAGA = −0.0007 (4971)	GTTA = −0.0104 (3718)
CGTA = −0.0057 (2955)	GCTA = −0.0069 (3236)	CGTA = 0.0014 (3747)	CTCA = −0.0095 (5707)
GATA = −0.0045 (4236)	GATA = −0.0059 (3368)	ACGC = 0.0019 (7405)	GTAA = −0.0093 (3437)
CTAC = −0.0034 (4526)	ATTG = −0.0058 (5479)	GATT = 0.0027 (8208)	GTGT = −0.0092 (5366)
CTAG = −0.0025 (3236)	CCTA = −0.0058 (2559)	TAGT = 0.0027 (3768)	AACC = −0.0084 (7436)
Highest 10 Pearson's correlation outliers and their motifs
TAGT = 0.5918 (3456)	TACC = 0.2893 (2881)	TAAC = 0.1998 (4402)	CGAA = 0.1921 (7254)
AAGG = 0.2979 (9354)	TAGC = 0.2644 (3270)	CAAG = 0.1738 (12,860)	ACGG = 0.1903 (5536)
CACT = 0.2686 (5633)	CTTG = 0.1926 (8830)	GCGT = 0.1731 (7475)	CAAG = 0.1680 (10,024)
TATC = 0.2670 (4379)	GTGT = 0.1847 (4722)	TATC = 0.1621 (5203)	TCAA = 0.1667 (8631)
TGGA = 0.2447 (13,432)	TATG = 0.1786 (3069)	AACA = 0.1608 (14,808)	CCCT = 0.1633 (9198)
CAAG = 0.2332 (10,781)	CGGA = 0.1735 (7791)	CCCG = 0.1565 (14,386)	AGGA = 0.1460 (11,653)
ATCC = 0.2281 (8329)	CAAG = 0.1729 (8732)	ATCG = 0.1526 (7561)	CTTA = 0.1453 (3943)
CCCT = 0.2256 (9530)	GTTT = 0.1697 (8514)	ACCA = 0.1421 (13,261)	CTTG = 0.1422 (10,212)
AATT = 0.2210 (10,274)	AATT = 0.1612 (7287)	CTAG = 0.1413 (4125)	GTAC = 0.1415 (3391)
TACG = 0.2165 (2974)	GGTA = 0.1547 (2923)	CCAA = 0.1409 (14,643)	TACC = 0.1394 (3586)

Sample sizes are given in parentheses.

3.2.1. GC content of the replicates

We show box plots of the intra-exon 4-mer motif-pair correlations as a function of both motif and mean exon GC content, for both D. melanogaster data sets (Fig. 5). We note in the wild-type data set that the median correlation shows a dependence on the GC content of the motif, with the median correlations having a range of 0.1–0.5.

FIG. 5.

Correlation (Pearson's) as a function of 4-mer motif and exon GC content in (A) wild-type, (B) wild-type-r2, (C) mutant-r2, and (D) mutant-r3 Drosophila melanogaster transcriptomes.

However, in the mutant-r2 data set, the median correlations in all cases are <0.1 and independent of the GC content of the motif or the exon. In the mutant-r3 data set, the range of correlations is consistent with zero (within the first and third quartiles).

3.3. FastQC and BamQC—Quality control analysis of D. melanogaster data set

The presence of raw unaligned read data for the above data sets allows us to compare these results with more traditional measures of quality control.

With this in mind, we ran FastQC (Andrews, 2010) on all the data sets we analyzed in this research. We note that relevant sets of reads for the IVT-Plasmids are not available.

We also ran Qualimap to analyze the source BAM files that the SAM files were directly derived from the D. melanogaster data set to see if there were any alignment or coverage issues in the mutant data (see Section C, Supplementary Table S2, and Supplementary Fig. S4 in the Supplementary Data).

In Figure 6, we plot the heat maps of the Sequencing per-tile Means of the IVT-PolyA and IVT-NoSel reads. These depict deviations from the average tile quality within areas of flowcells on the sequencing apparatus (the explicit values underlying the heat maps are referred to as the per-tile Means; Babraham Institute, 2015). Both data sets demonstrate a low deviation from the average time quality of for each tile.

FIG. 6.

Heat maps of the Sequencing per-tile Means, obtained from FastQC analysis for the two RNA selection Homo sapiens IVT RNA-Seq samples analyzed—sequencing tile information was not present in the IVT-Plasmids sample. (Top) IVT-PolyA, (bottom) IVT-NoSel.

However, in comparing the FastQC analysis results for the data sets, we found that the per-tile sequencing quality heat maps were radically different for the mutant and wild-type data sets (Fig. 7). In the wild-type replicates, these heat maps were plain and uniform, whereas both of the mutant replicates showed considerable variation from the Means. There were no differences in any of the other FastQC analysis parameters measured between the wild-type and mutant data sets.

FIG. 7.

Heat maps of the Sequencing per-tile Means, obtained from FastQC analysis for all Drosophila melanogaster RNA-Seq samples and their replicates analyzed. (Top) wild-type replicates, (Bottom) mutant replicates.

To investigate this effect further, we plotted the distribution of per-tile Means as both a histogram and a density plot (Fig. 8) for all the samples in all the data sets we have analyzed. In Figure 8, we observe that the distributions in the IVT and wild-type data sets show the least deviation from the Mean tile qualities. The mutant D. melanogaster data sets show considerable deviation from the Mean tile qualities—which match the low correlations seen across the mutant replicates. By computing intra-exon motif-pair correlations, we have observed that sequencing errors occurring in flowcell tiles result in widespread deviations from the uniformity of mapped reads across exons and are independent of sequence or GC content.

FIG. 8.

Visualization of the distributions of the Sequencing per-tile Means, obtained from FastQC analysis for all RNA-Seq samples analyzed (species Drosophila melanogaster and Homo sapiens). (A) Histogram. (B) Density plot. The distributions in the IVT data sets, and both wild-type replicates show that less of the tiles in the flowcells of the sequencing apparatus deviate from the per-tile Means. This is indicated by the tallest narrowest peaks, centered around 0, for the IVT-PolyA data set followed by IVT-NoSel, and then the first wild-type replicate. The second wild-type replicate was the third tallest peak and showed a reasonably good narrow distribution of the data about the per-tile Means, although lower than the IVT data sets and first wild-type replicate. The two mutant replicates, however, show very wide distributions about the per-tile Means in the flowcells, as depicted by the two shorter wider curves respectively. NB: The IVT-Plasmids sample did not carry sequencing tile information.

In the BamQC data, there is not much difference between the mutant and wild-type data in terms of the mapping quality, nucleotide frequencies, and GC content, although we note an increase in variance in the coverage of the mutant D. melanogaster data set as indicated by the slightly higher standard deviation.

4. DISCUSSION

We have analyzed RNA-Seq samples from two species: H. sapiens (MGC plasmids) and D. melanogaster (transcriptomes). The IVT data set was produced in a highly controlled manner—we have analyzed three H. sapiens samples from this, two were subjected to mRNA selection before RNA-Seq, and the other was not. The IVT samples contain fewer reads than the D. melanogaster samples in which the wild and mutant-r2 types comprise ∼12.9 and 15.0 M reads respectively, whereas the IVT-Only sample has 0.406 M reads, the IVT-PolyA sample has 0.397 M reads, and the IVT-Plasmids has 0.181 M reads. Although these are comparatively small numbers of reads, the IVT samples have the intra-exon bias within them, as indicated by coverage of reads mapping to the source MGC plasmids used, well quantified.

From the box and whisker plots of the two IVT samples that underwent mRNA selection RNA-Seq protocols (ribosomal depletion and Poly-A selection), we observe a dependence of intra-exon motif pair correlation on motif GC content. This effect is also seen in the wild-type replicates of the D. melanogaster samples, but not in the IVT-Plasmids sample which is not subjected to mRNA selection. The mutant D. melanogaster data have very low overall correlations, and hence, this pattern is not seen here. The median of the IVT-Plasmids correlations is higher than the IVT-PolyA and IVT-Only correlations and that the dependence on the motif GC content is stronger in the latter cases than the former. There is a much more noticeable difference in the median of the correlations between the wild-type and mutant Drosophila data.

4.1. Dependence of intra-exon correlation on GC content appears to be due to mRNA selection

We note that the intra-exon motif correlations as a function of both GC content parameters are much higher in the IVT-Plasmids mRNA selection free RNA-Seq sample than in the other RNA-Seq samples that we analyzed that did undergo mRNA selection: both ribosomal depletion (IVT-Only) and PolyA selection (IVT-PolyA and wild-type). Additionally, both the H. sapiens and wild-type D. melanogaster samples that underwent mRNA selection in the RNA-Seq process had slightly lower correlations than the H. sapiens IVT-Plasmids sample, which did not, suggesting that this is likely of technical origin. Importantly, all the samples from all of the data sets we analyzed were sequenced on the same platform—Illumina HiSeq 2000 (as detailed in Table 1).

As the dependence on overall GC concentration is not observed in the IVT-Plasmids control sample, but is observed in RNA-Seq samples that underwent mRNA selection, we can exclude platform-specific sequencing bias as a source of this effect in the IVT samples. This suggests that not only do mRNA selection protocols result in bias in the distribution of mapped RNA-Seq reads as Lahens et al (2014) demonstrated, but that mRNA selection is also responsible for the dependence of correlation on GC content—we have observed this in all the RNA-Seq samples that underwent mRNA selection across both H. sapiens and D. melanogaster species. Risso et al (2011) also noted motif-specific GC effects, manifesting as deviations from uniform read distribution, which they attribute these problematic motifs being underrepresented.

However, the GC effect we observe occurs in both the wild-type replicates, which have large sample sizes (Table 5), as well as in the IVT-PolyA and IVT-NoSel mRNA selection data sets, which have low sample sizes (due to the controlled way in which the samples were prepared; Tables 3 and 4). Furthermore, the numbers of counts for high GC content motifs, as indicated by the highest outliers (Tables 5–8), are of the same order as other motifs. The dependence of intra-exon correlations on the GC content of the motifs could be explained by mRNA selection methods, which are routinely employed in RNA-Seq experiments and are known to introduce bias (Cao et al, 2002; Lahens et al, 2014; Theurkauf et al, 2012).

5. CONCLUSION

In this work, we have applied our novel k-mer-based analysis method that allows an investigation into RNA-Seq read data at the deep exon level and quantifies sequence-specific deviations in the uniformity of the distribution of mapped reads to a reference genome. The method we have applied uses distributed computing to count reads overlapping 4-mer motifs exhaustively (AAAA through to GGGG) at different regions of the same exon. The assumption is that, when there is no deviation in the distribution of mapped reads, these counts should be correlated. We have demonstrated that the correlations we have computed correspond to mean tile data in the samples from the two data sets we analyzed.

We have shown that sequencing errors occurring in flowcell tiles result in widespread deviations from the uniformity of mapped reads across exons and are independent of sequence or GC content; hence, we propose that this be used when mean tile data are not available and therefore requires only the raw short-read data alone. We have also observed that extremely poor correlations are an indication of technical sources of bias, such as sequencing flowcell tile errors or batch effects. We note also the reduction in the above correlations when not using mRNA selection techniques in the IVT data sets.

This is potentially useful in assessing the quality of RNA-Seq data in public repositories where raw unaligned read data have not been deposited as the above approach can be carried out on a post hoc basis. In practice, if poor or unusual motif correlations are observed, this suggests bias resulting from flowcell errors or batch effects. We propose that this will help with data sets where FastQC cannot be done due to the lack of sequencing metadata such as read-id's and flowcell metadata. For this purposes, readers can use the code developed for this study [Hercules is available on GitHub (Alnasir, 2018) and at Zenodo (Alnasir and Shanahan, 2022)].

When the analysis is performed, plots such as Figure 6 (IVT-PolyA and IVT-NoSel plots) and Figure 7 (Wild-type plots, upper panels) indicate good quality sequencing data, whereas Figure 7 (Mutant-r and Mutant plots, lower panels) indicate potential bias and are a caveat with the data set. Furthermore, these flowcell tile plots can be augmented with plots of the distribution of tile-error means (Fig. 8), where tall narrow distributions suggest good quality data and short flatter distributions suggest bias and poor quality data.

The method and software we have developed are scalable and can, therefore, be applied to analyze the large data sets present in these vast data repositories—as of October 2021 (the time of writing), the SRA alone contains more than 57 peta bases (57.920 × 10¹⁵) of sequencing data, that is in excess of 7 petabytes (Sequence Read Archive, 2017). Application of this method would enable a consistent assessment of quality across such data sets.

Footnotes

ACKNOWLEDGMENTS

We thank Prof. Adrian Shepherd and Dr. Matthew Bashton for suggesting the data sets used in this study as a means to calibrate the method. The authors also wish to thank Eszter Ábrahám for proofreading the article.

AUTHORs' CONTRIBUTIONS

J.J.A. and H.P.S. conceived the study. J.J.A. developed the analysis code. H.P.S. supervised the research and contributed to the writing of the article. Both the authors read and approved the final article.

AUTHOR DISCLOSURE STATEMENT

The authors declare they have no conflicting financial interests.

FUNDING INFORMATION

J.A. was supported through direct funding from the Department of Computer Science at Royal Holloway, University of London.

SUPPLEMENTARY MATERIAL

References

Aerts

. GEO-GSE39781: RNA-Seq in wild-type and glass mutant eye-antennal discs in Drosophila melanogaster. 2012. Available from: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39781 [Last accessed: March 27, 2019].

Alnasir

, Shanahan

. Investigation into the annotation of protocol sequencing steps in the sequence read archive. GigaScience, 2015; 4:23.

Alnasir

, Shanahan

. A novel method to detect bias in short read NGS data. J Integr Bioinform, 2017; 14(3):20170025; doi: 10.1515/jib-2017-0025

Alnasir

. Source code and results data for Transcriptomics Analysis System (Hercules). 2018. Available from: https://github.com/jamie-alnasir/hercules [Last accessed: August 20, 2020].

Alnasir

, Shanahan

Source code and results data for Transcriptomics Analysis System (Hercules). Zenodo Archive, 2022; doi: 10.5281/zenodo.6582436

Alnasir

, Shanahan

. Transcriptomics: Quantifying non-uniform read distribution using mapreduce. Int J Found Comp Sci, 2018; 29(08):1355–1372.

Andrews

, Krueger

, Segonds-Pichon

, et al. FastQC: A quality control tool for high throughput sequence data. 2010:370.

Antonazzo

, Crosby

, Emmert

, et al. Flybase at 25: Looking to the future. Nucleic Acids Res, 2017; 45(Database Issue):D663–D671.

Babraham Institute. FastQC documentation—Per-tile sequencing quality. 2015. Available from: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3 Analysis Modules/12 Per Tile Sequence Quality.html [Last accessed: December 03, 2019].

10.

Baumeister

, Cocuzza

, Dam

, et al. A system for rapid DNA sequencing with fluorescent chain-terminating dideoxynucleotides. Science, 1987; 238(4825):336–341.

11.

Black

, Dueck

, Hayer

, et al. GSE50445: IVT-Seq reveals extreme bias in RNA-sequencing. 2012. Available from: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50445 [Last accessed: November 01, 2019].

12.

Bowers

, Jarosz

, Jones

, et al. Direct RNA sequencing. Nature, 2009; 461(7265):814–818.

13.

Campbell

, Efron

, Faghri

, et al. Big data: Astronomical or genomical? PLoS Biol, 2015; 13(7):1–11.

14.

Cao

, Chen

, Clark

, et al. Oligo (dt) primer generates a high frequency of truncated cDNAs through internal poly (a) priming during reverse transcription. Proc Natl Acad Sci USA, 2002; 99(9):6152–6156.

15.

Chen

, Chiang

, Hwang

, et al. Effects of gc bias in next-generation-sequencing data on de novo genome assembly. PLoS One, 2013; 8(4):e62856.

16.

Connell

, Dodd

, Heiner

, et al. Fluorescence detection in automated DNA sequence analysis. Nature, 1985; 321(6071):674–679.

17.

Cyranoski

China's bid to be a DNA superpower. Nature, 2016; 534(7608):462–463.

18.

Derge

, Feingold

, Gerhard

, et al. The completion of the mammalian gene collection (MGC). Genome Res, 2009; 19(12):2324–2333.

19.

Dickerson

, Liu

, Loraine

. Comparisons of computational methods for differential alternative splicing detection using RNA-Seq in plant systems. BMC Bioinform, 2014; 15(1):1.

20.

Dohm

, Lottaz

, Borodina

, et al. Substantial biases in ultra-short read data sets from high-throughput DNA sequencing. Nucleic Acids Res, 2008; 36(16):e105.

21.

Donaghey

, Pachter

, Rinn

, et al. Improving RNA-Seq expression estimates by correcting for fragment bias. Genome Biol, 2011; 12(3):R22.

22.

Edgar

Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res, 2002; 30(1):207–210.

23.

Gallagher

DNA project “to make UK world genetic research leader.” BBC; 2014. Available from: www.bbc.co.uk/news/health-28488313 [Last accessed: January 21, 2019].

24.

Garber

, Grabherr

, Guttman

, et al. Computational methods for transcriptome annotation and quantification using RNA-Seq. Nat Methods, 2011; 8(6):469–477.

25.

Gerstein

, Snyder

, Wang

RNA-Seq: A revolutionary tool for transcriptomics. Nat Rev Genet, 2009; 10(1):57–63.

26.

Gordon

, Hannon

. Fastx-toolkit. FASTQ/A short-reads preprocessing tools (unpublished). 2010. Available from: http://hannonlab.cshl.edu/fastx_toolkit [Last accessed: January 21, 2019].

27.

Hansen

, Brenner

, Dudoit

. Biases in Illumina transcriptome sequencing caused by random hexamer priming. Nucleic Acids Res, 2010; 38(12):e131.

28.

Harrison

, Memon

, Owen

, et al. Identifying the impact of G-quadruplexes on affymetrix 3’ arrays using cloud computing. J Integr Bioinform, 2010; 7(2):111.

29.

, Shendure

Next-generation DNA sequencing. Nat Biotechnol, 2008; 26(10):1135–1145.

30.

Kapranov

, Letovsky

, Lipson

, et al. Protocol dependence of sequencing-based gene expression measurements. PLoS One, 2011; 6(5):e19287.

31.

Kukurba

, Montgomery

. RNA sequencing and analysis. Cold Spring Harb Protoc, 2015; 2015(11):951–969.

32.

Lahens

, Kavakli

, Zhang

, et al. IVT-Seq reveals extreme bias in RNA sequencing. Genome Biol, 2014; 15(6):R86.

33.

Leinonen

, Shumway

, Sugawara

The sequence read archive. Nucleic Acids Res, 2011; 39(Database Issue):D19–D21.

34.

Mardis

ER.

Anticipating the 1,000 dollar genome. Genome Biol, 2006; 7(7):112.

35.

McCue

, Mortazavi

, Schaeffer

, et al. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods, 2008; 5(7):621–628.

36.

Meacham

, Boffelli

, Dhahbi

, et al. Identification and correction of systematic error in high-throughput sequence data. BMC Bioinform, 2011; 12(1):451.

37.

Metzker

ML.

Sequencing technologies—The next generation. Nat Rev Genet, 2010; 11(1):31–46.

38.

Nakazato

, Ohta

, Bono

. Experimental design-based functional mining and characterization of high-throughput sequencing data in the sequence read archive. PLoS One, 2013; 8(10):e77910.

39.

Naval-Sánchez

, Potier

, Haagen

, Sánchez

, et al. Comparative motif discovery combined with comparative transcriptomics yields accurate targetome and enhancer predictions. Genome Res, 2013; 23(1):74–88.

40.

Pachter

, Salzberg

, Trapnell

Tophat: Discovering splice junctions with RNA-Seq. Bioinformatics, 2009; 25(9):1105–1111.

41.

Park

, Shen

, Tokheim

, et al. Identifying differential alternative splicing events from RNA sequencing data using RNAseq-mats. In: Deep Sequencing Data Analysis. (Shomron N, ed.) Humana Press: Totowa, NJ, USA; 2013; pp. 171–179.

42.

Reuters. 1000 Genomes Project Releases Data from Pilot Projects on Path to Providing Database for 2,500 Human Genomes—Freely available data supporting next generation of human genetic research. 2015. Available from: www.reuters.com/article/us-usa-obama-precisionmedicine-idUSKBN0L313R20150130 [Last accessed: February 02, 2019].

43.

Risso

, Schwartz

, Sherlock

, et al. GC-content normalization for RNA-Seq data. BMC Bioinform, 2011; 12(1):480.

44.

Sakharkar

, Vincent

, Pandjassarame

Distributions of exons and introns in the human genome. In Silico Biology, 2004; 4(4):387–439.

45.

Sequence Read Archive. Overview of the Sequence Read Archive (SRA). 2017. Available from: https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi/ [Last accessed: January 3, 2019].

46.

Siva

UK gears up to decode 100 000 genomes from NHS patients. Lancet, 2015; 385(9963):103–104.

47.

Soper

, Steven

, et al. DNA Sequencing Using Fluorescence Detection. Topics in Fluorescence Spectroscopy. Springer Publishing, Boston, MA; 2003; pp. 1–66.

48.

Theurkauf

, Weng

, Zamore

, et al. Strand-specific libraries for high throughput RNA sequencing (RNA-Seq) prepared without poly (a) selection. Silence, 2012; 3(1):9.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.24 MB