%\VignetteDepends{ShortRead}
%\VignetteIndexEntry{Sequence Analysis : Lecture}
%\VignettePackage{SeattleIntro2010}

\documentclass{beamer}
\usepackage{BioconductorSlides}
\SweaveOpts{keep.source=TRUE}

%\AtBeginSection[]{%
%  \begin{frame}<beamer>{Outline}
%    \tableofcontents
%  \end{frame}}

\title{Sequence Analysis : An Introduction}
\author{Valerie Obenchain (\href{mailto:vobencha@fhcrc.org}{vobencha@fhcrc.org}) and \\ 
Martin Morgan (\href{mailto:mtmorgan@fhcrc.org}{mtmorgan@fhcrc.org})}
\institute[FHCRC]{Fred Hutchinson Cancer Research Center}
\date{9-10 December, 2010}

\begin{document}
<<setup, echo=FALSE>>=
options(width = 40)
@

\frame{\titlepage}

\section[Outline]{}
\frame{\tableofcontents}


\section{Work flow}

\subsection{Experiment}
\begin{frame}{Experiments}
  \begin{itemize}
  \item ChIP-seq
  \item Differential expression
  \item RNA-seq (alternate splicing)
  \item Metagenomic
  \item \ldots
  \end{itemize}
\end{frame}

\subsection{Technology}
\begin{frame}{Technology}
  Platforms
  \begin{itemize}
  \item Illumina / Genome Analyzer
  \item Roche / 454
  \item Applied Biosystems / SOLiD
  \item Complete Genomics
  \end{itemize}
\end{frame}

\subsection{Pre-processing}
\begin{frame}{Pre-processing}
  Vendor and third-party
  \begin{itemize}
  \item Image processing, base calling
  \item Machine quality assessment
  \item Alignment
  \end{itemize}
  Bioconductor
  \begin{itemize}
  \item Quality assessment and representation: \Biocpkg{ShortRead},
    \Biocpkg{GenomicRanges}
  \item Read remediation, trimming, primer removal, specialized
    manipulation: \Biocpkg{IRanges}, \Biocpkg{ShortRead},
    \Biocpkg{Biostrings}
  \item Specialized alignment tasks: \Biocpkg{Biostrings},
    \Biocpkg{BSgenome}
  \end{itemize}
\end{frame}

\subsection{Analysis}
\begin{frame}
  \frametitle{\Bioconductor{} Sequence Analysis Packages}
  ChIP-seq
  \begin{itemize}
  \item \Biocpkg{BayesPeak}, \Biocpkg{chipseq}, \Biocpkg{ChIPseqR},
    \Biocpkg{ChIPsim}, \Biocpkg{PICS}
  \end{itemize}
  RNA-seq and snRNA discovery
  \begin{itemize}
  \item \Biocpkg{Genominator}, \Biocpkg{rnaSeqMap},
    \Biocpkg{segmentSeq}
  \end{itemize}
  Metagenomics
  \begin{itemize}
  \item \Biocpkg{OTUbase}
  \end{itemize}
  Methyl-seq
  \begin{itemize}
  \item \Biocpkg{MEDIPS}, \Biocpkg{methVisual}
  \end{itemize}
  Related: genotyping
  \begin{itemize}
  \item \Biocpkg{GGtools}, \Biocpkg{VanillaICE} (variant calls from
    SNP arrays)
  \end{itemize}
\end{frame}

\subsection{Annotation and Integration}
\begin{frame}{Annotation and Integration}
  Annotation
  \begin{itemize}
  \item Genome coordinate / gene (and other) relationships,
    \Biocpkg{GenomicFeatures}, \Biocpkg{ChIPpeakAnno}
  \item Resources originally developed for microarray analysis --
    \Biocpkg{AnnotationDbi}, \Rpkg{org.*.db}, \Biocpkg{KEGG.db},
    \Biocpkg{GO.db}, \Biocpkg{Category}, \Biocpkg{GOstats}
  \end{itemize}
  Integration
  \begin{itemize}
  \item Digital and microarray differential expression
  \item RNAseq and gene ontology / pathway,  \Biocpkg{goseq}
  \item HapMap, 1000 genomes, UCSC, Sequence Read Archive, GEO,
    ArrayExpress, \Biocpkg{rtracklayer}, \Biocpkg{biomaRt},
    \Biocpkg{Rsamtools}, \Biocpkg{GEOquery}, \Biocpkg{SRAdb}
  \end{itemize}
\end{frame}



\begin{frame}[t]
  \frametitle{\Bioconductor{} Sequence Packages}
\begin{figure}[ht]
 \centering
  \includegraphics[width=4in,height=3in]{packageStack}
\end{figure}
\end{frame}



\section{Data I/O : Short Reads}

\subsection{ShortRead}
\subsubsection{Input and exploration}

\begin{frame}[fragile]
  \frametitle{\Biocpkg{ShortRead} data input}

<<library, results=hide>>=
library(SeattleIntro2010)
library(ShortRead)
fl <- system.file("extdata", "SRR002051.chrI-V.bam",
  package="SeattleIntro2010") 
aln <- readAligned(fl, type = "BAM")
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{The \Rclass{AlignedRead} class}
<<aln>>=
aln
table(strand(aln), useNA="always")
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Accessing reads, base quality, and other data}
<<sread>>=
head(sread(aln), 3) 
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Alphabet by cycle}
  Expectation: nucleotide use independent of cycle
<<abc>>=
alnp <- aln[strand(aln) == "+"]
abc <- alphabetByCycle(sread(alnp))
class(abc)
abc[1:6,1:4]
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Alphabet by cycle}
  \begin{columns}
    \column{.5\textwidth} 
    \Rfunction{matplot} takes a matrix and plots each column as a set
    of points
<<abc-plot,fig=TRUE,include=FALSE,prefix=FALSE>>=
tabc <- t(abc[1:4,])
matplot(tabc, type="l",
        lty=rep(1, 4))
@ 
    \column{.5\textwidth}
    \includegraphics[width=\textwidth,height=!]{abc-plot.pdf}
  \end{columns}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Quality by cycle}
  \begin{columns}
    \column{.45\textwidth} Encoded quality scores can be decoded to
    their numerical values and represented as a matrix. Calculating
    the average of the column means creates a vector of average
    quality scores across cycle.
<<quality-plot,fig=TRUE,include=FALSE,prefix=FALSE>>=
m <- as(quality(alnp), 
        "matrix")
plot(colMeans(m), type="b")
@ 
    \column{.5\textwidth}
    \includegraphics[width=\textwidth,height=!]{quality-plot.pdf}
  \end{columns}
\end{frame}

\subsubsection{Manipulation}
\begin{frame}[fragile]
  \frametitle{Recoding and updating}
  \begin{enumerate}
  \item Access the chromosome information
  \item Extract the chromosome number from the factor level
  \item Recode the chromosome number to roman (!), create new levels,
    and update the chromosome
  \item Update the \Rclass{AlignedRead}
  \end{enumerate}
<<chr-recode>>=
chrom <- chromosome(alnp)
i <- sub("S288C_([[:digit:]]+)", "\\1", levels(chrom))
levels(chrom) <- paste("chr", as.roman(i), sep="")
alnp <- renew(alnp, chromosome=chrom)
@
\end{frame}

\subsection{Rsamtools and GenomicRanges}
\begin{frame}{\Biocpkg{Rsamtools}}
  Enhancements to \Rclass{AlignedRead} and \Biocpkg{ShortRead}
  \begin{itemize}
  \item Easy to input arbitrary subset of reads
  \item Not necessary to read sequence, quality, identifier and other
    information when not necessary
  \item Reads can be aligned with indels and gaps
  \end{itemize}
\end{frame}

\begin{frame}{\software{samtools} and \Biocpkg{Rsamtools}}
  \href{http://samtools.sourceforge.net/}{\software{samtools}}
  \begin{itemize}
  \item Data format -- text (SAM) and binary (BAM)
  \item Tools to manipulate (e.g., merge), analyze (e.g., pileup) and
    view
  \item For developers -- bindings to other languages, e.g.,
    \href{http://picard.sourceforge.net/index.shtml}{Picard}
  \end{itemize}
  \Biocpkg{Rsamtools}
  \begin{itemize}
  \item Input and represent BAM files.
  \item High-level: \Rfunction{readAligned(..., type="BAM")};
    \Rfunction{readGappedAlignments}; \Rfunction{readPileup}
  \item Flexible: \Rfunction{scanBam}
  \end{itemize}
\end{frame}

\begin{frame}{Input}
  \Rclass{ScanBamParam}
  \begin{description}
  \item[which] \Rclass{GRanges} selecting reference, genome
    coordinates, strand.
  \item[flag] select paired / mapped / mate mapped reads
  \item[what] fields to retrieve, e.g., query name, reference name,
    strand, position, width, cigar
  \end{description}
\end{frame}

\begin{frame}{GenomicRanges : Gapped alignments}
  The \Rclass{GappedAlignments} class in \Biocpkg{GenomicRanges} 
  \begin{itemize}
  \item \Rfunction{readGappedAlignments} uses \Rfunction{scanBam}
  \item Genomic coordinates, `cigar', covered intervals
  \item Cigar: run length encoding; M (match), I, D (insertion,
    deletion), N (skipped), S, H (soft, hard clip), P (padding). E.g.,
    \texttt{35M}, \texttt{18M2I15M}
  \item Accessors, subsetting, narrowing, \Rfunction{pintersect},
    \Rfunction{coverage}, \ldots
  \end{itemize}
\end{frame}


\begin{frame}[fragile]
  \frametitle{Example}
<<scanBam-2, eval=FALSE>>=
## reads on chr III overlapping 100000-110000
which <- GRange("chrIII", IRanges(100000, 110000))
param <- ScanBamParam(which=which)
bf <- scanBam("path/to/bamfile", param=param)
@ 
  \Rfunction{scanBam} returns a nested list
  \begin{itemize}
  \item One element for each row of \Rfunction{GRanges}
  \item Nested elements correspond to \Rfunarg{what}
  \end{itemize}
\end{frame} 

\section{Other Data Examples}

\begin{frame}[fragile]
  \frametitle{454 Microbiome Pre-Processing}
<<bar, eval=FALSE>>=
library(ShortRead)
dir <- "/not/public"
bar <- read454(dir)              # Input
code <- narrow(sread(bar), 1, 8) # Extract bar code
aBar <- bar[code == "AAGCGCTT"]  # Subset one bar code
noBar <-                         # Remove bar code
    narrow(aBar, 11, width(aBar)) 
pcrPrimer <- "GGACTACCVGGGTATCTAAT"
trimmed <-                       # Remove primer
    trimLRPatterns(pcrPrimer, noBar, Lfixed=FALSE)
writeFastq(trimmed,              # Output
           file.path(dir, "trimmed.fastq"))
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Digital Gene Expression}
<<DGE, eval=FALSE>>=
library(GenomicFeatures)
bamFile <- "/path/to/file.bam"
aligns <- readGappedAlignments(bamFile)
## ... txdb: transcripts from UCSC 'knownGenes'
exonRanges <- exonsBy(txdb, "tx") 
## ... housekeeping
counts <- countOverlaps(exonRanges, aligns)
## ... normalization --> 'highScores' variable
txs <- transcripts(txdb, 
           vals=list(tx_id=names(highScores)), 
           columns=c("tx_id","gene_id"))
systematicNames <- elementMetadata(txs)[["gene_id"]]
@ 
\end{frame}


\section{Resources}

\begin{frame}{Resources}
  \Bioconductor{} Web site
  \begin{itemize}
  \item \url{http://bioconductor.org}
  \item `\href{http://bioconductor.org/docs/install/}{Installation}',
    `\href{http://bioconductor.org/packages/release/Software.html}{Software}',
    and `\href{http://bioconductor.org/docs/mailList.html}{Mailing
      lists}' links.
  \end{itemize}
  Help in \R{}
  \begin{itemize}
  \item \Rfunction{help.start()} to view a help browser.
  \item \text{help(package = "Biostrings")}
  \item \texttt{?readAligned}
  \item \Rfunction{browseVignettes("GenomicRanges")}
  \end{itemize}
\end{frame}


\end{document}

