ShortRead

Reading FASTQ files

The FASTQ file format is the standard way of representing raw (unaligned) next generation sequencing reads, particular for the Illumina platform. The format basically consists of 4 lines per read, with the lines containing

1. Read name (sometimes includes flowcell ID or other information).
2. Read nucleotides
3. Either empty or a repeat of line 1
4. Encoded read quality scores

Paired-end reads are stored in two separate files, where the reads are ordered the same (this is obviously fragile; what if reads are re-ordered in one file and not the other).

These files are read by readFastq() which produces an object of class ShortReadQ

fastqDir <- system.file("extdata", "E-MTAB-1147", package = "ShortRead")
fastqPath <- list.files(fastqDir, pattern = ".fastq.gz$", full = TRUE)[1]
reads <- readFastq(fastqPath)
reads

## class: ShortReadQ
## length: 20000 reads; width: 72 cycles

Here we directly point the function to the file path. A paradigm which is often used in Bioconductor is to first put the file path into an object which represents a specific file type and then read it; see

fqFile <- FastqFile(fastqPath)
fqFile

## class: FastqFile 
## path: /Library/Frameworks/R.framework/Versio.../ERR127302_1_subset.fastq.gz
## isOpen: FALSE

reads <- readFastq(fqFile)

This appears to make little sense in this situation, but for really big files it makes sense to access them in chunks, see below for a BAM file example.

The ShortReadQ class is very similar to a DNAStringSet but it has two sets of strings: one for the read nucleotides and one for the base qualities. They are accessed as

sread(reads)[1:2]

##   A DNAStringSet instance of length 2
##     width seq
## [1]    72 GTCTGCTGTATCTGTGTCGGCTGTCTCGCGG...CAATGAAGGCCTGGAATGTCACTACCCCCAG
## [2]    72 CTAGGGCAATCTTTGCAGCAATGAATGCCAA...GTGGCTTTTGAGGCCAGAGCAGACCTTCGGG

quality(reads)[1:2]

## class: FastqQuality
## quality:
##   A BStringSet instance of length 2
##     width seq
## [1]    72 HHHHHHHHHHHHHHHHHHHHEBDBB?B:BBG...FEFBDBD@DDECEE3>:?;@@@>?=BAB?##
## [2]    72 IIIIHIIIGIIIIIIIHIIIIEGBGHIIIIH...IHIIIHIIIIIGIIIEGIIGBGE@DDGGGIG

id(reads)[1:2]

##   A BStringSet instance of length 2
##     width seq
## [1]    53 ERR127302.8493430 HWI-EAS350_0441:1:34:16191:2123#0/1
## [2]    53 ERR127302.21406531 HWI-EAS350_0441:1:88:9330:2587#0/1

A word on quality scores

Note how the quality scores are listed as characters. You can convert them into standard 0-40 integer quality scores by

as(quality(reads), "matrix")[1:2,1:10]

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]   39   39   39   39   39   39   39   39   39    39
## [2,]   40   40   40   40   39   40   40   40   38    40

In this conversion, each letter is matched to an integer between 0 and 40. This matching is known as the “encoding” of the quality scores and there has been different ways to do this encoding. Unfortunately, it is not stored in the FASTQ file which encoding is used, so you have to know or guess the encoding. The ShortRead package does this for you.

These numbers are supposed to related to the probability that the reported base is different from the template fragment (ie. a sequence error). One should be aware that this probabilistic interpretation is not always true; methods such as “quality-remapping” helps to ensure this.

Reading alignment files

In the early days of next generation sequencing, there was no standardized alignment output format. different aligners produced different output file, including Bowtie and MAQ. Later on, the SAM / BAM format was introduced and this is now the standard alignment output. ShortRead contains tools for reading these older alignment formats through the readAligned() function (the type argument support options such as type="Bowtie" and type="MAQMap" and type="MAQMapShort").

The package has some very old support for parsing BAM files, but use Rsamtools and GenomicAlignments for this task instead.