Quality assessment

FASTQ files are often used for basic quality assessment, often to augment the purely technical QA that might be provided by the sequencing center with QA relevant to overall experimental design. A QA report is generated by creating a vector of paths to FASTQ files

fls <- dir("/path/to", "*fastq$", full=TRUE)

collecting statistics over the files

qaSummary <- qa(fls, type="fastq")

and creating and viewing a report

browseURL(report(qaSummary))

By default, the report is based on a sample of 1M reads.

These QA facilities are easily augmented by writing custom functions for reads sampled from files, or by exploiting the elements of the object returned from qa(), e.g., for an analysis of ArrayExpress experiment E-MTAB-1147:

qaSummary
## class: FastqQA(10)
## QA elements (access with qa[["elt"]]):
##   readCounts: data.frame(16 3)
##   baseCalls: data.frame(16 5)
##   readQualityScore: data.frame(8192 4)
##   baseQuality: data.frame(1504 3)
##   alignQuality: data.frame(16 3)
##   frequentSequences: data.frame(800 4)
##   sequenceDistribution: data.frame(1953 4)
##   perCycle: list(2)
##     baseCall: data.frame(5681 4)
##     quality: data.frame(44246 5)
##   perTile: list(2)
##     readCounts: data.frame(0 4)
##     medianReadQualityScore: data.frame(0 4)
##   adapterContamination: data.frame(16 1)

For instance, the count of reads in each lane is summarized in the readCounts element, and can be displayed as

head(qaSummary[["readCounts"]])
##                          read filter aligned
## ERR127302_1.fastq.gz 29741852     NA      NA
## ERR127302_2.fastq.gz 29741852     NA      NA
## ERR127303_1.fastq.gz 32665567     NA      NA
## ERR127303_2.fastq.gz 32665567     NA      NA
## ERR127304_1.fastq.gz 31876181     NA      NA
## ERR127304_2.fastq.gz 31876181     NA      NA
head(qaSummary[["baseCalls"]])
##                             A        C        G        T     N
## ERR127302_1.fastq.gz 16439860 19641395 19547421 16335620 35704
## ERR127302_2.fastq.gz 16238041 20020655 19608896 16060661 71747
## ERR127303_1.fastq.gz 16826258 19204659 19448727 16507994 12362
## ERR127303_2.fastq.gz 16426991 19822132 19374419 16324978 51480
## ERR127304_1.fastq.gz 16279217 19740457 19879137 16089405 11784
## ERR127304_2.fastq.gz 15984998 20297064 19812474 15853510 51954

The readCounts element contains a data frame with 1 row and 3 columns (these dimensions are indicated in the parenthetical annotation of readCounts in the output of qaSummary). The rows represent different lanes. The columns indicated the number of reads, the number of reads surviving the Solexa filtering criteria, and the number of reads aligned to the reference genome for the lane. The baseCalls element summarizes the base calls in the unfiltered reads.

The functions that produce the report tables and graphics are internal to the package. They can be accessed by calling ShortRead:::functionName where functionName is one of the functions listed below, organized by report section.

• Run Summary: .ppnCount, .df2a, .laneLbl, .plotReadQuality
• Read Distribution: .plotReadOccurrences, .freqSequences
• Cycle Specific: .plotCycleBaseCall, .plotCycleQuality
• Tile Performance: .atQuantile, .colorkeyNames, .plotTileLocalCoords, .tileGeometry, .plotTileCounts, .plotTileQualityScore
• Alignment: .plotAlignQuality
• Multiple Alignment: .plotMultipleAlignmentCount
• Depth of Coverage: .plotDepthOfCoverage
• Adapter Contamination: .ppnCount

Filtering and trimming

It is straightforward to create filters to eliminate reads or to trim reads based on diverse characteristics. The basic structure is to open a FASTQ file, iterate through chunks of the file, performing filtering or trimming steps, and appending the filtered data to a new file.

myFilterAndTrim <-
    function(fl, destination=sprintf("%s_subset", fl))
{
    ## open input stream
    stream <- open(FastqStreamer(fl))
    on.exit(close(stream))
    repeat {
        ## input chunk
        fq <- yield(stream)
        if (length(fq) == 0)
            break
        ## trim and filter, e.g., reads cannot contain 'N'...
        fq <- fq[nFilter()(fq)]  # see ?srFilter for pre-defined filters
        ## trim as soon as 2 of 5 nucleotides has quality encoding less
        ## than "4" (phred score 20)
        fq <- trimTailw(fq, 2, "4", 2)
        ## drop reads that are less than 36nt
        fq <- fq[width(fq) >= 36]
        ## append to destination
        writeFastq(fq, destination, "a")
    }
}

This is memory efficient and flexible. Care must be taken to coordinate pairs of FASTQ files representing paired-end reads, to preserve order.

Data I/O

ShortRead provides a variety of methods to read data into R, in addition to readAligned.

## location of file
exptPath <- system.file("extdata", package="ShortRead")
sp <- SolexaPath(exptPath)
pattern <- "s_2_export.txt"
fl <- file.path(analysisPath(sp), pattern)
strsplit(readLines(fl, n=1), "\t")
## [[1]]
##  [1] "HWI-EAS88"                           "3"                                  
##  [3] "2"                                   "1"                                  
##  [5] "451"                                 "945"                                
##  [7] ""                                    ""                                   
##  [9] "CCAGAGCCCCCCGCTCACTCCTGAACCAGTCTCTC" "YQMIMIMMLMMIGIGMFICMFFFIMMHIIHAAGAH"
## [11] "NM"                                  ""                                   
## [13] ""                                    ""                                   
## [15] ""                                    ""                                   
## [17] ""                                    ""                                   
## [19] ""                                    ""                                   
## [21] ""                                    "N"
length(readLines(fl))
## [1] 1000

colClasses <- rep(list(NULL), 21)
colClasses[9:10] <- c("DNAString", "BString")
names(colClasses)[9:10] <- c("read", "quality")

cols <- readXStringColumns(analysisPath(sp), pattern, colClasses)
cols
## $read
## DNAStringSet object of length 1000:
##        width seq
##    [1]    35 CCAGAGCCCCCCGCTCACTCCTGAACCAGTCTCTC
##    [2]    35 AGCCTCCCTCTTTCTGAATATACGGCAGAGCTGTT
##    [3]    35 ACCAAAAACACCACATACACGAGCAACACACGTAC
##    [4]    35 AATCGGAAGAGCTCGTATGCCGGCTTCTGCTTGGA
##    [5]    35 AAAGATAAACTCTAGGCCACCTCCTCCTTCTTCTA
##    ...   ... ...
##  [996]    35 GTGGCAGCGGTGAGGCGGCGGGGGGGGGTTGTTTG
##  [997]    35 GTCGGAGGTCAGCAAGCTGTAGTCGGTGTAAAGCT
##  [998]    35 GTCATAAATTGGACAGTGTGGCTCCAGTATTCTCA
##  [999]    35 ATCTACATTAAGGTCAATTACAATGATAAATAAAA
## [1000]    35 TTCTCAGCCATTCAGTATTCCTCAGGTGAAAATTC
## 
## $quality
## BStringSet object of length 1000:
##        width seq
##    [1]    35 YQMIMIMMLMMIGIGMFICMFFFIMMHIIHAAGAH
##    [2]    35 ZXZUYXZQYYXUZXYZYYZZXXZZIMFHXQSUPPO
##    [3]    35 LGDHLILLLLLLLIGFLLALDIFDILLHFIAECAE
##    [4]    35 JJYYIYVSYYYYYYYYSDYYWVUYYNNVSVQQELQ
##    [5]    35 LLLILIIIDLLHLLLLLLLLLLLALLLLHLLLLEL
##    ...   ... ...
##  [996]    35 ZZZZZZZYZZYUYZYUYZKYUDUZIYYODJGUGAA
##  [997]    35 ZZZZZZZZZZZZZZZZZZYZZYXXZYSSXXUUHHQ
##  [998]    35 ZZZZZZZZZZZZZZZYZZZZYZZZZYZZXZUUUUS
##  [999]    35 ZZZZZZZZZZZYXZYZYZZYZYZZXKZSYXUUNUN
## [1000]    35 ZZZZZZZZZZZZZZYZZZZZZZZYYSYSZXUUUUU

object.size(cols$read)
## 51032 bytes
object.size(as.character(cols$read))
## 102128 bytes

Sorting

Short reads can be sorted using srsort, or the permutation required to bring the short read into lexicographic order can be determined using srorder. These functions are different from sort and order because the result is independent of the locale, and they operate quickly on DNAStringSet and BStringSet objects.

The function srduplicated identifies duplicate reads. This function returns a logical vector, similar to duplicated. The negation of the result from srduplicated is useful to create a collection of unique reads. An experimental scenario where this might be useful is when the sample preparation involved PCR. In this case, replicate reads may be due to artifacts of sample preparation, rather than differential representation of sequence in the sample prior to PCR.

Summarizing read occurrence

The tables function summarizes read occurrences, for instance,

tbls <- tables(fq)
names(tbls)
## [1] "top"          "distribution"
tbls$top[1:5]
## CTATTCTCTACAAACCACAAAGACATTGGAACACTATACCTATTATTCGGCGCATGAGCTGGAGTCCTAGGC 
##                                                                        7 
## GTTTGGTCTAGGGTGTAGCCTGAGAATAGGGGAAATCAGTGAATGAAGCCTCCTATGATGGCAAATACAGCT 
##                                                                        7 
## CGATAACGTTGTAGATGTGGTCGTTACCTAGAAGGTTGCCTGGCTGGCCCAGCTCGGCTCGAATAAGGAGGC 
##                                                                        6 
## CTAGCATTTACCATCTCACTTCTAGGAATACTAGTATATCGCTCACACCTCATATCCTCCCTACTATGCCTA 
##                                                                        6 
## CACGAGCATATTTCACCTCCGCTACCATAATCATCGCTATCCCCACCGGCGTCAAAGTATTTAGCTGACTCG 
##                                                                        5
head(tbls$distribution)
##   nOccurrences nReads
## 1            1  19291
## 2            2    247
## 3            3     34
## 4            4     18
## 5            5      3
## 6            6      2

The top component returned by tables is a list tallying the most commonly occurring sequences in the short reads. Knowledgeable readers will recognize the top-occurring read as a close match to one of the manufacturer adapters.

The distribution component returned by tables is a data frame that summarizes how many reads (e.g., 19291) are represented exactly 1 times.

Finding near matches to short sequences

Facilities exist for finding reads that are near matches to specific sequences, e.g., manufacturer adapter or primer sequences. srdistance reports the edit distance between each read and a reference sequence. srdistance is implemented to work efficiently for reference sequences whose length is of the same order as the reads themselves (10’s to 100’s of bases). To find reads close to the most common read in the example above, one might say

dist <- srdistance(sread(fq), names(tbls$top)[1])[[1]]
table(dist)[1:10]
## dist
##  0  4  6 10 14 18 20 21 31 32 
##  7  1  3  1  3  1  4  1  3 11

‘Near’ matches can be filtered, e.g.,

fqSubset <- fq[dist>4]

A different strategy can be used to tally or eliminate reads that consist predominantly of a single nucleotide. alphabetFrequency calculates the frequency of each nucleotide (in DNA strings) or letter (for other string sets) in each read. Thus one could identify and eliminate reads with more than 30 adenine nucleotides with

countA <- alphabetFrequency(sread(fq))[,"A"]
fqNoPolyA <- fq[countA < 30]

alphabetFrequency, which simply counts nucleotides, is much faster than srdistance, which performs full pairwise alignment of each read to the subject.

Users wanting to use R for whole-genome alignments or more flexible pairwise alignment are encouraged to investigate the Biostrings and pwalign packages, especially the PDict class and matchPDict and pairwiseAlignment functions.