@@ -1,19 +1,24 @@

 Package: ShortRead

 Type: Package

-Title: Classes and methods for high-throughput short-read sequencing data.

-Version: 1.21.1

+Title: FASTQ input and manipulation

+Version: 1.21.2

 Author: Martin Morgan, Michael Lawrence, Simon Anders

 Maintainer: Bioconductor Package Maintainer <maintainer@bioconductor.org>

-Description: Base classes, functions, and methods for representation of

-             high-throughput, short-read sequencing data.

+Description: This package implements sampling, iteration, and input of

+    FASTQ files. The package includes functions for filtering and

+    trimming reads, and for generating a quality assessment

+    report. Data are represented as DNAStringSet-derived objects, and

+    easily manipulated for a diversity of purposes.  The package also

+    contains legacy support for early single-end, ungapped alignment

+    formats.

 License: Artistic-2.0

 LazyLoad: yes

-Depends: methods, BiocGenerics (>= 0.1.0), IRanges (>= 1.19.34),

-         GenomicRanges (>= 1.13.43), Biostrings (>= 2.29.18),

-         lattice, Rsamtools (>= 1.13.1)

-Imports: Biobase, hwriter, zlibbioc, latticeExtra

+Depends: BiocGenerics (>= 0.1.0), IRanges (>= 1.19.34),

+    GenomicRanges (>= 1.13.43), Biostrings (>= 2.29.18),

+    Rsamtools (>= 1.13.1)

+Imports: Biobase, hwriter, methods, zlibbioc, lattice, latticeExtra

 Enhances: Rmpi, parallel

-Suggests: biomaRt, RUnit, GenomicFeatures, yeastNagalakshmi

+Suggests: BiocStyle, RUnit, biomaRt, GenomicFeatures, yeastNagalakshmi

 LinkingTo: IRanges, XVector, Biostrings

 biocViews: DataImport, Sequencing, HighThroughputSequencing, 

            QualityControl

@@ -71,14 +71,14 @@

 setMethod(qa, "ShortReadQ", .qa_ShortReadQ)

 .qa_fastq_lane <-

-    function(dirPath, pattern, ..., sample=TRUE, type="fastq", 

+    function(dirPath, ..., sample=TRUE, type="fastq", 

         verbose=FALSE)

 {

     if (verbose)

         message("qa 'fastq' pattern:", pattern)

     if (sample) {

-        samp <- FastqSampler(file.path(dirPath, pattern), ...)

-        qa <- qa(yield(samp), pattern, ..., verbose=verbose)

+        samp <- FastqSampler(dirPath, ...)

+        qa <- qa(yield(samp), basename(dirPath), ..., verbose=verbose)

         close(samp)

         elts <- .srlist(qa)

         elts$readCounts$read <- samp$status()[["total"]]

@@ -95,8 +95,7 @@ setMethod(qa, "ShortReadQ", .qa_ShortReadQ)

 {

     fls <- .file_names(dirPath, pattern)

     lst <-

-        srapply(basename(fls), .qa_fastq_lane,

-                dirPath=dirPath, type=type, ...,

+        srapply(fls, .qa_fastq_lane, type=type, ...,

                 reduce=.reduce(1), verbose=verbose, 

                 USE.NAMES=TRUE)

     lst <- do.call(rbind, lst)

@@ -108,7 +108,7 @@ uniqueFilter <-

 }

 occurrenceFilter <-

-    function(min=1L, max=1L, withSread=c(TRUE, FALSE, NA),

+    function(min=1L, max=1L, withSread=c(NA, TRUE, FALSE),

              duplicates=c("head", "tail", "sample", "none"),

              .name=.occurrenceName(min, max, withSread,

                  duplicates))

old mode 100644

new mode 100755

new file mode 100644

Binary files /dev/null and b/inst/extdata/E-MTAB-1147/ERR127302_1_subset.fastq.gz differ

new file mode 100644

Binary files /dev/null and b/inst/extdata/E-MTAB-1147/ERR127302_2_subset.fastq.gz differ

@@ -4,7 +4,7 @@

 \alias{append,AlignedDataFrame,AlignedDataFrame-method}

 \title{

-  "AlignedDataFrame" representing alignment annotations as a data frame

+  (Legacy) "AlignedDataFrame" representing alignment annotations as a data frame

 }

 \description{

@@ -1,7 +1,7 @@

 \name{AlignedDataFrame}

 \alias{AlignedDataFrame}

-\title{AlignedDataFrame constructor}

+\title{(Legacy) AlignedDataFrame constructor}

 \description{

   Construct an \code{AlignedDataFrame} from a data frame and its

@@ -27,7 +27,7 @@

 \alias{detail,AlignedRead-method}

 \alias{show,AlignedRead-method}

-\title{"AlignedRead" class for aligned short reads}

+\title{(Legacy) "AlignedRead" class for aligned short reads}

 \description{

@@ -1,7 +1,7 @@

 \name{AlignedRead}

 \alias{AlignedRead}

-\title{Construct objects of class "AlignedRead"}

+\title{(Legacy) Construct objects of class "AlignedRead"}

 \description{

   This function constructs objects of

@@ -4,7 +4,7 @@

 \alias{report,BowtieQA-method}

 \alias{report_html,BowtieQA-method}

-\title{Quality assessment summaries from Bowtie files}

+\title{(Legacy) Quality assessment summaries from Bowtie files}

 \description{

@@ -7,7 +7,7 @@

 \alias{show,ExperimentPath-method}

 \alias{detail,ExperimentPath-method}

-\title{"ExperimentPath" class representing a file hierarchy of data

+\title{(Legacy) "ExperimentPath" class representing a file hierarchy of data

   files}

 \description{

@@ -2,7 +2,6 @@

 \docType{class}

 % Class:

-\alias{class:GappedReads}

 \alias{GappedReads-class}

 \alias{GappedReads}

@@ -22,7 +21,7 @@

 \alias{narrow,GappedReads-method}

-\title{GappedReads objects}

+\title{(Legacy) GappedReads objects}

 \description{

   The GappedReads class extends the \link[GenomicRanges]{GAlignments}

@@ -21,7 +21,7 @@

 \alias{show,Intensity-method}

 \alias{show,IntensityMeasure-method}

-\title{"Intensity", "IntensityInfo", and "IntensityMeasure" base

+\title{(Legacy) "Intensity", "IntensityInfo", and "IntensityMeasure" base

   classes for short read image intensities}

 \description{

@@ -5,7 +5,7 @@

 \alias{report,MAQMapQA-method}

 \alias{report_html,MAQMapQA-method}

-\title{Quality assessment summaries from MAQ map files}

+\title{(Legacy) Quality assessment summaries from MAQ map files}

 \description{

@@ -23,7 +23,7 @@

 %

 \alias{RochePath}

-\title{"RochePath" class representing a Roche (454) experiment location}

+\title{(Legacy) "RochePath" class representing a Roche (454) experiment location}

 \description{

@@ -3,7 +3,7 @@

 \alias{RocheSet-class}

 \alias{RocheSet}

-\title{Roche (454) experiment-wide data container}

+\title{(Legacy) Roche (454) experiment-wide data container}

 \description{

   This class is meant to coordinate all data in a Roche (454)

@@ -2,7 +2,7 @@

 \docType{class}

 \alias{RtaIntensity-class}

-\title{Class "RtaIntensity"}

+\title{(Legacy) Class "RtaIntensity"}

 \description{

   Subclass of \code{\linkS4class{Intensity}} for representing image

@@ -1,7 +1,7 @@

 \name{RtaIntensity}

 \alias{RtaIntensity}

-\title{Construct objects of class "RtaIntensity"}

+\title{(Legacy) Construct objects of class "RtaIntensity"}

 \description{

@@ -10,7 +10,8 @@

 \alias{show,SRSet-method}

 \alias{detail,SRSet-method}

-\title{A base class for Roche experiment-wide data}

+\title{(Legacy) A base class for Roche experiment-wide data}

+

 \description{

   This class coordinates phenotype (sample) and sequence data, primarily

   as used on the Roche platform.

@@ -8,13 +8,17 @@

 \docType{package}

 \title{

-  Base classes and methods for high-throughput short-read sequencing data.

+  FASTQ input and manipulation.

 }

 \description{

-  Base classes, functions, and methods for representation of

-  high-throughput, short-read sequencing data.

+  This package implements sampling, iteration, and input of FASTQ

+  files. The package includes functions for filtering and trimming

+  reads, and for generating a quality assessment report. Data are

+  represented as DNAStringSet-derived objects, and easily manipulated

+  for a diversity of purposes.  The package also contains legacy support

+  for early single-end, ungapped alignment formats.

 }

 \details{

@@ -8,7 +8,8 @@

 \alias{report_html,SolexaRealignQA-method}

 \alias{show,SolexaExportQA-method}

-\title{Quality assessment summaries from Solexa export and realign files}

+\title{(Legacy) Quality assessment summaries from Solexa export and

+  realign files}

 \description{

@@ -2,7 +2,8 @@

 \alias{SolexaIntensity}

 \alias{SolexaIntensityInfo}

-\title{Construct objects of class "SolexaIntensity" and "SolexaIntensityInfo"}

+\title{(Legacy) Construct objects of class "SolexaIntensity" and

+"SolexaIntensityInfo"}

 \description{

@@ -18,7 +18,7 @@

 \alias{readBaseQuality,SolexaPath-method}

 \alias{readAligned,SolexaPath-method}

-\title{"SolexaPath" class representing a standard output file hierarchy}

+\title{(Legacy) "SolexaPath" class representing a standard output file hierarchy}

 \description{

   Solexa produces a hierarchy of output files. The content of the

@@ -12,7 +12,7 @@

 % transforming methods

 \alias{readAligned,SolexaSet-method}

-\title{"SolexaSet" coordinating Solexa output locations with sample annotations}

+\title{(Legacy) "SolexaSet" coordinating Solexa output locations with sample annotations}

 \description{

   This class coordinates the file hierarchy produced by the Solexa

@@ -72,16 +72,16 @@ x = numeric(1000)

 x[sample(1000, 100)] <- abs(rnorm(100))

 df <- data.frame(x = c(x, -x), pos = seq(1, 1e5, length.out=1000),

                  group = rep(c("positive", "negative"), each=1000))

-cv <- xyplot(x ~ pos, df, group=group, type="s",

-             col=col, main="yeast chrI:1 - 2e5",

-             ylab="Coverage", xlab="Coordinate",

-             scales=list(y=list(tck=c(1,0)),

-                         x=list(rot=45, tck=c(1,0), tick.number=20)),

-             panel=function(...) {

-                     panel.xyplot(...)

-                     panel.grid(h=-1, v=20)

-                     panel.abline(a=0, b=0, col="grey")

-             })

+cv <- lattice::xyplot(x ~ pos, df, group=group, type="s",

+    col=col, main="yeast chrI:1 - 2e5",

+    ylab="Coverage", xlab="Coordinate",

+    scales=list(y=list(tck=c(1,0)),

+                x=list(rot=45, tck=c(1,0), tick.number=20)),

+    panel=function(...) {

+            lattice::panel.xyplot(...)

+            lattice::panel.grid(h=-1, v=20)

+            lattice::panel.abline(a=0, b=0, col="grey")

+    })

 s <- SpTrellis(cv)

 s

 zi(s)

@@ -24,7 +24,7 @@

 \alias{laneDescription}

 \alias{laneNames}

-\title{Accessors for ShortRead classes}

+\title{(Legacy) Accessors for ShortRead classes}

 \description{

@@ -3,7 +3,7 @@

 \alias{readAligned}

 \alias{readAligned,character-method}

-\title{Read aligned reads and their quality scores into R representations}

+\title{(Legacy) Read aligned reads and their quality scores into R representations}

 \description{

@@ -3,7 +3,7 @@

 \alias{readBaseQuality}

 \alias{readBaseQuality,character-method}

-\title{Read short reads and their quality scores into R representations}

+\title{(Legacy) Read short reads and their quality scores into R representations}

 \description{

@@ -1,7 +1,7 @@

 \name{readBfaToc}

 \alias{readBfaToc}

-\title{Get a list of the sequences in a Maq .bfa file}

+\title{(Legacy) Get a list of the sequences in a Maq .bfa file}

 \description{

 As \code{\link{coverage}} needs to know the lengths of the reference sequences,

@@ -3,7 +3,7 @@

 \alias{readIntensities}

 \alias{readIntensities,character-method}

-\title{Read Illumina image intensity files}

+\title{(Legacy) Read Illumina image intensity files}

 \description{

@@ -3,7 +3,7 @@

 \alias{readPrb}

 \alias{readPrb,character-method}

-\title{Read Solexa prb files as fastq-style quality scores}

+\title{(Legacy) Read Solexa prb files as fastq-style quality scores}

 \description{

@@ -3,7 +3,7 @@

 \alias{readQseq}

 \alias{readQseq,character-method}

-\title{Read Solexa qseq files as fastq-style quality scores}

+\title{(Legacy) Read Solexa qseq files as fastq-style quality scores}

 \description{

@@ -35,12 +35,8 @@ compose(filt, ..., .name)

 idFilter(regex=character(0), fixed=FALSE, exclude=FALSE,

          .name="idFilter")

-chromosomeFilter(regex=character(0), fixed=FALSE, exclude=FALSE,

-                 .name="ChromosomeFilter")

-positionFilter(min=-Inf, max=Inf, .name="PositionFilter")

-strandFilter(strandLevels=character(0), .name="StrandFilter")

 occurrenceFilter(min=1L, max=1L,

-                 withSread=c(TRUE, FALSE, NA),

+                 withSread=c(NA, TRUE, FALSE),

                  duplicates=c("head", "tail", "sample", "none"),

                  .name=.occurrenceName(min, max, withSread,

                                        duplicates))

@@ -50,6 +46,15 @@ polynFilter(threshold=0L, nuc=c("A", "C", "T", "G", "other"),

 dustyFilter(threshold=Inf, batchSize=NA, .name="DustyFilter")

 srdistanceFilter(subject=character(0), threshold=0L,

                  .name="SRDistanceFilter")

+

+##

+## legacy filters for ungapped alignments

+##

+

+chromosomeFilter(regex=character(0), fixed=FALSE, exclude=FALSE,

+                 .name="ChromosomeFilter")

+positionFilter(min=-Inf, max=Inf, .name="PositionFilter")

+strandFilter(strandLevels=character(0), .name="StrandFilter")

 alignQualityFilter(threshold=0L, .name="AlignQualityFilter")

 alignDataFilter(expr=expression(), .name="AlignDataFilter")

 }

@@ -174,7 +179,7 @@ alignDataFilter(expr=expression(), .name="AlignDataFilter")

   treated: \code{TRUE} to include the sread, chromosome, strand, and

   position when determining occurrence, \code{FALSE} to include

   chromosome, strand, and position, and \code{NA} to include only

-  sread. The default is \code{withSread=TRUE}. \code{duplicates}

+  sread. The default is \code{withSread=NA}. \code{duplicates}

   determines how reads with more than \code{max} reads are

   treated. \code{head} selects the first \code{max} reads of each set of

   duplicates, \code{tail} the last \code{max} reads, and \code{sample} a

@@ -63,7 +63,7 @@ showMethods("tables")

 sp <- SolexaPath(system.file("extdata", package="ShortRead"))

 aln <- readAligned(sp)

 tables(sread(aln), n=6)

-xyplot(log10(nReads)~log10(nOccurrences),

+lattice::xyplot(log10(nReads)~log10(nOccurrences),

        tables(sread(aln))$distribution)

 }

 \keyword{manip}

@@ -1,394 +1,170 @@

 %\VignetteIndexEntry{An introduction to ShortRead}

-%\VignetteDepends{}

+%\VignetteDepends{BiocStyle}

 %\VignetteKeywords{Short read, I/0, quality assessment}

 %\VignettePackage{ShortRead}

 \documentclass[]{article}

-\usepackage{times}

-\usepackage{hyperref}

-

-\newcommand{\Rfunction}[1]{{\texttt{#1}}}

-\newcommand{\Robject}[1]{{\texttt{#1}}}

-\newcommand{\Rpackage}[1]{{\textit{#1}}}

-\newcommand{\Rmethod}[1]{{\texttt{#1}}}

-\newcommand{\Rfunarg}[1]{{\texttt{#1}}}

-\newcommand{\Rclass}[1]{{\textit{#1}}}

-\newcommand{\Rcode}[1]{{\texttt{#1}}}

-

-\newcommand{\software}[1]{\textsf{#1}}

-\newcommand{\R}{\software{R}}

-\newcommand{\ShortRead}{\Rpackage{ShortRead}}

+<<style, eval=TRUE, echo=FALSE, results=tex>>=

+BiocStyle::latex()

+@ 

-\newcommand{\ELAND}{\software{ELAND}}

-\newcommand{\MAQ}{\software{MAQ}}

-\newcommand{\Bowtie}{\software{Bowtie}}

+\newcommand{\ShortRead}{\Biocpkg{ShortRead}}

 \title{An Introduction to \Rpackage{ShortRead}}

 \author{Martin Morgan}

-\date{Modified: 28 September 2010. Compiled: \today}

+\date{Modified: 21 October, 2013. Compiled: \today}

 \begin{document}

 \maketitle

-<<options,echo=FALSE>>=

-options(width=60)

-@ 

-

 <<preliminaries>>=

 library("ShortRead")

 @ 

-The \Rpackage{ShortRead} package aims to provide key functionality for

-input, quality assurance, and basic manipulation of `short read' DNA

-sequences such as those produced by Solexa, 454, and related

-technologies, including flexible import of common short read data

-formats. This vignette introduces key functionality.

-

-Support is most fully developed for Solexa; contributions from the

-community are welcome.

-

-\section{A first workflow}

-

-This section walks through a simple work flow. It outlines the

-hierarchy of files produced by Solexa. It then illustrates a common

-way for reading short read data into \R{}.

-

-\subsection{\Rclass{SolexaPath}: navigating Solexa output}

-

-\Rclass{SolexaPath} provides functionality to navigate files produced

-by Solexa Genome Analyzer pipeline software. A typical way to start a

-\ShortRead{} session is to point to the root of the output file

-hierarchy. The \ShortRead{} package includes a very small subset of

-files emulating this hierarchy. The root is found at

-<<SolexaPath-root>>=

-exptPath <- system.file("extdata", package="ShortRead")

-@ 

-%% 

-Usually \Rcode{exptPath} would be a location on the users' file system. Key

-components of the hierarchy are parsed into \R{} with

-<<SolexaPat>>=

-sp <- SolexaPath(exptPath)

-sp

-@ 

-%% 

-\Rfunction{SolexaPath} scans the directory hierarchy to identifying

-useful directories. For instance, image intensity files are in the

-`Firecrest' directory, while summary and alignment files are in the

-analysis directory

-<<firecrest>>=

-imageAnalysisPath(sp)

-analysisPath(sp)

-@ 

-%% 

-Most functionality in \ShortRead{} uses \Rcode{baseCallPath} or

-\Rcode{analysisPath}. Solexa documentation provides details of file

-content. \Rfunction{SolexaPath} accepts additional arguments that

-allow individual file paths to be specified.

- 

-Many functions for Solexa data input `know' where appropriate files

-are located. Specifying \Rcode{sp} is often sufficient for identifying

-the desired directory path. Examples of this are illustrated below,

-with for instance \Rfunction{readAligned} and \Rfunction{readFastq}.

-

-Displaying an object, e.g., \Robject{sp}, provides hints at how to

-access information in the object, e.g., \Rfunction{analysisPath}. This

-is a convention in \ShortRead{}.

-

-\subsection{\Rfunction{readAligned}: reading aligned data into \R{}}

-

-Solexa \texttt{s\_N\_export.txt} files (\texttt{\_N\_} is a

-placeholder for the lane identifier) represent one place to start

-working the short read data in \R{}. These files result from running

-ANALYSIS eland\_extended in the Solexa Genome Analyzer. The files

-contain information on all reads, including alignment information for

-those reads successfully aligned to the genome.  \ShortRead{} parses

-additional alignment files, including \MAQ{} binary and text

-(\texttt{mapview}) files and \Bowtie{} text files; consult the help page

-for \Rfunction{readAligned} for details. \ShortRead{} flexibly parses

-many other Solexa files; aligned reads represent just one entry point.

-

-To read a single \texttt{s\_N\_export.txt} file into \R{}, for instance

-from lane 2, use the command

-<<readAligned-simple>>=

-aln <- readAligned(sp, "s_2_export.txt")

-aln

-@ 

-%% 

-This illustrates the convention used for identifying files for input

-into \R{} and used by \ShortRead{}. The function takes a directory

-path and a pattern (as a regular expression, similar to the \R{}

-function \Rfunction{list.files}) of file names to match in the

-directory. Usually, all files matching the pattern are read into a

-\emph{single} \R{} object; this behavior is desirable for several of

-the input functions in \ShortRead{}. In the present case the usual

-expectation is that a single \texttt{s\_N\_export.txt} file will be

-read into a single \R{} object, so the \Rfunarg{pattern} argument will

-identify a single file.

-

-\subsubsection{Input of other aligned read files}

-

-\ELAND{} software provides access to much interesting data, in

-addition to alignments, but if the interest is in aligned reads then

-input may come from any of a number of different software

-packages. Many of these alignments can be input with

-\Rpackage{ShortRead}.

-

-\Bowtie{} is a very fast aligner, taking a few tens of minutes to

-align entire lanes of reads to reference genomes. Use

-\Rfunction{readAligned} with the \Rfunarg{type="Bowtie"} argument to

-input alignments.  Reading \Bowtie{} output using

-\Rfunction{readAligned} produces the same class of object as reading

-\ELAND{} output.  Like \ELAND{}, \Bowtie{} provides information on

-short read, quality, chromosome, position, and strand; there is no

-information on alignment quality avaiable from \Bowtie{}. \ELAND{} and

-\Bowtie{} provide very different auxiliary information.  Consult the

-\Rfunction{readAligned} and help page \Bowtie{} manual for additional

-detail.

-

-\MAQ{} is another poplar aligner. \Rpackage{ShortRead} can input

-\MAQ{} binary or text formats (see the arguments

-\Rfunarg{type="MAQMapShort"}, \Rfunarg{"MAQMap"}, and

-\Rfunarg{"MAQMapview"}). As with \Bowtie{}, \MAQ{} provides essential

-information about reads and their aligments, plus additional

-information that differs somewhat from the additional information

-provided by \ELAND{}.

-

-Alignment information may come in a variety of different text-based

-formats. Not all of these will be supported by

-\Rpackage{ShortRead}. There are a number of tools available to input

-this into \R{}. 

-

-A basic strategy is involves two passes over the data, followed by

-synthesis of results into an \Rclass{AlignedRead} object.  First,

-input alignment data using functions such as

-\Rfunction{read.table}. Use the \Rfunarg{colClasses} argument to

-`mask-out' (i.e., avoid importing) DNA and quality sequences. Next,

-use \Rfunction{readXStringColumns} or \Rfunction{readFastq} to import

-the short read and quality information. Finally, use the alignment

-data and reads as arugments to the \Rfunction{AlignedRead} function to

-synthesize the input. The following illustrates use of

-\Rfunction{readXStringColumns} and \Rfunction{readFastq}. These

-functions receive further attention below.

-

-\subsubsection{Cautions}

-

-There are several confusing areas of input. (1) Some alignment

-programs and genome resources start numbering nucleotides of the

-subject sequence at 0, whereas others start at 1. (2) Some alignment

-programs report matches on the minus strand in terms of the

-`left-most' position of the read (i.e., the location of the 3' end of

-the aligned read), whereas other report `five-prime''matches (i.e., in

-terms of the 5' end of the read), regardless of whether the alignment

-is on the plus or minus strand. (3) Some alignment programs reverse

-complement the sequence of reads aligned to the minus strand. (4) Base

-qualities are sometimes encoded as character strings, but the encoding

-differs between `fastq' and `solexa fastq'. It seems that all

-combinations of these choices are common `in the wild'.

-

-The help page for \Rfunction{readAligned} attempts to be explicit

-about how reads are formatted. Briefly:

-\begin{itemize}

-\item Subject sequence nucleotides are numbered starting at 1, rather

-  than zero. \Rfunction{readAligned} adjusts the coordinate system of

-  input reads if necessary (e.g., reading \MAQ{} alignments).

-\item Alignments on the minus strand are reported in `left-most'

-  coordinates systems.

-\item \ELAND{} and \Bowtie{} alignments on the minus strand are not reverse

-  complemented.

-\item Character-encoded base quality scores are intrepreted as the

-  default for the software package being parsed, e.g., as `Solexa

-  fastq' for \ELAND{}. The object returned by \Rfunction{quality}

-  applied to an \Rclass{AlignedRead} object is either

-  \Rclass{FastqQuality} or \Rclass{SFastqQuality}.

-\end{itemize}

-Alignment programs sometimes offer the opportunity to custommize

-output; such customization needs to be accomodated when reads are

-input using \Rpackage{ShortRead}.

-

-\subsubsection{Filtering input}

-

-Downstream analysis may often want to use a well-defined subset of

-reads. These can be selected with the \Rfunarg{filter} argument of

-\Rfunction{readAligned}. There are built-in filters, for instance to

-remove all reads containing an \texttt{N} nucleotide, to select just

-those reads that map to the genome file \texttt{chr5.fa}, to select

-reads on the \texttt{+} strand, or to `level the playing field' by

-selecting only a single read for any chromosome, position and strand:

-<<filter-egs>>=

-nfilt <- nFilter()

-cfilt <- chromosomeFilter('chr5.fa')

-sfilt <- strandFilter("+")

-ofilt <- occurrenceFilter(withSread=FALSE)

+The \Rpackage{ShortRead} package provides functionality for working

+with FASTQ files from high throughput sequence analysis. The package

+also contains functions for legacy (single-end, ungapped) aligned

+reads; for working with BAM files, please see the \Biocpkg{Rsamtools},

+\Biocpkg{GenomicRanges}, and related packages.

+

+\section{Sample data}

+

+Sample FASTQ data are derived from ArrayExpress record

+\href{http://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-1147/}{E-MTAB-1147}.

+Paired-end FASTQ files were retrieved and then sampled to 20,000

+records with

+<<sample, eval=FALSE>>=

+sampler <- FastqSampler('E-MTAB-1147/fastq/ERR127302_1.fastq.gz', 20000)

+set.seed(123); ERR127302_1 <- yield(sampler)

+sampler <- FastqSampler('E-MTAB-1147/fastq/ERR127302_2.fastq.gz', 20000)

+set.seed(123); ERR127302_2 <- yield(sampler)

+@ 

+

+\section{Functionality}

+

+Functionality is summarized in Table~\ref{tab:fastq}. 

+

+\begin{table}

+  \centering

+  \begin{tabular}{lll}

+    \hline

+    \multicolumn{3}{l}{Input} \\

+    & \Rfunction{FastqStreamer} & Iterate through FASTQ files in chunks \\

+    & \Rfunction{FastqSampler} & Draw random samples from FASTQ files \\

+    & \Rfunction{readFastq} & Read an entire FASTQ file into memory \\

+    & \Rfunction{writeFastq} & Write FASTQ objects to a connection (file) \\

+    \multicolumn{3}{l}{Sequence and quality summary} \\

+    & \Rfunction{alphabetFrequency} & Nucleotide or quality score use per read\\

+    & \Rfunction{alphabetByCycle} & Nucleotide or quality score use by cycle\\

+    & \Rfunction{alphabetScore} & Whole-read quality summary\\

+    & \Rfunction{encoding} & Character / `phred' score mapping \\

+    \multicolumn{3}{l}{Quality assessment} \\

+    & \Rfunction{qa} & Visit FASTQ files to collect QA statistics \\

+    & \Rfunction{report} & Generate a quality assessment report \\

+    \multicolumn{3}{l}{Filtering and trimming} \\

+    & \Rfunction{srFilter} & Pre-defined and bespoke filters \\

+    & \Rfunction{trimTails}, etc. & Trim low-quality nucleotides \\

+    & \Rfunction{narrow} & Remove leading / trailing nucleotides \\

+    & \Rfunction{tables} & Summarize read occurrence \\

+    & \Rfunction{srduplicated}, etc. & Identify duplicate reads \\

+    & \Rfunction{filterFastq} & Filter reads from one file to another\\

+    \multicolumn{3}{l}{Parallel evaluation} \\

+    & \Rfunction{srapply} & \Rfunction{lapply}-like parallel evaluation \\

+    \hline

+  \end{tabular}

+  \caption{Key functions for working with FASTQ files}

+  \label{tab:fastq}

+\end{table}

+

+\paragraph{Input} FASTQ files are large so processing involves

+iteration in `chunks' using \Rfunction{FastqStreamer}

+<<stream, eval=FALSE>>=

+strm <- FastqStreamer("a.fastq.gz")

+repeat {

+    fq <- yield(strm)

+    if (length(fq) == 0)

+        break

+    ## process chunk

+}

 @

-%% 

-Here we select only those reads that map to \texttt{chr5.fa}:

-<<readAligned-filter>>=

-chr5 <- readAligned(sp, "s_2_export.txt", filter=cfilt)

-@ 

-%% 

-Filters can be `composed' to act in unison, e.g., selecting only reads

-mapping to \texttt{chr5.fa} and on the \texttt{+} strand:

-<<readAligned-compose-filter>>=

-filt <- compose(cfilt, sfilt)

-chr5plus <- readAligned(sp, "s_2_export.txt", filter=filt)

-@ 

-%% 

-Filters can subset aligned reads at other stages in the work flow,

-using a paradigm like the following:

-<<AlignedRead-filter>>=

-chr5 <- aln[cfilt(aln)]

-@ 

-%% 

-Users can easily create their own filter by writing a function that

-accepts an object of class \Rcode{AlignedRead}, and returns a logical

-vector indicating which reads in the object pass the filter. See the

-example on the \Rcode{srFilter} help page for details, and for

-information about additional built-in filters.

-

-\subsection{Exploring \ShortRead{} objects}

-

-\Robject{aln} is an object of \Rclass{\Sexpr{class(aln)}} class. It

-contains short reads and their (calibrated) qualities:

-<<aln-sread-quality>>=

-sread(aln)

-quality(aln)

-@ 

-

-The short reads are stored as a \Rclass{DNAStringSet} class.  This

-class is defined in \Rpackage{Biostrings}. It represents DNA sequence

-data relatively efficiently.  There are a number of very useful

-methods defined for \Rclass{DNAStringSet}. Some of these methods are

-illustrated in this vignette. Other methods are described in the help

-pages and vignettes of the \Rpackage{Biostrings} and

-\Rpackage{IRanges} packages.

-

-Qualities are represented as \Sexpr{class(quality(aln))}-class

-objects. The qualities in the \Robject{aln} object returned by

-\Rfunction{readAligned} are of class \Rclass{BStringSet}. The

-\Rclass{BStringSet} class is also defined in \Rpackage{Biostrings},

-and shares many methods with those of \Rclass{DNAStringSet}.

-

-The \Robject{aln} object contains additional information about

-alignments. Some of this additional information is expected from any

-alignment, whether generated by Solexa or other software.  For

-example, \Robject{aln} contains the particular sequence within a

-target (e.g., chromosomes in a genome assembly), the position (e.g.,

-base pair coordinate), and strand to which the alignment was made, and

-the quality of the alignment. The display of \Robject{aln} suggests

-how to access this information. For instance, the strand to which

-alignments are made can be extracted (as a factor with three levels

-and possibly \Rcode{NA}; the level \Rcode{"*"} corresponds to reads

-for which strand alignment is intrinsically not meaningful, whereas

-\Rcode{NA} represents the traditional concept of information not

-available, e.g., because the read did not align at all) and tabulated

-using familiar \R{} functions.

-<<chromosomes>>=

-whichStrand <- strand(aln)

-class(whichStrand)

-levels(whichStrand)

-table(whichStrand, useNA="ifany")

-@ 

-%% 

-This shows that about 

-%% 

-\Sexpr{format(100*sum(is.na(whichStrand))/ length(whichStrand))}\%{}

-%% 

-of reads were not aligned (level \Rcode{NA}).

-

-The \Robject{aln} object contains information in addition to that

-expected of all alignments. This information is accessible using

-\Rfunction{alignData}:

-<<alignData>>=

-alignData(aln)

-@ 

-%% 

-Users familiar with the \Rclass{ExpressionSet} class in

-\Rpackage{Biobase} will recognize this as an

-\Rclass{AnnotatedDataFrame}-like object, containing a data frame with

-rows for each short read. The \Rclass{AlignedDataFrame} contains

-additional meta data about the meaning of each column. For instance,

-data extracted from the Solexa export file includes:

-<<varMetadata>>=

-varMetadata(alignData(aln))

+or drawing a random sample from the file

+<<sampler, eval=FALSE>>=

+sampler <- FastqSampler("a.fastq.gz")

+fq <- yield(sampler)

+@ 

+\noindent The default size for both streams and samples is 1M records;

+this volume of data fits easily into memory. Small FASTQ files can be

+read in to memory in their entirety using \Rfunction{readFastq}; we do

+this for our sample data

+<<readFastq>>=

+fl <- system.file(package="ShortRead", "extdata", "E-MTAB-1147", 

+                  "ERR127302_1_subset.fastq.gz")

+fq <- readFastq(fl)

+@ 

+

+The result of data input is an instance of class \Rclass{ShortReadQ}

+(Table~\ref{tab:classes}). 

+\begin{table}

+  \centering

+  \begin{tabular}{ll}

+    \hline

+    \Rclass{DNAStringSet} & (\Biocpkg{Biostrings}) Short read sequences \\

+    \Rclass{FastqQuality}, etc. & Quality encodings \\

+    \Rclass{ShortReadQ} & Reads, quality scores, and ids \\

+    \hline

+  \end{tabular}

+  \caption{Primary data types in the \Biocpkg{ShortRead} package}

+  \label{tab:classes}

+\end{table}

+This class contains reads, their quality scores, and optionally the id

+of the read.

+<<ShortReadQ>>=

+fq

+fq[1:5]

+head(sread(fq), 3)

+head(quality(fq), 3)

 @ 

-%% 

-Guides to the precise meaning of this data are on the help page for

-the \Rclass{AlignedRead} class, and in the manufacturer manuals. 

-

-Simple information about the alignments can be found by querying this