@@ -125,8 +125,9 @@ example demonstrates:

 library(gmapR)

 if (!suppressWarnings(require(BSgenome.Dmelanogaster.UCSC.dm3))) {

-  source("http://bioconductor.org/biocLite.R")

-  biocLite("BSgenome.Dmelanogaster.UCSC.dm3")

+  if (!requireNamespace("BiocManager", quietly=TRUE))

+      install.packages("BiocManager")

+  BiocManager::install("BSgenome.Dmelanogaster.UCSC.dm3")

   library(BSgenome.Dmelanogaster.UCSC.dm3)

 }

@@ -348,8 +349,9 @@ object:

 <<create_hg19_GmapGenome, eval=FALSE>>=

 if (!suppressWarnings(require(BSgenome.Hsapiens.UCSC.hg19))) {

-  source("http://bioconductor.org/biocLite.R")

-  biocLite("BSgenome.Hsapiens.UCSC.hg19")

+  if (!requireNamespace("BiocManager", quietly=TRUE))

+      install.packages("BiocManager")

+  BiocManager::install("BSgenome.Hsapiens.UCSC.hg19")

   library(BSgenome.Hsapiens.UCSC.hg19)

 }

 gmapGenome <- GmapGenome(genome=Hsapiens,

@@ -308,10 +308,11 @@ param <- BamTallyParam(genome,

                        primary_only = FALSE,

                        min_depth = 0L, variant_strand = 1L,

                        ignore_query_Ns = TRUE,

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

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

+                       min_base_quality = 23L)

 tallies <-bam_tally(bam_file,

                     param)

-variantSummary(tallies, high_base_quality = 23L)

+variantSummary(tallies)

 @ 

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

@@ -165,7 +165,7 @@ 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, vals = list(gene_id = eg))

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

 roi <- range(tx) + 1e6

 strand(roi) <- "*"

 @ 

@@ -308,7 +308,7 @@ param <- BamTallyParam(genome,

                        primary_only = FALSE,

                        min_depth = 0L, variant_strand = 1L,

                        ignore_query_Ns = TRUE,

-                       indels = FALSE, include_soft_clips = 1L)

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

 tallies <-bam_tally(bam_file,

                     param)

 variantSummary(tallies, high_base_quality = 23L)

@@ -308,7 +308,7 @@ param <- BamTallyParam(genome,

                        primary_only = FALSE,

                        min_depth = 0L, variant_strand = 1L,

                        ignore_query_Ns = TRUE,

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

+                       indels = FALSE, include_soft_clips = 1L)

 tallies <-bam_tally(bam_file,

                     param)

 variantSummary(tallies, high_base_quality = 23L)

@@ -308,7 +308,7 @@ param <- BamTallyParam(genome,

                        primary_only = FALSE,

                        min_depth = 0L, variant_strand = 1L,

                        ignore_query_Ns = TRUE,

-                       indels = FALSE, include_soft_clips = 1L)

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

 tallies <-bam_tally(bam_file,

                     param)

 variantSummary(tallies, high_base_quality = 23L)

@@ -308,7 +308,7 @@ param <- BamTallyParam(genome,

                        primary_only = FALSE,

                        min_depth = 0L, variant_strand = 1L,

                        ignore_query_Ns = TRUE,

-                       indels = FALSE, include_soft_clips = 1L, count_xs = TRUE)

+                       indels = FALSE, include_soft_clips = 1L)

 tallies <-bam_tally(bam_file,

                     param)

 variantSummary(tallies, high_base_quality = 23L)

@@ -308,7 +308,7 @@ param <- BamTallyParam(genome,

                        primary_only = FALSE,

                        min_depth = 0L, variant_strand = 1L,

                        ignore_query_Ns = TRUE,

-                       indels = FALSE, include_soft_clips = 1L)

+                       indels = FALSE, include_soft_clips = 1L, count_xs = TRUE)

 tallies <-bam_tally(bam_file,

                     param)

 variantSummary(tallies, high_base_quality = 23L)

@@ -308,7 +308,7 @@ param <- BamTallyParam(genome,

                        primary_only = FALSE,

                        min_depth = 0L, variant_strand = 1L,

                        ignore_query_Ns = TRUE,

-                       indels = FALSE)

+                       indels = FALSE, include_soft_clips = 1L)

 tallies <-bam_tally(bam_file,

                     param)

 variantSummary(tallies, high_base_quality = 23L)

@@ -182,9 +182,12 @@ library("gmapR")

 p53Seq <- getSeq(BSgenome.Hsapiens.UCSC.hg19::Hsapiens, roi,

                  as.character = FALSE) 

 names(p53Seq) <- "TP53"

-gmapGenome <- GmapGenome(genome = p53Seq, name = "TP53_demo", create = TRUE, 

+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>>=

@@ -294,11 +297,8 @@ 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>>=

-library(gmapR)

-

 bam_file <- system.file("extdata/H1993.analyzed.bam", 

                         package="LungCancerLines", mustWork=TRUE)

-

 breaks <- c(0L, 15L, 60L, 75L)

 bqual <- 56L

 mapq <- 13L

@@ -215,7 +215,7 @@ gsnapParam <- GsnapParam(genome=gmapGenome,

                          splicing="knownGene",

                          indel_penalty=1L,

                          distant_splice_penalty=1L,

-                         clip_overlaps=TRUE)

+                         clip_overlap=TRUE)

 @ 

 Now we are ready to align.

@@ -137,9 +137,9 @@ gmapGenomeDirectory <- GmapGenomeDirectory(gmapGenomePath, create = TRUE)

 ##path: /reshpcfs/home/coryba/projects/gmapR2/testGenome 

 gmapGenome <- GmapGenome(genome=Dmelanogaster,

-                         directory = gmapGenomeDirectory,

-                         name = "dm3",

-                         create = TRUE,

+                         directory=gmapGenomeDirectory,

+                         name="dm3",

+                         create=TRUE,

                          k = 12L)

 ##> gmapGenome

 ##GmapGenome object

@@ -167,6 +167,7 @@ eg <- org.Hs.eg.db::org.Hs.egSYMBOL2EG[["TP53"]]

 txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

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

 roi <- range(tx) + 1e6

+strand(roi) <- "*"

 @ 

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

@@ -181,7 +182,14 @@ library("gmapR")

 p53Seq <- getSeq(BSgenome.Hsapiens.UCSC.hg19::Hsapiens, roi,

                  as.character = FALSE) 

 names(p53Seq) <- "TP53"

-gmapGenome <- GmapGenome(genome = p53Seq, name = "TP53_demo", create = TRUE, k = 12L)

+gmapGenome <- GmapGenome(genome = p53Seq, name = "TP53_demo", 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

@@ -199,15 +207,15 @@ behavior.

 <<create_gsnapParam>>=

 ##specify how GSNAP should behave using a GsnapParam object

-gsnapParam <- GsnapParam(genome = gmapGenome,

-                         unique_only = FALSE,

-                         max_mismatches = NULL,

-                         suboptimal_levels = 0, mode = "standard",

-                         npaths = 10,

-                         novelsplicing = FALSE, splicing = NULL, 

-                         nthreads = 1,

-                         batch = 2L,

-                         gunzip=TRUE)

+gsnapParam <- GsnapParam(genome=gmapGenome,

+                         unique_only=FALSE,

+                         suboptimal_levels=2L, 

+                         npaths=10L,

+                         novelsplicing=TRUE,

+                         splicing="knownGene",

+                         indel_penalty=1L,

+                         distant_splice_penalty=1L,

+                         clip_overlaps=TRUE)

 @ 

 Now we are ready to align.

@@ -216,6 +224,7 @@ Now we are ready to align.

 gsnapOutput <- gsnap(input_a=fastqs["H1993.first"],

                      input_b=fastqs["H1993.last"],

                      params=gsnapParam)

+

 ##gsnapOutput

 ##An object of class "GsnapOutput"

 ##Slot "path":

@@ -282,7 +282,7 @@ 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 \code{VRanges} object

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

 describing putative genetic variants in the sample.

 <<run_bamtally, eval=FALSE>>=

 library(gmapR)

@@ -279,6 +279,11 @@ demonstrate the use of \software{bam\_tally} through the

 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 \code{VRanges} object

+describing putative genetic variants in the sample.

 <<run_bamtally, eval=FALSE>>=

 library(gmapR)

@@ -288,8 +293,7 @@ bam_file <- system.file("extdata/H1993.analyzed.bam",

 breaks <- c(0L, 15L, 60L, 75L)

 bqual <- 56L

 mapq <- 13L

-param <- BamTallyParam(genome, read_pos_breaks = breaks,

-                       high_base_quality = bqual,

+param <- BamTallyParam(genome,

                        minimum_mapq = mapq,

                        concordant_only = FALSE, unique_only = FALSE,

                        primary_only = FALSE,

@@ -298,6 +302,7 @@ param <- BamTallyParam(genome, read_pos_breaks = breaks,

                        indels = FALSE)

 tallies <-bam_tally(bam_file,

                     param)

+variantSummary(tallies, high_base_quality = 23L)

 @ 

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

@@ -277,7 +277,7 @@ 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 nucleotide cycle for each allele.

+scores, and read position for each allele.

 <<run_bamtally, eval=FALSE>>=

 library(gmapR)

@@ -288,7 +288,7 @@ bam_file <- system.file("extdata/H1993.analyzed.bam",

 breaks <- c(0L, 15L, 60L, 75L)

 bqual <- 56L

 mapq <- 13L

-param <- BamTallyParam(genome, cycle_breaks = breaks,

+param <- BamTallyParam(genome, read_pos_breaks = breaks,

                        high_base_quality = bqual,

                        minimum_mapq = mapq,

                        concordant_only = FALSE, unique_only = FALSE,

@@ -181,7 +181,7 @@ library("gmapR")

 p53Seq <- getSeq(BSgenome.Hsapiens.UCSC.hg19::Hsapiens, roi,

                  as.character = FALSE) 

 names(p53Seq) <- "TP53"

-gmapGenome <- GmapGenome(genome = p53Seq, name = "TP53_demo", create = TRUE)

+gmapGenome <- GmapGenome(genome = p53Seq, name = "TP53_demo", create = TRUE, k = 12L)

 @ 

 The data package \software{LungCancerLines} contains fastqs of reads

@@ -139,7 +139,8 @@ gmapGenomeDirectory <- GmapGenomeDirectory(gmapGenomePath, create = TRUE)

 gmapGenome <- GmapGenome(genome=Dmelanogaster,

                          directory = gmapGenomeDirectory,

                          name = "dm3",

-                         create = TRUE)

+                         create = TRUE,

+                         k = 12L)

 ##> gmapGenome

 ##GmapGenome object

 ##genome: dm3 

@@ -169,7 +169,7 @@ roi <- range(tx) + 1e6

 @ 

 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

+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.)

@@ -273,7 +273,7 @@ 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

+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 nucleotide cycle for each allele.

@@ -176,6 +176,7 @@ 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"

@@ -130,7 +130,7 @@ if (!suppressWarnings(require(BSgenome.Dmelanogaster.UCSC.dm3))) {

   library(BSgenome.Dmelanogaster.UCSC.dm3)

 }

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

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

 gmapGenomeDirectory <- GmapGenomeDirectory(gmapGenomePath, create = TRUE)

 ##> gmapGenomeDirectory

 ##GmapGenomeDirectory object

@@ -151,41 +151,51 @@ gmapGenome <- GmapGenome(genome=Dmelanogaster,

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 The \software{GSNAP} algorithm incorporates biological knowledge to

-provide accurate alignments, particularly for RNA-seq data. The

-following example synthesizes some input reads and aligns them using

-\software{GSNAP}.

-

-If you have created the D. melanogaster \Rclass{GmapGenome} package as

-detailed above and installed it, you can begin this section as

-follows:

+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, vals = list(gene_id = eg))

+roi <- range(tx) + 1e6

+@ 

-<<<load_GmapGenome, eval=FALSE>>=

-library("GmapGenome.Dmelanogaster.UCSC.dm3")

-gmapGenome <- GmapGenome.Dmelanogaster.UCSC.dm3

+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")

+p53Seq <- getSeq(BSgenome.Hsapiens.UCSC.hg19::Hsapiens, roi,

+                 as.character = FALSE) 

+names(p53Seq) <- "TP53"

+gmapGenome <- GmapGenome(genome = p53Seq, name = "TP53_demo", create = TRUE)

 @ 

-Otherwise, you can create the \Rcode{gmapGenome} object as detailed in the

-initial section of this vignette.

+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}.

-<<align_with_gsnap, eval=FALSE>>=

-##synthesize paired-end FASTA input to align against the genome

-library("BSgenome.Dmelanogaster.UCSC.dm3")

-dnaStringSet <- Dmelanogaster$chr4

-makeTestFasta <- function(seqStarts, pairNum) {

-  seqEnds <- seqStarts + 74

-  nucleotides <- sapply(seq_along(seqStarts),

-                        function(i) as.character(subseq(dnaStringSet, seqStarts[i], seqEnds[i]))

-                        )

-  fastaHeaders <- paste0(">", "fakeRead", seq_along(nucleotides), "#0/", pairNum)

-  fastaLines <- paste(fastaHeaders, nucleotides, sep="\n")

-  tmpFasta <- paste0(tempfile(), ".", pairNum, ".fa")

-  writeLines(fastaLines, con=tmpFasta)

-  tmpFasta

-}

-seqStarts <- seq(1, 100000, by = 10000) + length(dnaStringSet) * 2 / 3

-tmpFasta1 <- makeTestFasta(seqStarts, '1')

-tmpFasta2 <- makeTestFasta(seqStarts + 125, '2')

+<<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,

@@ -194,10 +204,15 @@ gsnapParam <- GsnapParam(genome = gmapGenome,

                          npaths = 10,

                          novelsplicing = FALSE, splicing = NULL, 

                          nthreads = 1,

-                         batch = 2L)

+                         batch = 2L,

+                         gunzip=TRUE)

+@ 

-gsnapOutput <- gsnap(input_a=tmpFasta1,

-                     input_b=tmpFasta2,

+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"

@@ -244,29 +259,41 @@ alignments that \software{GSNAP} outputs along with other utilities.

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 Running the \Rmethod{bam\_tally} method will return a \Rclass{GRanges}

-of information per nucleotide. Below is a basic example. See the

-documentation for the \Rmethod{bam\_tally} method for details.

+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 nucleotide cycle for each allele.

 <<run_bamtally, eval=FALSE>>=

 library(gmapR)

-library(GmapGenome.Hsapiens.UCSC.hg19)

-genome <- GmapGenome.Hsapiens.UCSC.hg19

-bam_file <- file.path(system.file(package = "gmapR", mustWork=TRUE),

-                      "extdata/bam_tally_test.bam")

+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, cycle_breaks = breaks,

-                       high_quality_cutoff = bqual,

+                       high_base_quality = bqual,

                        minimum_mapq = mapq,

                        concordant_only = FALSE, unique_only = FALSE,

                        primary_only = FALSE,

                        min_depth = 0L, variant_strand = 1L,

                        ignore_query_Ns = TRUE,

                        indels = FALSE)

-

 tallies <-bam_tally(bam_file,

                     param)

 @ 

@@ -330,26 +357,28 @@ 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:

+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

+%\section{Aligning with SNP Tolerance}

+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

-<<<load_GmapGenomeHsapiens>>=

-library("GmapGenome.Hsapiens.UCSC.hg19")

-@ 

+%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()

@@ -155,8 +155,9 @@ provide accurate alignments, particularly for RNA-seq data. The

 following example synthesizes some input reads and aligns them using

 \software{GSNAP}.

-If you've created the D. melanogaster \Rclass{GmapGenome} package as detailed above and installed

-it, you can begin this section as follows:

+If you have created the D. melanogaster \Rclass{GmapGenome} package as

+detailed above and installed it, you can begin this section as

+follows:

 <<<load_GmapGenome, eval=FALSE>>=

 library("GmapGenome.Dmelanogaster.UCSC.dm3")

@@ -189,11 +190,11 @@ tmpFasta2 <- makeTestFasta(seqStarts + 125, '2')

 gsnapParam <- GsnapParam(genome = gmapGenome,

                          unique_only = FALSE,

                          max_mismatches = NULL,

-                         suboptimal_levels = 0L, mode = "standard",

-                         npaths = 10L,

+                         suboptimal_levels = 0, mode = "standard",

+                         npaths = 10,

                          novelsplicing = FALSE, splicing = NULL, 

-                         nthreads = 1L,

-                         batch = "2")

+                         nthreads = 1,

+                         batch = 2L)

 gsnapOutput <- gsnap(input_a=tmpFasta1,

                      input_b=tmpFasta2,

@@ -289,7 +289,7 @@ makeGmapGenomePackage(gmapGenome=gmapGenome,

                       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

@@ -295,6 +295,61 @@ 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()

 @ 

@@ -1,3 +1,4 @@

+

 % \VignetteIndexEntry{gmapR}

 % \VignetteKeywords{gmap,gsnap}

 % \VignettePackage{gmapR}

@@ -256,9 +257,9 @@ bam_file <- file.path(system.file(package = "gmapR", mustWork=TRUE),

 breaks <- c(0L, 15L, 60L, 75L)

 bqual <- 56L

 mapq <- 13L

-param <- BamTallyParam(cycle_breaks = as.integer(breaks),

-                       high_quality_cutoff = as.integer(bqual),

-                       minimum_mapq = as.integer(mapq),

+param <- BamTallyParam(genome, cycle_breaks = breaks,

+                       high_quality_cutoff = bqual,

+                       minimum_mapq = mapq,

                        concordant_only = FALSE, unique_only = FALSE,

                        primary_only = FALSE,

                        min_depth = 0L, variant_strand = 1L,

@@ -266,7 +267,6 @@ param <- BamTallyParam(cycle_breaks = as.integer(breaks),

                        indels = FALSE)

 tallies <-bam_tally(bam_file,

-                    genome,

                     param)

 @ 

new file mode 100644

@@ -0,0 +1,301 @@

+% \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(), "testGenome")

+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)

+##> 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. The

+following example synthesizes some input reads and aligns them using

+\software{GSNAP}.

+

+If you've created the D. melanogaster \Rclass{GmapGenome} package as detailed above and installed

+it, you can begin this section as follows:

+

+<<<load_GmapGenome, eval=FALSE>>=

+library("GmapGenome.Dmelanogaster.UCSC.dm3")

+gmapGenome <- GmapGenome.Dmelanogaster.UCSC.dm3

+@ 

+

+Otherwise, you can create the \Rcode{gmapGenome} object as detailed in the

+initial section of this vignette.

+

+<<align_with_gsnap, eval=FALSE>>=

+##synthesize paired-end FASTA input to align against the genome

+library("BSgenome.Dmelanogaster.UCSC.dm3")

+dnaStringSet <- Dmelanogaster$chr4

+makeTestFasta <- function(seqStarts, pairNum) {

+  seqEnds <- seqStarts + 74

+  nucleotides <- sapply(seq_along(seqStarts),

+                        function(i) as.character(subseq(dnaStringSet, seqStarts[i], seqEnds[i]))

+                        )

+  fastaHeaders <- paste0(">", "fakeRead", seq_along(nucleotides), "#0/", pairNum)

+  fastaLines <- paste(fastaHeaders, nucleotides, sep="\n")

+  tmpFasta <- paste0(tempfile(), ".", pairNum, ".fa")

+  writeLines(fastaLines, con=tmpFasta)

+  tmpFasta

+}

+seqStarts <- seq(1, 100000, by = 10000) + length(dnaStringSet) * 2 / 3

+tmpFasta1 <- makeTestFasta(seqStarts, '1')

+tmpFasta2 <- makeTestFasta(seqStarts + 125, '2')

+

+##specify how GSNAP should behave using a GsnapParam object

+gsnapParam <- GsnapParam(genome = gmapGenome,

+                         unique_only = FALSE,

+                         max_mismatches = NULL,