%\VignetteDepends{IRanges, GenomicRanges, Biostrings, BSgenome, GenomicFeatures, BSgenome.Scerevisiae.UCSC.sacCer2, SeattleIntro2010}
%\VignetteIndexEntry{DNA strings and ranges}
%\VignettePackage{SeattleIntro2010}

\documentclass[9pt]{beamer}

\usepackage{BioconductorSlides}
\SweaveOpts{keep.source=TRUE, eps=FALSE, width=9, height=3}

\AtBeginSection[]
{
  \begin{frame}<beamer>{Outline}
    \tableofcontents[currentsection,currentsubsection]
  \end{frame}
}

\title{Working with DNA strings and ranges}
\author{Herv\'e Pag\`es}
\institute[FHCRC]{Fred Hutchinson Cancer Research Center}
\date{9-10 December, 2010}

\begin{document}

<<setup, echo=FALSE>>=
options(width = 64)
library(GenomicRanges)
library(Biostrings)
library(SeattleIntro2010)
plotGRanges <- function(x, main = deparse(substitute(x)), sep = 0.5, ...) 
{
  height <- 1
  bins <- disjointBins(IRanges(start(x), end(x) + 1))
  plot.new()
  plot.window(c(0, 50), c(0, max(bins)*(height + sep)))
  ybottom <- bins * (sep + height) - height
  col <- ifelse(as.vector(strand(x) == "+"), "blue",
         ifelse(as.vector(strand(x) == "-"), "red", "purple"))
  rect(start(x) - 0.5, ybottom, end(x) + 0.5, ybottom + height, col = col, ...)
  title(main)
  axis(1)
}
@

\frame{\titlepage}

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}

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

\begin{frame}[fragile]
  \frametitle{Sequence infrastructure packages}
  \begin{block}{\Biocpkg{IRanges}}
    Lowest-level containers \& tools (\emph{not} Biology-specific)
  \end{block}

  \begin{block}{\Biocpkg{GenomicRanges}}
    Representation \& manipulation of genomic intervals
  \end{block}

  \begin{block}{\Biocpkg{Biostrings}}
    Representation \& manipulation of long biological sequences
    (DNA, RNA or amino-acids)
  \end{block}

  \begin{block}{\Biocpkg{BSgenome}}
    High level containers for full genomes
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Some of the most important containers}
  \begin{block}{In the \Biocpkg{IRanges} package}
    \begin{itemize}
      \item \Rclass{IRanges}
      \item \Rclass{Rle}
      \item \Rclass{Views}
      \item \Rclass{DataFrame}
    \end{itemize}
  \end{block}

  \begin{block}{In the \Biocpkg{GenomicRanges} package}
    \begin{itemize}
      \item \Rclass{GRanges}
      \item \Rclass{GRangesList}
    \end{itemize}
  \end{block}

  \begin{block}{In the \Biocpkg{Biostrings} package}
    \begin{itemize}
      \item \Rclass{DNAStringSet}
    \end{itemize}
  \end{block}

  \begin{block}{In the \Biocpkg{BSgenome} package}
    \begin{itemize}
      \item \Rclass{BSgenome}
    \end{itemize}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{\Rclass{Rle} objects}
  \begin{block}{Issue}
    \begin{itemize}
      \item Chromosomes can be hundreds of million of base
            pairs long, making them hard to manage in computer
            memory
      \item Fortunately, coverage vectors tend to follow an
            integer step function
    \end{itemize}
  \end{block}

  \begin{block}{Solution}
    \begin{itemize}
      \item Run-length encoding (RLE) is a common compression
            technique for storing long sequences with lengthy
            repeats
      \item An RLE couples values with run lengths, e.g. the
            vector 0, 0, 0, 1, 1, 2 would be represented as
            (3) 0's, (2) 1's, and (1) 2
      \item The \Biocpkg{IRanges} package uses the
            \Rclass{Rle} and \Rclass{RleList} classes to house
            coverage vectors
    \end{itemize}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{\Rclass{Views} objects}
  \begin{block}{Issue}
    \begin{itemize}
      \item Chromosomes can be hundreds of million of base
            pairs long, making subsequence selection inefficient
    \end{itemize}
  \end{block}

  \begin{block}{Solution}
    \begin{itemize}
      \item Store the original sequence using a pass-by-reference
            semantic
      \item Associate ranges with the sequence to select subsequence
      \item Example:
        \begin{itemize}
          \item 7007-letter sequence:
            \Rcode{$<<$SNIP-3000$>>$AGATTCA$<<$SNIP-4000$>>$}
          \item View range: \Rcode{[3001, 3007]}
          \item $=>$ 7-letter subsequence: \Rcode{AGATTCA}
        \end{itemize}
    \end{itemize}
  \end{block}
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Genomic ranges}

\begin{frame}[fragile]
  \frametitle{Naive representation for intervals with data}
  \begin{block}{Data characteristics}
    \begin{itemize}
      \item Genomic coordinates consist of chromosome, position,
        and potentially strand information
      \item May have additional values, such as GC content or
        alignment coverage
    \end{itemize}
  \end{block}

  \begin{block}{\Rclass{data.frame} approach}
{\small
<<naive-GRanges-1>>=
chr <- c("chr1", "chr2", "chr1")
strand <- c("+", "+", "-")
start <- c(3L, 4L, 1L)
end <- c(7L, 5L, 3L)
naive <- data.frame(chr = chr, strand = strand,
                    start = start, end = end)
@
}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{\Rclass{GRanges} construction}
  \begin{block}{\Rclass{GRanges} constructor}
    \begin{itemize}
      \item Instances are created using the \Rfunction{GRanges}
        constructor
      \item Starts and ends are wrapped in an \Rclass{IRanges}
        constructor
      \item Chromosome/contig supplied to \Rfunarg{seqnames}
        argument
      \item Underlying sequence lengths can be suppled to
        \Rfunarg{seqlengths} argument
    \end{itemize}
  \end{block}

  \begin{block}{\Rclass{GRanges} example}
{\small
<<GRanges-construction>>=
bioc <- GRanges(seqnames = chr,
                ranges = IRanges(start = start, end = end),
                strand = strand,
                seqlengths = c(chr1 = 24, chr2 = 18))
@
}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{\Rclass{GRanges} display}
  \begin{block}{\Rclass{GRanges} \Rfunction{show} method}
{\small
<<GRanges-display>>=
bioc
@
}
  \end{block}

  \begin{block}{Note}
    \begin{itemize}
      \item Optional interval data would appear to the right of
        \Rcode{|} divider
    \end{itemize}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Some common operations on \Rclass{GRanges} objects}
  \begin{block}{}
    \begin{itemize}
      \item \Rfunction{reduce}
      \item \Rfunction{gaps}
      \item \Rfunction{coverage}
      \item \Rfunction{intersect} \& \Rfunction{pintersect}
      \item \Rfunction{findOverlaps} \& \Rfunction{countOverlaps}
    \end{itemize}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Creating a new \Rclass{GRanges} object}
  \begin{block}{We will use the following \Rclass{GRanges} objects to
                illustrate some of those operations}
{\small
<<ranges-construct>>=
ir <- IRanges(c(1, 8, 14, 15, 19, 34, 40), 
              width=c(12, 6,  6, 15, 6, 2, 7))
strand <- rep(c("+", "-"), c(4, 3))
gr <- GRanges(seqnames = "chr1", ranges = ir, strand = strand,
              seqlengths = c(chr1 = 50))
gr
@
}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{}
<<ranges-plot, fig=TRUE, echo=FALSE>>=
plotGRanges(gr)
@
\newline
<<ranges-reduce-plot, echo=FALSE, fig=TRUE>>=
plotGRanges(reduce(gr))
@
\newline
<<ranges-gaps-plot, echo=FALSE, fig=TRUE>>=
plotGRanges(gaps(gr))
@
{\footnotesize
\newline
  blue = positive strand, red = negative strand
}
\end{frame}

\begin{frame}[fragile]
  \frametitle{\Rclass{GRanges} subsetting}
  \begin{block}{}
{\small
<<ranges-subset>>=
gr[strand(gr) == "-"]
@
}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Finding overlapping ranges}
  \begin{block}{}
    Given 2 objects containing ranges, \Rfunction{findOverlap} (or its
    variant \Rfunction{countOverlaps}) finds (or counts) the overlaps
    between the 2 objects
  \end{block}
{\small
<<ranges-overlap>>=
ol <- findOverlaps(gr, reduce(gr))
as.matrix(ol)
countOverlaps(reduce(gr), gr)
@
}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Coverage of the underlying sequence}
  \begin{block}{}
    \begin{itemize}
      \item \Rfunction{coverage} counts the number of ranges over
      each position in the underlying sequence
      \item Subset by strand to get stranded coverage
    \end{itemize}
  \end{block}
  \begin{block}{}
{\small
<<ranges-coverage>>=
cvg <- coverage(gr)
@
}
<<ranges-coverage-plot, height=4, echo=FALSE, fig=TRUE>>=
cvg <- coverage(gr)
plotGRanges(gr, main="coverage(gr)")
cvg <- as.integer(cvg[["chr1"]])
mat <- cbind(seq_along(cvg)-0.5, cvg)
d <- diff(cvg) != 0
mat <- rbind(cbind(mat[d,1]+1, mat[d,2]), mat)
mat <- mat[order(mat[,1]),]
lines(mat, col="black", lwd=4)
axis(2)
@
  \end{block}
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Long biological sequences}

\begin{frame}[fragile]
  \frametitle{DNA sequences}
  \begin{block}{}
{\scriptsize
<<DNAString>>=
library(Biostrings)
data(yeastSEQCHR1)
c(class(yeastSEQCHR1), nchar(yeastSEQCHR1))
yeast1 <- DNAString(yeastSEQCHR1)
yeast1
IUPAC_CODE_MAP
data(srPhiX174)
length(srPhiX174)
head(srPhiX174, 3)
@
}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Basic string utilities}
  \begin{block}{Subsequence selection}
    \Rfunction{subseq}, \Rfunction{Views}
  \end{block}

  \begin{block}{Letter frequencies}
    \Rfunction{alphabetFrequency},
    \Rfunction{dinucleotideFrequency},
    \Rfunction{trinucleotideFrequency},
    \Rfunction{oligonucleotideFrequency},
    \Rfunction{letterFrequencyInSlidingView},
    \Rfunction{uniqueLetters}
  \end{block}

  \begin{block}{Letter consensus}
    \Rfunction{consensusMatrix}, \Rfunction{consensusString}
  \end{block}

  \begin{block}{Letter transformation}
    \Rfunction{reverse}, \Rfunction{complement},
    \Rfunction{reverseComplement}, \Rfunction{translate},
    \Rfunction{chartr}
  \end{block}

  \begin{block}{I/O}
    \Rfunction{read.DNAStringSet}, \Rfunction{read.RNAStringSet},
    \Rfunction{read.AAStringSet}, \Rfunction{read.BStringSet},
    \Rfunction{write.XStringSet}, \Rfunction{save.XStringSet}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{String matching/alignment utilities}
  \begin{block}{\Rfunction{matchPDict}}
    \Rfunction{matchPDict}, \Rfunction{countPDict},
    \Rfunction{whichPDict}, \Rfunction{vmatchPDict},
    \Rfunction{vcountPDict}, \Rfunction{vwhichPDict}
  \end{block}

  \begin{block}{\Rfunction{vmatchPattern}}
    \Rfunction{matchPattern}, \Rfunction{countPattern},
    \Rfunction{vmatchPattern}, \Rfunction{vcountPattern}
  \end{block}

  \begin{block}{\Rfunction{pairwiseAlignment}}
    \Rfunction{pairwiseAlignment}, \Rfunction{stringDist}
  \end{block}

  \begin{block}{\Rfunction{matchPWM}}
    \Rfunction{matchPWM}, \Rfunction{countPWM}
  \end{block}

  \begin{block}{OTHER}
    \Rfunction{matchLRPatterns}, \Rfunction{trimLRPatterns},
    \Rfunction{matchProbePair}, \Rfunction{findPalindromes},
    \Rfunction{findComplementedPalindromes}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Letter frequencies}
  \begin{block}{Single-letter frequencies}
{\footnotesize
<<alphabetFrequency>>=
alphabetFrequency(yeast1, baseOnly=TRUE)
@
}
  \end{block}

  \begin{block}{Multi-letter frequencies}
{\footnotesize
<<dinucleotideFrequency>>=
dinucleotideFrequency(yeast1)
head(trinucleotideFrequency(yeast1), 12)
@
}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Basic transformations}
  \begin{block}{}
{\small
\begin{verbatim}
> x
  21-letter "DNAString" instance
seq: TCAACGTTGAATAGCGTACCG
> reverseComplement(x)
  21-letter "DNAString" instance
seq: CGGTACGCTATTCAACGTTGA
> translate(x)
  7-letter "AAString" instance
seq: STLNSVP
\end{verbatim}
}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{}
  \begin{block}{Consensus matrix}
{\footnotesize
<<consensusMatrix>>=
snippet <- subseq(head(sort(srPhiX174), 5), 1, 10)
consensusMatrix(snippet, baseOnly=TRUE)
@
}
  \end{block}

  \begin{block}{Consensus string}
{\footnotesize
<<consensusString>>=
consensusString(snippet)
consensusString(snippet, ambiguityMap = "N", threshold = 0.5)
@
}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Pattern matching}
  \begin{block}{Count the matches only}
{\footnotesize
<<countPDict>>=
data(phiX174Phage)
genome <- phiX174Phage[["NEB03"]]
negPhiX174 <- reverseComplement(srPhiX174)
posCounts <- countPDict(PDict(srPhiX174), genome)
negCounts <- countPDict(PDict(negPhiX174), genome)
table(posCounts, negCounts)
@
}
  \end{block}

  \begin{block}{Find the locations of the matches}
{\footnotesize
<<matchPDict>>=
matchPDict(PDict(srPhiX174[posCounts > 0]), genome)
@
}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Pairwise alignments}
  \begin{block}{Alignment scores}
{\footnotesize
<<pairwiseAlignments-1>>=
data(phiX174Phage)
posScore <- pairwiseAlignment(srPhiX174, genome,
              type = "global-local", scoreOnly = TRUE)
negScore <- pairwiseAlignment(negPhiX174, genome,
              type = "global-local", scoreOnly = TRUE)
cutoff <- max(pmin.int(posScore, negScore))
@
}
  \end{block}

  \begin{block}{Alignments}
{\footnotesize
<<matchPDict>>=
pairwiseAlignment(srPhiX174[posScore > cutoff], genome,
                  type = "global-local")
@
}
  \end{block}
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Full genomes}

\begin{frame}[fragile]
  \frametitle{\textit{BSgenome data packages}}
  \begin{block}{}
    \begin{itemize}
      \item One genome per package
      \item Full genome sequences stored in \Biocpkg{Biostrings} containers
      \item 14 organisms / 24 packages in the current release (BioC 2.7)
      \item Most (but not all) packages contain sequences with builtin masks
      \item Naming convention: BSgenome.\textit{Organism}.\textit{Provider}.\textit{BuildVersion}
    \end{itemize}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Available genomes}
  \begin{block}{}
    Use the \Rfunction{available.genomes} function (from the
    \Biocpkg{BSgenome} software package) to get the list:
{\scriptsize
<<available.genomes>>=
library(BSgenome)
available.genomes()
@
}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Making your own \textit{BSgenome data package}}
  \begin{block}{}
    \begin{itemize}
      \item It's easy to make your own package if your organism is not
            supported
      \item Documented in the BSgenomeForge vignette from the
            \Biocpkg{BSgenome} software package
    \end{itemize}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Working with a genome}
  \begin{block}{}
{\scriptsize
\begin{verbatim}
> library(BSgenome.Celegans.UCSC.ce2)
> Celegans
Worm genome
| 
| organism: Caenorhabditis elegans (Worm)
| provider: UCSC
| provider version: ce2
| release date: Mar. 2004
| release name: WormBase v. WS120
| 
| single sequences (see '?seqnames'):
|   chrI    chrII   chrIII  chrIV   chrV    chrX    chrM    
| 
| multiple sequences (see '?mseqnames'):
|   upstream1000  upstream2000  upstream5000  
| 
| (use the '$' or '[[' operator to access a given sequence)
> Celegans$chrI
  15080483-letter "DNAString" instance
seq: GCCTAAGCCTAAGCCTAAGCCTAAGCCTAAG...CTTAGGCTTAGGCTTAGGTTTAGGCTTAGGC
\end{verbatim}
}
  Note that a given chromosome sequence is loaded from disk
  to memory at the time the user makes reference to it (like
  in the above example), and is unloaded (i.e. the memory is
  freed) when there are no more references to it.
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Working with a genome (continued)}
  \begin{block}{}
{\small
\begin{verbatim}
> seqlengths(Celegans)
    chrI    chrII   chrIII    chrIV     chrV     chrX     chrM 
15080483 15279308 13783313 17493791 20922231 17718849    13794 
\end{verbatim}
}
  \end{block}
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Resources}

\begin{frame}[fragile]
  \frametitle{Resources}
  \begin{block}{\Bioconductor{} web site}
    \begin{itemize}
      \item \url{http://bioconductor.org}
      \item `Help' link
    \end{itemize}
  \end{block}

  \begin{block}{Help in \R{}}
    \begin{itemize}
      \item \Rcode{help.start()} to view a help browser.
      \item \Rcode{help(package = "Biostrings")}
      \item \Rcode{?findOverlaps}
      \item \Rcode{browseVignettes("GenomicRanges")}
    \end{itemize}
  \end{block}
\end{frame}


\end{document}
