% \VignetteIndexEntry{gmapR}
 % \VignetteKeywords{gmap,gsnap}
 % \VignettePackage{gmapR}
 \documentclass[10pt]{article}
 
 \usepackage{times}
 \usepackage{hyperref}
 \usepackage{Sweave}
 
 \textwidth=6.5in
 \textheight=8.5in
 % \parskip=.3cm
 \oddsidemargin=-.1in
 \evensidemargin=-.1in
 \headheight=-.3in
 
 \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}}
 
 \title{gmapR: Use the GMAP Suite of Tools in R}
 \author{Michael Lawrence, Cory Barr}
 \date{\today}
 
 \begin{document}
 
 \maketitle
 \tableofcontents
 \newpage
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Introduction}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 The \Rpackage{gmapR} packages provides users with a way to access
 \software{GSNAP}, \software{bam\_tally}, and other utilities from the
 \software{GMAP} suite of tools within an R session. In this vignette,
 we briefly look at the \software{GMAP} suite of tools available
 through the \Rpackage{gmapR} package and work through an example.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{What is \software{GMAP}, \software{GSNAP}, and \software{bam\_tally}?}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 The \software{GMAP} suite offers useful tools for the following:
 
 \begin{itemize}
 
   \item Genomic mapping: Given a cDNA, find where it best aligns to
     an entire genome
 
   \item Genomic alignment: Given a cDNA and a genomic segment,
     provide a nucleotide-level correspondence for the exons of the
     cDNA to the genomic segment
 
   \item Summarization via coverage plus reference and allele
     nucleotide polymorphism counts for an aligned set of sequencing
     reads over a given genomic location
 
 \end{itemize}
 
 \software{GMAP} (Genomic Mapping and Alignment Program) is
 particularly suited to relatively long mRNA and EST sequences such as
 those that are obtained from Roche 454 or Pacific Biosciences
 sequencing technologies. (At present, only \software{GSNAP} is
 available through the \Rpackage{gmapR}. \software{GMAP} integration is
 scheduled for the near future.)
 
 \software{GSNAP} (Genomic Short-read Nucleotide Alignment Program)
 also provides users with genomic mapping and alignment capabilities,
 but is optimized to handle issues that arise when dealing with the
 alignment of short reads generated from sequencing technologies such
 as those from Illumina/Solexa or ABI/SOLiD. \software{GSNAP} offers
 the following functionality, as mentioned in
 \href{http://bioinformatics.oxfordjournals.org/content/26/7/873.full}{Fast
   and SNP-tolerant detection of complex variants and splicing in short
   reads} by Thomas D. Wu and Serban Nacu:
 
 \begin{itemize}
   
   \item fast detection of complex variants and splicing in short
     reads, based on a successively constrained search process of
     merging and filtering position lists from a genomic index
 
   \item alignment of both single- and paired-end reads as short as 14
     nt and of arbitrarily long length
     
   \item detection of short- and long-distance splicing, including
     interchromosomal splicing, in individual reads, using
     probabilistic models or a database of known splice sites
     
   \item SNP-tolerant alignment to a reference space of all possible
     combinations of major and minor alleles
 
   \item alignment of reads from bisulfite-treated DNA for the study of
     methylation state
     
 \end{itemize}
     
 \software{bam\_tally} provides users with the coverage as well as
 counts of both reference alleles, alternative alleles, and indels over
 a genomic region.
 
 For more detailed information on the \software{GMAP} suite of tools
 including a detailed explication of algorithmic specifics, see
 \url{http://research-pub.gene.com/gmap/}.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Create a \Rclass{GmapGenome} Object}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 To align reads with \software{GMAP} or \software{GSNAP}, or to use
 \software{bam\_tally}, you will need to either obtain or create a
 \Rclass{GmapGenome} object. \Rclass{GmapGenome} objects can be created
 from FASTA files or \Rclass{BSgenome} objects, as the following
 example demonstrates:
 
 <<create_GmapGenome, eval=FALSE>>=
 library(gmapR)
 
 if (!suppressWarnings(require(BSgenome.Dmelanogaster.UCSC.dm3))) {
   source("http://bioconductor.org/biocLite.R")
   biocLite("BSgenome.Dmelanogaster.UCSC.dm3")
   library(BSgenome.Dmelanogaster.UCSC.dm3)
 }
 

 gmapGenomePath <- file.path(getwd(), "flyGenome")

 gmapGenomeDirectory <- GmapGenomeDirectory(gmapGenomePath, create = TRUE)
 ##> gmapGenomeDirectory
 ##GmapGenomeDirectory object
 ##path: /reshpcfs/home/coryba/projects/gmapR2/testGenome 
 
 gmapGenome <- GmapGenome(genome=Dmelanogaster,

                          directory=gmapGenomeDirectory,
                          name="dm3",
                          create=TRUE,

                          k = 12L)

 ##> gmapGenome
 ##GmapGenome object
 ##genome: dm3 
 ##directory: /reshpcfs/home/coryba/projects/gmapR2/testGenome 
 @
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Aligning with \software{GSNAP}}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 The \software{GSNAP} algorithm incorporates biological knowledge to

 provide accurate alignments, particularly for RNA-seq data. In this
 section, we will align reads from an RNA-seq experiment provided in a
 fastq file to a selected region of the human genome. In this example,
 we will align to the region of the human genome containing TP53 plus
 an additional one megabase on each side of this gene.
 
 First we need to obtain the desired region of interest from the genome:
 
 <<get_TP53_coordinates>>=
 library("org.Hs.eg.db")
 library("TxDb.Hsapiens.UCSC.hg19.knownGene")
 eg <- org.Hs.eg.db::org.Hs.egSYMBOL2EG[["TP53"]]
 txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

 tx <- transcripts(txdb, filter = list(gene_id = eg))

 roi <- range(tx) + 1e6

 strand(roi) <- "*"

 @ 

 Next we get the genetic sequence and use it to create a GmapGenome

 object. (Please note that the TP53\_demo GmapGenome object is used by

 many examples and tests in the gmapR and VariationTools packages. If
 the object has been created before, its subsequent creation will be
 instantaneous.)
 
 <<get_TP53_seq>>=
 library("BSgenome.Hsapiens.UCSC.hg19")

 library("gmapR")

 p53Seq <- getSeq(BSgenome.Hsapiens.UCSC.hg19::Hsapiens, roi,
                  as.character = FALSE) 
 names(p53Seq) <- "TP53"

 gmapGenome <- GmapGenome(genome = p53Seq, 
                          name = paste0("TP53_demo_", 
                            packageVersion("TxDb.Hsapiens.UCSC.hg19.knownGene")), 
                          create = TRUE, 

                          k = 12L)

 @

 
 We add the known transcripts (splice sites) to the genome index:
 <<set_TP53_splicesites>>=
 exons <- gmapR:::subsetRegion(exonsBy(txdb), roi, "TP53")
 spliceSites(gmapGenome, "knownGene") <- exons

 @ 
 

 The data package \software{LungCancerLines} contains fastqs of reads
 obtained from sequencing H1993 and H2073 cell lines. We will use these
 fastqs to demonstrate \software{GSNAP} alignment with \software{gmapR}.

 <<get_lung_cancer_fastqs>>=
 library("LungCancerLines")
 fastqs <- LungCancerFastqFiles()
 @ 

 \software{GSNAP} is highly configurable. Users create
 \Rclass{GsnapParam} objects to specify desired \software{GSNAP}
 behavior. 
 
 <<create_gsnapParam>>=

 ##specify how GSNAP should behave using a GsnapParam object

 gsnapParam <- GsnapParam(genome=gmapGenome,
                          unique_only=FALSE,
                          suboptimal_levels=2L, 
                          npaths=10L,
                          novelsplicing=TRUE,
                          splicing="knownGene",
                          indel_penalty=1L,
                          distant_splice_penalty=1L,

                          clip_overlap=TRUE)

 @ 

 Now we are ready to align.
 
 <<align_with_gsnap, eval=FALSE>>=
 gsnapOutput <- gsnap(input_a=fastqs["H1993.first"],
                      input_b=fastqs["H1993.last"],

                      params=gsnapParam)

 ##gsnapOutput
 ##An object of class "GsnapOutput"
 ##Slot "path":
 ##[1] "/local/Rtmporwsvr"
 ##
 ##Slot "param":
 ##A GsnapParams object
 ##genome: dm3 (/reshpcfs/home/coryba/projects/gmapR2/testGenome)
 ##part: NULL
 ##batch: 2
 ##max_mismatches: NULL
 ##suboptimal_levels: 0
 ##snps: NULL
 ##mode: standard
 ##nthreads: 1
 ##novelsplicing: FALSE
 ##splicing: NULL
 ##npaths: 10
 ##quiet_if_excessive: FALSE
 ##nofails: FALSE
 ##split_output: TRUE
 ##extra: list()
 ##
 ##Slot "version":
 ## [1] NA NA NA NA NA NA NA NA NA NA NA NA
 ##
 @ 
 
 The \Rcode{gsnapOutput} object will be of the class
 \Rclass{GsnapOutput}. It will provide access to the BAM files of
 alignments that \software{GSNAP} outputs along with other utilities.
 
 <<get_gsnap_bam_files, eval=FALSE>>=
 ##> dir(path(gsnapOutput), full.names=TRUE, pattern="\\.bam$")
 ##[1] "/local/Rtmporwsvr/file1cbc73503e9.1.nomapping.bam"        
 ##[2] "/local/Rtmporwsvr/file1cbc73503e9.1.unpaired_mult.bam"    
 ##[3] "/local/Rtmporwsvr/file1cbc73503e9.1.unpaired_transloc.bam"
 ##[4] "/local/Rtmporwsvr/file1cbc73503e9.1.unpaired_uniq.bam"    
 @ 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Using \Rmethod{bam\_tally}}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 Running the \Rmethod{bam\_tally} method will return a \Rclass{GRanges}

 of information per nucleotide. Below is an example demonstrating how
 to find variants in the TP53 gene of the human genome. See the
 documentation for the \Rmethod{bam\_tally} method for more details.
 
 \software{gmapR} provides access to a demo genome for examples. This
 genome encompasses the TP53 gene along with a 1-megabase flanking
 region on each side.
 <<TP53Genome_accessor>>=
 genome <- TP53Genome()
 @ 
 
 The \software{LungCancerLines} R package contains a BAM file of reads
 aligned to the TP53 region of the human genome. We'll use this file to

 demonstrate the use of \software{bam\_tally} through the

 \software{gmapR} package. The resulting data structure will contain
 the needed information such as number of alternative alleles, quality

 scores, and read position for each allele.

 The call to \Rfunction{bam\_tally} returns an opaque pointer to a
 C-level data structure. We anticipate having multiple means of
 summarizing these data into R data structures. For now, there is one:

 \Rfunction{variantSummary}, which returns a \Rcode{VRanges} object

 describing putative genetic variants in the sample.

 <<run_bamtally, eval=FALSE>>=

 bam_file <- system.file("extdata/H1993.analyzed.bam", 
                         package="LungCancerLines", mustWork=TRUE)

 breaks <- c(0L, 15L, 60L, 75L)
 bqual <- 56L
 mapq <- 13L

 param <- BamTallyParam(genome,

                        minimum_mapq = mapq,

                        concordant_only = FALSE, unique_only = FALSE,
                        primary_only = FALSE,
                        min_depth = 0L, variant_strand = 1L,
                        ignore_query_Ns = TRUE,

                        indels = FALSE, include_soft_clips = 1L, xs=TRUE)

 tallies <-bam_tally(bam_file,
                     param)

 variantSummary(tallies, high_base_quality = 23L)

 @ 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Creating a \Rclass{GmapGenome} Package}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 After creating a \Rclass{GmapGenome} object, you might want to
 distribute it for collaboration or version it for reproducible
 research. The function \Rcode{makeGmapGenomePackage} allows you to do
 this. Continuing on with the D. melanogaster example from above, here
 is how to archive the D. melanogaster \Rclass{GmapGenome} object in a
 \Rclass{GmapGenome} package:
 
 <<create_GmapGenomePackge, eval=FALSE>>=
 makeGmapGenomePackage(gmapGenome=gmapGenome,
                       version="1.0.0",
                       maintainer="<your.name@somewhere.com>",
                       author="Your Name",
                       destDir="myDestDir",
                       license="Artistic-2.0",
                       pkgName="GmapGenome.Dmelanogaster.UCSC.dm3")

 @

 
 After creating the package, you can run \Rcode{R CMD INSTALL
   myDestDir} to install it, or run \Rcode{R CMD build myDestDir} to
 create a distribution of the package.
 

 Many users with be working with the human genome. Many of the examples
 used in \software{gmapR} make use of a particular build of the human
 genome. As such, creating a \Rclass{GmapGenome} of hg19 is
 recommended. Here is one way to create it, using a \Rclass{BSgenome}
 object:
 
 <<create_hg19_GmapGenome, eval=FALSE>>=
 if (!suppressWarnings(require(BSgenome.Hsapiens.UCSC.hg19))) {
   source("http://bioconductor.org/biocLite.R")
   biocLite("BSgenome.Hsapiens.UCSC.hg19")
   library(BSgenome.Hsapiens.UCSC.hg19)
 }
 gmapGenome <- GmapGenome(genome=Hsapiens,
                          directory = "Hsapiens",
                          name = "hg19",
                          create = TRUE)
 destDir <- "HsapiensGmapGenome"
 pkgName <- "GmapGenome.Hsapiens.UCSC.hg19"
 makeGmapGenomePackage(gmapGenome=gmapGenome,
                       version="1.0.0",
                       maintainer="<your.name@somewhere.com>",
                       author="Your Name",
                       destDir=destDir,
                       license="Artistic-2.0",
                       pkgName="GmapGenome.Hsapiens.UCSC.hg19")
 @ 
 
 After running the above code, you should be able to run \Rcode{R CMD
   INSTALL GmapGenome.Hsapiens.UCSC.hg19} in the appropriate directory
 from the command line, submit GmapGenome.Hsapiens.UCSC.hg19 to a
 repository, etc. After installation,
 \Rcode{library("GmapGenome.Hsapiens.UCSC.hg19")} will load a
 \Rclass{GmapGenome} object that has the same name as the package.
 
 
 

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %\section{Aligning with SNP Tolerance}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %Both \software{GSNAP} and \software{GMAP} have the ability to perform
 %alignment without penalizing against a set of known SNPs. In this
 %section, we work through an example of aligning with this SNP
 %tolerance. By using SNP-tolerance, we'll perform alignments without
 %penalizing reads that match SNPs contained in the provided SNPs. In
 %this case, we'll use dbSNP to provide the SNPs.
 %
 %This example will use hg19. If you have not installed
 %\Rcode{GmapGenome.Hsapiens.UCSC.hg19}, please refer to the section
 %\ref{create_hg19_GmapGenome}.
 %
 %First we load the Hsapiens \Rclass{GmapGenome} object:
 %
 %<<<load_GmapGenomeHsapiens>>=
 %library("GmapGenome.Hsapiens.UCSC.hg19")
 %@ 

 <<SessionInfo>>=
 sessionInfo()
 @ 
 \end{document}