### R code from vignette source 'vignettes/gmapR/inst/doc/gmapR.Rnw' ################################################### ### code chunk number 1: 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) ## ##> gmapGenome ## ##GmapGenome object ## ##genome: dm3 ## ##directory: /reshpcfs/home/coryba/projects/gmapR2/testGenome ################################################### ### code chunk number 2: 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 ################################################### ### code chunk number 3: 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 = "TP53_demo", create = TRUE) ################################################### ### code chunk number 4: get_lung_cancer_fastqs ################################################### library("LungCancerLines") fastqs <- LungCancerFastqFiles() ################################################### ### code chunk number 5: 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) ################################################### ### code chunk number 6: 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 ## ## ################################################### ### code chunk number 7: 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" ################################################### ### code chunk number 8: TP53Genome_accessor ################################################### genome <- TP53Genome() ################################################### ### code chunk number 9: 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 ## param <- BamTallyParam(genome, cycle_breaks = breaks, ## 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) ################################################### ### code chunk number 10: create_GmapGenomePackge (eval = FALSE) ################################################### ## makeGmapGenomePackage(gmapGenome=gmapGenome, ## version="1.0.0", ## maintainer="", ## author="Your Name", ## destDir="myDestDir", ## license="Artistic-2.0", ## pkgName="GmapGenome.Dmelanogaster.UCSC.dm3") ################################################### ### code chunk number 11: 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="", ## author="Your Name", ## destDir=destDir, ## license="Artistic-2.0", ## pkgName="GmapGenome.Hsapiens.UCSC.hg19") ################################################### ### code chunk number 12: SessionInfo ################################################### sessionInfo()