%\VignetteDepends{ShortRead,rtracklayer,GenomicFeatures,Rsamtools,org.Sc.sgd.db}
%\VignetteIndexEntry{Using the GenomicFeatures package}
%\VignettePackage{SeattleIntro2010}

\documentclass{beamer}

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

\title{Using the \Biocpkg{GenomicFeatures} package}
\author{Marc Carlson }
\institute[FHCRC]{Fred Hutchinson Cancer Research Center}
\date{December 10th 2010}

\begin{document}

<<setup, echo=FALSE>>=
options(width = 64)
@

\frame{\titlepage}

\section{Introduction}

\begin{frame}[t]
  \frametitle{\Bioconductor{} Annotation Packages: a bigger picture}
\begin{figure}[ht]
 \centering
  \includegraphics[width=4in,height=3in]{databaseTypes}
\end{figure}
\end{frame}

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

\begin{frame}[t]
  \frametitle{\Bioconductor{} Sequence Annotations}
  \begin{block}{\Biocpkg{The GenomicFeatures package}}
    \begin{itemize}
      \item Transcript information stored in SQLite databases
      \item These databases attempt to represent the relational 
        nature of splicing data correctly
      \item Transcript information exposed to user via the
        \Biocpkg{GenomicRanges} infrastructure
      \item Parsers are available to construct these databases from biomaRt 
        or the UCSC genome browser.  Also there is a generic parser for 
        custom jobs.
    \end{itemize}
  \end{block}
\end{frame}


\section{GenomicFeatures}

\begin{frame}[t,fragile]
  \frametitle{\Biocpkg{GenomicFeatures} transcript sources}
  \begin{block}{Constructors}
    \Rfunction{makeTranscriptDbFromBiomart},
    \Rfunction{makeTranscriptDbFromUCSC}
  \end{block}

{\footnotesize
<<loadGenomicFeatures>>=
library(GenomicFeatures)
nrow(supportedUCSCtables())
head(supportedUCSCtables(), 10)
@
}
\end{frame}

\begin{frame}[t,fragile]
  \frametitle{\Rclass{TranscriptDb} basics}
  \begin{block}{Making a \Rclass{TranscriptDb} object}
{\small
<<makeUCSC, eval=FALSE>>=
mm9KG <-
  makeTranscriptDbFromUCSC(genome = "mm9",
                           tablename = "knownGene")
@
}
  \end{block}
  \begin{block}{Saving and Loading}
{\small
<<saveFeatures, eval=FALSE>>=
saveFeatures(mm9KG, file="mm9_knownGene.sqlite")
@
<<loadFeatures>>=
txdb <-
  loadFeatures(system.file("extdata", "mm9_knownGene.sqlite",
                           package = "SeattleIntro2010"))
@
}
  \end{block}
\end{frame}

\begin{frame}[t,fragile]
  \frametitle{\Rclass{TranscriptDb} class}
{\small
<<TranscriptDb>>=
txdb
@
}
\end{frame}

\begin{frame}{\Rclass{TranscriptDb} schema}
  \begin{figure}
    \includegraphics[width=4.25in]{transcriptschema.png}
  \end{figure}
\end{frame}

\begin{frame}[t,fragile]
  \frametitle{Ungrouped transcript-related information}
  \begin{block}{Extractors}
    \Rfunction{transcripts}, \Rfunction{exons}, \Rfunction{cds}
  \end{block}

{\footnotesize
<<transcripts>>=
tx <- transcripts(txdb)
length(tx)
head(tx, 5)
@
}
\end{frame}



\begin{frame}[t,fragile]
  \frametitle{Grouped transcript-related information}
  \begin{block}{Extractors}
    \Rfunction{transcriptsBy}, \Rfunction{exonsBy},
    \Rfunction{cdsBy}, \Rfunction{intronsByTranscript},
    \Rfunction{fiveUTRsByTranscript},
    \Rfunction{threeUTRsByTranscript}
  \end{block}

{\footnotesize
<<exonsBy>>=
txExons <- exonsBy(txdb)
txIntrons <- intronsByTranscript(txdb)
txExons[6]
@
}
\end{frame}


\begin{frame}[t,fragile]
  \frametitle{Standard GRanges accessors can be used}
  \begin{block}{Accessors}
    \Rmethod{names}, \Rmethod{length},
    \Rmethod{elementMetaData}, \Rmethod{width},
    \Rmethod{ranges},\Rmethod{start},\Rmethod{end}
  \end{block}
  \begin{block}{Examples:}
{\footnotesize
<<accessors>>=
head(start(tx))
head(ranges(txExons), n=1)
head(elementMetadata(tx), n=2)
@
}
  \end{block}
\end{frame}




%EXERCISE 1
\begin{frame}[fragile]
  \frametitle{GenomicFeatures exercise 1}
  Load the transciptDB object above and then gather the annotations
for transcripts as grouped by their genes.
\end{frame}

\begin{frame}[fragile]
\end{frame}


\begin{frame}[fragile]
  \frametitle{GenomicFeatures exercise 1 solution}
{\small
<<exer1, results=hide>>=
library(GenomicFeatures)
txdb <-
  loadFeatures(system.file("extdata", "mm9_knownGene.sqlite",
                           package = "SeattleIntro2010"))

myTranscripts <- transcriptsBy(txdb, by="gene")
@
}
\end{frame}


%% TODO: Add a slide that demonstrates %in% and match in a GRanges
%% context? - low priority and likely that Herve will cover this.

\begin{frame}[t,fragile]
  \frametitle{Overlapping with transcripts}
  \begin{block}{Methods}
    \Rfunction{findOverlaps}, \Rfunction{countOverlaps},
    \Rfunction{match}, \Rfunction{\%in\%},
    \Rfunction{subsetByOverlaps}
  \end{block}
\vskip12pt
  \begin{block}{Usage and help:}
{\footnotesize
<<findOverlaps-help, eval=FALSE>>=
findOverlaps(query, subject, maxgap = 0L, minoverlap = 1L,
             type = c("any", "start", "end"),
             select = c("all", "first"))
@
}
  \end{block}
  \begin{block}{Example}
{\footnotesize
<<countOverlaps>>=
grngs <- GRanges("chr1", gaps(ranges(txIntrons[[7]])), "-")
countOverlaps(grngs, tx)
@
}
  \end{block}
\end{frame}



\section{Using Annotations and Genomic Ranges together}

\begin{frame}[fragile]
  \frametitle{Using Annotations and Genomic Ranges}
  Let's begin by learning what we need to extract data for chromosomes 3 and 5
  \begin{block}{Example}
{\footnotesize
<<subselect, results=hide>>=
library(GenomicFeatures)
txdb <- makeTranscriptDbFromUCSC(genome="sacCer2", tablename="sgdGene")
saveFeatures(txdb, "sacCer2_sgdGene.sqlite")
txdb <- loadFeatures(system.file("extdata", "sacCer2_sgdGene.sqlite",
                 package = "SeattleIntro2010"))
t <- transcripts(txdb)
seqlengths(t)
## length of chr3 and 5 is 316617 and 576869
chr3Len <- seqlengths(t)["chrIII"]
chr5Len <- seqlengths(t)["chrV"]
## Define a GRanges object that describes the sequences we want to read in:
gr <- GRanges(seqnames = Rle(c("chrIII","chrV")),
              ranges = IRanges(1,width=c(chr3Len,chr5Len)),
              seqlengths = c(chr3Len, chr5Len))
## And of course we also need a bam file
fl <- system.file("extdata", "SRR002051.chrI-V.bam",
                 package = "SeattleIntro2010")
@
}
  \end{block}  
\end{frame}


% \begin{frame}[fragile]
%   \frametitle{Using Annotations and Genomic Ranges}
%   Setting up to use scanBam
%   \begin{block}{Example}
% {\footnotesize
% <<scanBam, results=hide>>=
% ## now to extract our sequences
% library(ShortRead)
% sbp <- ScanBamParam(which=gr)
% ## look at the headers will also show you the chromosome lengths
% scanBamHeader(fl) ## It's a good idea to verify this.
% ## sort the Bam file
% sort <- sortBam(fl, "SRR002051.chrI-V.sorted")
% ## make an index file
% index <- indexBam(sort, "SRR002051.chrI-V.sorted")
% ## run scanBam
% bam <- scanBam(fl, index=names(index), param=sbp)
% @
% }
%   \end{block}  
% \end{frame}


% \begin{frame}[fragile]
%   \frametitle{Using Annotations and Genomic Ranges}
%   Putting together the GRanges object by hand is complicated.
%   \begin{block}{Example}
% {\footnotesize
% <<scanBam2, results=hide>>=
% ## lets make bam into a GRanges object
% grbam <-
%   GRanges(seqnames =
%           Rle(c("chrIII","chrV"), 
%               c(length(bam[[1]][["rname"]]), length(bam[[2]][["rname"]]))),
%           ranges =
%           IRanges(start = c(bam[[1]][["pos"]], bam[[2]][["pos"]]), 
%                   width = c(bam[[1]][["qwidth"]], bam[[2]][["qwidth"]]) ),
%           seqlengths = c(chr3Len, chr5Len),
%           strand = Rle(strand(c(as.character(bam[[1]][["strand"]]), 
%             as.character(bam[[2]][["strand"]])))),
%           qual = c(bam[[1]][["qual"]], bam[[2]][["qual"]]),
%           seq = c(bam[[1]][["seq"]], bam[[2]][["seq"]]))
% @
% }
%   \end{block}
% \end{frame}


\begin{frame}[fragile]
  \frametitle{Using Annotations and Genomic Ranges}
  \begin{itemize}
    \item \Rfunction{readGappedAlignments} calls 
      \Rfunction{readBamGappedAlignments}
    \item This is a good demo of the "..." argument
    \item This also shows some of the complexity of the R help system
  \end{itemize}
  \begin{block}{Example}
{\footnotesize
<<readGappedAlignments1, results=hide>>=
## A convenient way:
library(ShortRead)
grGappedBam <- readGappedAlignments(fl)
grbam <- grg(grGappedBam)
grbam <- grbam[seqnames(grbam) %in% c("chrIII","chrV")]
seqlengths(grbam) <- c(chr3Len, chr5Len)

## Look at the help here
# ?readGappedAlignments
## and here
# ?readBamGappedAlignments
@
}
  \end{block}  
\end{frame}


\begin{frame}[fragile]
  \frametitle{Using Annotations and Genomic Ranges}
  \begin{itemize}
    \item \Rfunction{readGappedAlignments} calls 
      \Rfunction{readBamGappedAlignments}
    \item This is a good demo of the "..." argument
    \item This also shows some of the complexity of the R help system
  \end{itemize}
  \begin{block}{Example}
{\footnotesize
<<readGappedAlignments2, results=hide>>=
## The most convenient way
grGappedBam <- readGappedAlignments(fl, index=fl, 
                                    which=gr)
grbam <- grg(grGappedBam)
seqlengths(grbam) <- c(chr3Len, chr5Len)
## I saved this for you
save(grbam, file= "bamReadsChr3andChr5.rda")
@
}
  \end{block}  
\end{frame}



%% Demonstration of how we can get annotations for our data - find for
%% Ura3(YEL021W) and a neighboring gene TIM9 (YEL020W-A).  This is what
%% they will do in the exercise (but in the reverse direction.  I have to
%% show this also in reverse Starting with a range or two (that I know in
%% advance will be interesting and then using findOverlaps to match it up
%% with info.)

\begin{frame}[fragile]
  \frametitle{Practical examples: Ura3 and TIM9}
  Load in the Data and Annotations  
  \begin{block}{Example}
{\footnotesize
<<GetAnnotations, results=hide>>=
load(system.file("data", "bamReadsChr3andChr5.rda",
                 package = "SeattleIntro2010"))
## 1st get the annotations for all exons grouped by gene
library(GenomicFeatures)
txdb <- loadFeatures(system.file("extdata", "sacCer2_sgdGene.sqlite",
                 package = "SeattleIntro2010"))

txs <- transcriptsBy(txdb, by="gene")
@
}
  \end{block}
\end{frame}



%EXERCISE 2
\begin{frame}[fragile]
  \frametitle{GenomicFeatures exercise 2}
  Why did we use \Rfunction{transcriptsBy} instead of \Rfunction{transcripts}.
\end{frame}

\begin{frame}[fragile]
\end{frame}

\begin{frame}[fragile]
  \frametitle{GenomicFeatures exercise 2 solution}
  Because we want to consider the ranges associated with each gene,
  and the \Rfunction{transcripts} function would not distinguish
  between two actual transcripts that were associated with a single
  gene.
\end{frame}





\begin{frame}[fragile]
  \frametitle{Practical examples: Ura3 and TIM9}
  Subset the annotations down to just two genes so that we can overlap with 
  just those ranges.
  \begin{block}{Example}
{\footnotesize
<<GetAnnotations2, results=hide>>=
## Check if the gene IDs "YEL021W" & "YEL020W-A" are present
c("YEL021W","YEL020W-A") %in% names(txs)

## Lets extract only two genes
ind = names(txs) %in% c("YEL021W","YEL020W-A")
txpair = txs[ind]

## drop strand info. (so we will match both + and - strands)
sgrbam <- grbam
strand(sgrbam) <- rep("*",length(strand(sgrbam)))
@
}
  \end{block}
  
\end{frame}


\begin{frame}[fragile]
  \frametitle{Practical examples: Ura3 and TIM9}
  Use findOverlaps to see which things match
  \begin{block}{Example}
{\footnotesize
<<GetAnnotations3, results=hide>>=
OL <- findOverlaps(sgrbam, txpair)
head(matchMatrix(OL))
ind <- matchMatrix(OL)[,"query"]

## Subset to just the data that overlapped
subData <- grbam[ind]
seqlengths(subData) <- chr5Len
@
}
  \end{block}
  
\end{frame}


\begin{frame}[fragile]
  \frametitle{Practical examples: Ura3 and TIM9}
  We can also do this process in reverse (starting with ranges)
  \begin{block}{Example}
{\footnotesize
<<GetAnnotations4, results=hide>>=
##starting with the data object
subData

## call findOverlaps()
OLrev <- findOverlaps(subData, txs)
matchMatrix(OLrev)

## Then extract the genes that overlapped
ind <- unique(matchMatrix(OLrev)[,"subject"])
txs[ind]
@
}
  \end{block}
  
\end{frame}




%EXERCISE 3
\begin{frame}[fragile]
  \frametitle{GenomicFeatures exercise 3}
  Use the code below to load the results from the DE data from Martins
  talk earlier in the day. Use your gene centric annotation skills to
  find any gene names for the top 6 gene IDs, and then use your new
  skills to find the ranges that correspond to these genes.
{\small
<<exer3a, results=hide>>=
load(system.file("data", "nbTopTable.rda",
                 package = "SeattleIntro2010"))
@
}  
\end{frame}

\begin{frame}[fragile]
\end{frame}

%%TODO: find out what the results are from Martin (object) and then be
%%sure I can load that here and finish the exercise!
\begin{frame}[fragile]
  \frametitle{GenomicFeatures exercise 3 solution}
{\small
<<exer3b, results=hide>>=
## 1st lets see if any of these genes is even named
ids = head(nbTopTable)[,1]
library(org.Sc.sgd.db)
mget(ids, org.Sc.sgdGENENAME, ifnotfound=NA)

## Now lets get the ranges
ind = names(txs) %in% ids
txSub = txs[ind]
@
}
\end{frame}





\section{visualizing results}
%% Demonstration of how we can display the part of our data that maps to our gene of interest in a genome Browser

\begin{frame}[fragile]
  \frametitle{Visualizing the Data}
  \begin{itemize}
    \item Working with rtracklayer often involves wrapping things by using the 
  \Rfunction{GRangesForUCSCGenome} wrapper
    \item Then you can launch a genome browser from R
  \end{itemize}
  \begin{block}{Example}
{\footnotesize
<<BrowseData1, results=hide>>=
library(rtracklayer)
## construct a range that describes where you want to look
range <- GRangesForUCSCGenome(genome="sacCer2", chrom="chrV", 
                              ranges=ranges(subData),
                              strand = strand(subData))
# browseGenome(range = range)
@
}
  \end{block}
\end{frame}




\begin{frame}[fragile]
  \frametitle{Visualizing the Data}
  \begin{itemize}
    \item You can also save files as .bed .gff or .wig files
    \item Those files can then be imported into the genome browser
  \end{itemize}
  \begin{block}{Example}
{\footnotesize
<<BrowseData2, results=hide>>=
## You can export your tracks (bed, gff and wig are supported)
export( range, "GeneRange.bed")
## We also saved this for you 
loadedBedFile <- system.file("extdata", "GeneRange.bed",
                 package = "SeattleIntro2010")
loadedBedFile
## Can upload this into genome browser 
@
}
  \end{block}
\end{frame}




\begin{frame}[fragile]
  \frametitle{Visualizing the Data}
  Lets look at another gene (the last one had two genes that overlapped)
  \begin{block}{Example}
{\footnotesize
<<BrowseData3, results=hide>>=
ind = names(txs) %in% "YEL022W"
tx = txs[ind]
OL <- findOverlaps(grbam, tx)
ind <- matchMatrix(OL)[,"query"]
subDat <- grbam[ind]
@
}
  \end{block}
  
\end{frame}





\begin{frame}[fragile]
  \frametitle{Visualizing the Data} 
  \Rfunction{coverage} returns an rle with values computed to fit within the 
  entire chromosome.  Therefore we have to adjust for that.
  \begin{block}{Example}
{\footnotesize
<<BrowseData4, results=hide>>=
## a shift solution
cov <- coverage(subData, shift=list(chrV=1-min(start(subData))))
x <- cov[["chrV"]]

## or a slice solution
x <- slice(coverage(subData), 1L)[["chrV"]][[1]]
@
}
  \end{block}
  
\end{frame}




\begin{frame}[fragile]
  \frametitle{Visualizing the Data}
  \begin{itemize}
    \item This just defines a plotting function so we can look at the data
    \item Both the start and the runValues of the rle have to be plotted
  \end{itemize}
  \begin{block}{Example}
    {\footnotesize
<<BrowseData5, results=hide>>=
plotCoverage <- function(x, xlab = "Position", ylab = "Coverage",
 ...) {
 plot(c(start(x)), c(runValue(x)), type = "s", col = "blue", 
      xlab = xlab, ylab = ylab, ...)
}

## then call out plot function
plotCoverage(x)
@
}
  \end{block}
  
\end{frame}



\begin{frame}[fragile]
  \frametitle{Visualizing the Data}  
  Of course we also could have also plotted this using UCSC's genome browser
  \begin{block}{Example}
    {\footnotesize
<<BrowseData6, results=hide>>=
sd <- GRangesForUCSCGenome(genome="sacCer2", chrom="chrV", 
                           ranges=ranges(subData),
                           strand = strand(subData))
export(sd, "subData.bed")
@
}
  \end{block}
  
\end{frame}






%EXERCISE 4
\begin{frame}[fragile]
  \frametitle{GenomicFeatures exercise 4}
  Use the data provided to look up data reads that map to
  LEU2(YCL018W) and NFS1 (YCL017C) from chromosome 3.  Then either
  plot them or display the data in a browser.  What do you notice
  about these two genes?
\end{frame}

\begin{frame}[fragile]
\end{frame}


\begin{frame}[fragile]
  \frametitle{GenomicFeatures exercise 4 solution}
{\small
<<exer4, results=hide>>=
## Get annots
ind = names(txs) %in% c("YCL018W","YCL017C")
txpair = txs[ind]

## Overlap and subset
OL <- findOverlaps(grbam, txpair)
ind <- matchMatrix(OL)[,"query"]
subData <- grbam[ind]
seqlengths(subData) <- chr3Len

## Save as a track:
range <- GRangesForUCSCGenome(genome="sacCer2", chrom="chrIII", 
                              ranges=ranges(subData),
                              strand = strand(subData))
export( range, "GeneRange2.bed")

## or just plot
x <- slice(coverage(subData), 1L)[["chrIII"]][[1]]
plotCoverage(x)
@
}
\end{frame}



\section{Getting other resources from rtracklayer}

\begin{frame}[fragile]
  \frametitle{Retrieving  other data from rtracklayer}
  \begin{block}{Next release watch for FeatureDb objects}
    But for now we can use rtracklayer
  \end{block}
{\small
<<retrieveTracksAsdataframes, results=hide>>=
## create a session object
session <- browserSession()
## then ucscGenomes will list the genomes
head(ucscGenomes())
## Choose one, and assign it into your session object
genome(session) <- "sacCer2"
## NOW you can list the trackNames
head(trackNames(session))
@
}
\end{frame}



\begin{frame}[fragile]
  \frametitle{Retrieving other data from rtracklayer}
  To complete the process, we must create a query and then call \Rfunction{getTable}
{\small
<<retrieveTracksAsdataframes, results=hide>>=
## this information is used to create a query object
query <- ucscTableQuery(session, "oreganno")
## and for each track, you can list the tables it contains
tableNames(query)
## This information can be used to refine/modify your query object
tableName(query) <- "oreganno"
## Which is all then used to getTable to retrieve the data
result <- getTable(query)
head(result)
@
}
\end{frame}




\begin{frame}[fragile]
  \frametitle{Retrieving other data from rtracklayer}
  The final step is that you have to take whatever data you got from
  your query and put it into a GRanges object that can be used with
  findOverlaps etc.
{\small
<<retrieveTracksAsdataframes, results=hide>>=
ranges <- GRanges(seqnames = result[["chrom"]],
                  ranges = IRanges(start = result[["chromStart"]],
                                     end = result[["chromEnd"]]),
                  strand = result[["strand"]],
                  name = result[["name"]])
@
}
\end{frame}


\begin{frame}[fragile]
  \frametitle{GenomicFeatures exercise 5}
  Use rtracklayer to get the track for ensembl genes for yeast. Then
  use this information along with the reads that mapped to the about
  the NFS1(YCL017C) gene to determine the coordinates of any ensembl
  genes are in the same location as the NFS1 gene.
\end{frame}

\begin{frame}[fragile]
\end{frame}

\begin{frame}[fragile]
  \frametitle{GenomicFeatures exercise 5 solution}
{\small
<<exer5, results=hide>>=
query <- ucscTableQuery(session, "ensGene")
tableNames(query)
tableName(query) <- "ensGene"
ensGene <- getTable(query)

## then make that into a GRanges object.
ensRanges <- GRanges(seqnames = ensGene[["chrom"]],
                     ranges = IRanges(start = ensGene[["txStart"]],
                                        end = ensGene[["txEnd"]]),
                     strand = ensGene[["strand"]],
                     ensId = ensGene[["name"]])

@
}
\end{frame}





\end{document}




%% TODO: pretty up some of the older slides (split data off to it's own slide when appropriate and shrink it so it fits better!)







