%\VignetteDepends{ShortRead, GenomicFeatures, edgeR}
%\VignetteIndexEntry{High-Throughput Sequence Annotation Lab}
%\VignettePackage{IWB2011}

\documentclass{article}

\usepackage{BioconductorVignette}
\SweaveOpts{keep.source=TRUE}

\title{High Throughput Sequence Manipulation and Annotation}
\author{Martin Morgan \\ Fred Hutchinson Cancer Research Center}
\date{19-21 January, 2011, edited 31 December, 2010}

\begin{document}

<<setup, echo=FALSE>>=
options(width = 70)
jnk <- capture.output(library(locfit))
library(IWB2011)
@ 

\maketitle

\tableofcontents

\section{Introduction}

This lab uses data from Nagalakshmi et al.~\cite{pmid18451266}. This
is the publication that defined the term `RNA-seq'; it is primarily
methodological. The authors investigated two different approaches to
generating DNA from poly(A) RNA, via `random hexamer' (RH), and
`oligo(dT)' (dT). There were both biological and technical
replicates. Sequencing used the Illumina GAI platform, so there are
relative few reads ($<$ 5 million per lane) of short (33bp) reads.

\section{Data input and exploration}

\subsection{Sequences}

Attach the \Rpackage{ShortRead} and \Rpackage{IWB2011} packages, and
locate the file path to the sample data set of the \Rpackage{IWB2011}
package.
<<sequence-fastqFile>>=
library(ShortRead)
library(IWB2011)
fastqFile <- system.file("extdata", "SRR002051.reads1-50k.fastq",
                         package="IWB2011")
@ 
Hint: if \Robject{fastqFile} is \Rcode{""}, then the
\Rfunction{system.file} argument is incorrect.

Input the fastq sequences, and explore the resulting object.
<<sequence-readFastq>>=
fq <- readFastq(fastqFile)
fq
head(sread(fq))
head(quality(fq))
@ 

As a simple illustration of tasks that can be performed, use the
\Rfunction{alphabetByCycle}, \Rfunction{t} (matrix transposition), and
\Rfunction{matplot} (`matrix plot') functions to visualize how
nucleotide use changes with cycle.
<<sequence-matplot, fig=TRUE>>=
abc <- alphabetByCycle(sread(fq))
matplot(t(abc[1:4,]), type="l", xlab="Cycle", ylab="Count", lwd=2)
@ 

As a second illustration, coerce quality scores to a matrix using
\Rfunction{as}. Then use \Rfunction{apply} and \Rfunction{mean} to
summarize the mean base quality per read. Create a subset of `high
quality' reads with mean quality greater than 20, using the
\Rfunction{[} operator and a (vectorized) logical test.
<<sequence-quality>>=
m <- as(quality(fq), "matrix")
baseQualPerRead <- apply(m, 1, mean)
fq[baseQualPerRead > 20]
@ 

\subsection{Gapped alignments}

Attach the \Rpackage{GenomicRanges} package, determine the path to the
sample BAM file, and read in the alignments using the
\Rfunction{readGappedAlignments} function. Use \Rfunction{head} to
view the first half dozen records.
<<gapped-file>>=
library(GenomicRanges)
bamFile <- system.file("extdata", "SRR002051.chrI-V.bam",
                       package="IWB2011")
galn <- readGappedAlignments(bamFile)
head(galn)
@ 
Be sure to verify that the file path has been specified correctly!

Use the \Rfunction{cigar} accessor to obtain the CIGAR code of each
read. Use the \Rfunction{table}, \Rfunction{sort} and \Rfunction{tail}
functions to summarize the most commonly occurring CIGARs.
<<gapped-input>>=
tail(sort(table(cigar(galn))))
@ 

\section{Annotation: transcripts}

Sequence-based annotations are generally more complicated than the
simple gene-level annotations from microarray data. The following
script downloads and saves a \href{http://genome.ucsc.edu/}{UCSC
  genome browser} track describing annotated genes in yeast
<<annotation>>=
readScript("create-txdb.R")
@ 
%% 
Load this file into your current session, and find out more
information about it with the following commands:
<<annotation-load>>=
library(GenomicFeatures)
txdbFile <- system.file("extdata", "sacCer2_sgdGene.sqlite",
                        package="IWB2011")
txdb <- loadFeatures(txdbFile)
txdb
@ 

Query the \Robject{txdb} object for transcripts using the
\Rfunction{transcripts} function and \Rcode{columns="gene\_id"}
argument. Summarize the number of transcripts associated with each
gene symbol by: (a) extracting the \Rcode{gene\_id} column using
\Rfunction{values} and column selectin; and (b) using
\Rfunction{table} twice, first to count the number of times each gene
identifier occurs, and then the number of times any gene identifier
occurs once, twice, \ldots.
<<annotation-transcripts>>=
tscript <- transcripts(txdb, columns="gene_id")
head(tscript)
geneIds <- as(values(tscript)$gene_id, "character")
table(table(geneIds))
## following line changes how 'gene_id' is stored, internally
values(tscript)$gene_id <- geneIds
@ 

For simplicity in subsequent analyses, select the subset of genes with
exactly one transcript. Do this using \Rfunction{table},
\Rfunction(names), and logical subsetting to identify candidate genes,
and \Rfunction{\%in\%} to select the corresponding rows of
\Robject{tscript}. Confirm that the resulting subset is consistent
with the tabulation performed in the previous exercise.
<<annotation-unique-tscript>>=
x <- table(geneIds)
ugeneIds <- names(x)[x==1]
idx <- values(tscript)$gene_id %in% ugeneIds
tscript1 <- tscript[idx]
length(tscript1)
@ 
\section{Counting and measuring}

The goal is to count the number of times a read `hits' a
transcript. We need to make decisions about what constitutes a `hit',
and what to do with reads that hit more than one transcript. We define
a hit as any overlap between the aligned read and a transcript, and we
discard reads that align to more than one transcript. The RNAseq
protocol used in this experiment does not distinguish between strands,
so in preparation for further analysis we indicate that the strand of
the transcript is unimportant
<<count-strand-blat>>=
strand(tscript1) <- "*"
@ 

Use the \Rpackage{IRanges} package and \Rfunction{countOverlaps}
function to count the number of times each read aligns to a
transcript. Subset the reads to include only those that align to
exactly one transcript. Use \Rfunction{countOverlaps} again, but with
arguments reversed, to count the number of times each transcript
overlaps a read.
<<count-union>>=
library(IRanges)
hits <- countOverlaps(galn, tscript1)
table(hits)
galn1 <- galn[hits==1]
counts <- countOverlaps(tscript1, galn1)
@ 

Explore the \Robject{counts} object to ensure that the results are
reasonable. For instance, non-zero counts should only involve
chromosomes I-V (because the subset of reads we are using are from
these chromosomes). Visualize the non-zero counts as a histogram; use
the \Rfunction{asinh} (inverse hyperbolic sine) function for a
log-like transformation of the data.
<<count-sanity-check, fig=TRUE>>=
table(seqnames(tscript1[counts!=0]))
hist(asinh(counts[counts!=0]))
@ 

\section{Differential representation}

This section is optional; greater detail on differential expression
will be provided tomorrow.

Code presented in the appendix was used to create an object
summarizing the entire data set.  Load the data set into \R{}. The
object is a kind of \Robject{data.frame}, designed for large data. A
feature is that the columns can be annotated with additional
information. Access this information with the
\Rfunction{elementMetadata} function. Coerce \Robject{tscriptCounts}
to a standard \Robject{data.frame} and summarize the
\Rfunction{asinh}-transformed counts across samples using
\Rfunction{splom} function from the \Rpackage{lattice} package. Use
the \Rfunarg{pch} argument to set the plot character to ``.''. Can you
identify, `by eye', the technical and biological replicates within
each protocol?
<<differential-representation-explore, fig=TRUE>>=
data("tscriptCounts")
head(tscriptCounts)
elementMetadata(tscriptCounts)
df <- as(tscriptCounts, "data.frame")
print(splom(asinh(df), pch="."))
@ 

Create a subset of the data that includes only biological
replicates. Do this by creating an index variable based on the
`Replicate' column of the \Rfunction{elementMetadata} of
\Robject{tscriptCounts}, and using logical subsetting. Remove any gene
for which no transcripts were found by using \Rfunction{rowSums} to
sum over the number of reads per row, and creating a logical index
variable that is \Rcode{TRUE} only when the row contains some
reads. Use this to subset \Robject{tscriptCounts1}.
<<differential-representation-select>>=
idx <- elementMetadata(tscriptCounts)$Replicate
tscriptCounts1 <- tscriptCounts[,idx != "Technical"]
ridx <- rowSums(as(tscriptCounts1, "data.frame")) != 0
tscriptCounts1 <- tscriptCounts1[ridx,]
@ 

We use the \Rpackage{edgeR} package to assess differential
representation. This package requires counts (created by coercing
\Robject{tscriptCounts1} to a \Robject{data.frame}) and a variable to
indicate groups to compare (from the `Protocol' column of the
\Rfunction{elementMetadata}), presented to the \Rfunction{DGEList}
function. \Rpackage{edgeR} implements methods for normalizing library
representations (with the \Rfunction{calcNormFactors} function),
estimating dispersion parameters of the negative binomial model
(\Rfunction{estimateCommonDisp}), performing exact tests of
statistical significance (\Rfunction{exactTest}), and presenting a
`top table' of differentially represented tags
(\Rfunction{topTags}). There are important statistical decisions at
each step of the analysis; evaluate the following command as a `fast
track' through this work flow.
<<differential-representation>>=
library(edgeR)
df <- as(tscriptCounts1, "data.frame")
grps <- elementMetadata(tscriptCounts1)$Protocol
dge <- DGEList(df, group=grps)
dge <- calcNormFactors(dge)
dge <- estimateCommonDisp(dge)
etest <- exactTest(dge)
ttags <- topTags(etest, n=10, adjust.method="BH")
@ 

Explore these results. For instance, use \Rfunction{hist} to visualize
the distribution of (unadjusted) $P$ values available in the
\Robject{etest} object.
<<differential-representation-hist, fig=TRUE>>=
hist(etest$table$p.value, 21)
@ 

\section{Annotation: biological function}

We can find out about the biological features identified as
differentially represented. For instance, use \Rfunction{rownames} to
extract the most differentially represented transcripts in
\Robject{ttags}. Use these to query \Robject{tscript} for the genomic
locations of these transcripts.
<<annotation-biology>>=
geneIds <- rownames(ttags$table)
idx <- values(tscript)$gene_id %in% geneIds
(topScripts <- tscript[idx])
@ 

The \Rpackage{BSgenome.*} packages provide DNA sequences of model
organisms. Attach the \Rpackage{BSgenome.Scerevisiae.UCSC.sacCer2}
package and display a summary of its content using the organism part
of the package name. Use \Rfunction{getSeq} to retrieve the sequences
corresponding to these regions, using the \Rfunarg{as.character=FALSE}
argument to return the results in a \Rclass{DNAStringSet} object. Use
the \Rfunction{alphabetFrequency} function with the
\Rfunarg{baseOnly=TRUE} argument from \Rpackage{Biostrings} to
summary nucleotide use. What is the GC content of the top transcripts?
<<annotation-BSgenome>>=
library(BSgenome.Scerevisiae.UCSC.sacCer2)
Scerevisiae
(topSeqs <- getSeq(Scerevisiae, topScripts, as.character=FALSE))
(abc <- alphabetFrequency(topSeqs, baseOnly=TRUE))
rowSums(abc[,c("G", "C")]) / rowSums(abc)
@ 

The \Rpackage{org.*} packages represent a snapshot of organism-centric
annotation information collated before each \Bioconductor{}
release. Attach the \emph{Saccharomyces cerivisiea} annotation package
created from SGD annotations. Use \Rfunction{mget} to look up the
GENENAME corresponding to gene ids of differentially represented
transcripts; only a couple of these transcripts have names, use
\Rfunction{Filter} and \Rfunction{Negate} to select those with names
<<annotation-org>>=
library(org.Sc.sgd.db)
(namedIds <- Filter(Negate(is.na), mget(geneIds, org.Sc.sgdGENENAME)))
@ 

There are many other annotation resources available in
\Bioconductor{}. Explore, on your own, the \Rpackage{biomaRt} and
\Rpackage{rtracklayer} packages.

\bibliographystyle{plain}
\bibliography{IWB2011}

\appendix

\section{Appendix: counting overlaps}

This script counts reads overlapping transcripts across six samples.
<<appendix-differential-representation-create>>=
readScript("create-tscriptCounts.R")
@ 

\end{document}

%% Local Variables:
%% compile-command: "R CMD make -f devmakefile Annotation-Lab.pdf"
%% End:
