### R code from vignette source 'SigFuge.Rnw' ################################################### ### code chunk number 1: style ################################################### BiocStyle::latex() ################################################### ### code chunk number 2: SigFuge.Rnw:49-56 ################################################### library(SigFuge) geneAnnot <- GRanges(seqnames=Rle("chr9", 5), ranges=IRanges(start=c(21967751,21968574,21970901,21974403,21994138), end=c(21968241,21968770,21971207,21975132,21994490)), strand=Rle(strand("-"), 5)) geneAnnot ################################################### ### code chunk number 3: SigFuge.Rnw:61-73 ################################################### library(org.Hs.eg.db) library(TxDb.Hsapiens.UCSC.hg19.knownGene) genesym <- c("CDKN2A") geneid <- select(org.Hs.eg.db, keys=genesym, keytype="SYMBOL", columns="ENTREZID") txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene ex <- exonsBy(txdb, "gene") geneAnnot <- ex[[geneid$ENTREZID[1]]] geneAnnot geneAnnot <- reduce(geneAnnot) geneAnnot ################################################### ### code chunk number 4: SigFuge.Rnw:78-99 ################################################### library(Rsamtools) library(prebsdata) geneAnnot <- GRanges(seqnames=Rle("20", 9), ranges=IRanges( start=c(49551405,49552685,49557402,49558568,49562274, 49562384,49565166,49571723,49574900), end=c(49551773,49552799,49557492,49558663,49562299, 49562460,49565199,49571822,49575060)), strand=Rle(strand("-"), 9)) bam_file1 <- system.file(file.path("sample_bam_files", "input1.bam"), package="prebsdata") bam_file2 <- system.file(file.path("sample_bam_files", "input2.bam"), package="prebsdata") sorted1 <- sortBam(bam_file1, tempfile()) indexBam(sorted1) sorted2 <- sortBam(bam_file2, tempfile()) indexBam(sorted2) bam_files <- c(sorted1, sorted2) ################################################### ### code chunk number 5: SigFuge.Rnw:104-120 ################################################### calcInfo <- function(x) { info <- apply(x[["seq"]], 2, function(y) { y <- y[c("A", "C", "G", "T"), , drop=FALSE] cvg <- colSums(y) }) info } param <- ApplyPileupsParam(which=geneAnnot, what=c("seq", "qual"), yieldBy="position", yieldAll=TRUE) fls <- PileupFiles(bam_files, param=param) res <- applyPileups(fls, calcInfo, param=param) geneDepth <- t(do.call(cbind,res)) colnames(geneDepth) <- c("Sample1", "Sample2") geneDepth[500:505, ] ################################################### ### code chunk number 6: SigFuge.Rnw:131-133 ################################################### data(geneAnnot) data(geneDepth) ################################################### ### code chunk number 7: SigFuge.Rnw:139-144 ################################################### genename <- "CDKN2A" mdata <- geneDepth[,101:150] SFfigure(data=mdata, locusname=genename, annot=geneAnnot, lplots=1:3, savestr=genename, titlestr="CDKN2A locus, LUSC samples") ################################################### ### code chunk number 8: SigFuge.Rnw:179-183 (eval = FALSE) ################################################### ## genename <- "CDKN2A" ## SFfigure(data=mdata, locusname=genename, ## annot=geneAnnot, lplots=4, savestr ="CDKN2A_individual_curves", ## titlestr="CDKN2A locus, LUSC samples") ################################################### ### code chunk number 9: SigFuge.Rnw:189-191 ################################################### output <- SFpval(mdata) output@pvalnorm ################################################### ### code chunk number 10: SigFuge.Rnw:199-202 ################################################### norm <- SFnormalize(mdata) labels <- SFlabels(norm) labels[1:10] ################################################### ### code chunk number 11: SigFuge.Rnw:207-219 (eval = FALSE) ################################################### ## exp.data <- mdata[,-which(norm$flag == 1)] ## labels <- labels[-which(norm$flag == 1)] ## ## SFfigure(exp.data, locusname=genename, annot=geneAnnot, ## data.labels=labels, flag=0, lplots=2, ## savestr="Raw_CDKN2A", titlestr="Raw coverage", ## pval=0) ## ## SFfigure(norm$data.norm, locusname=genename, annot=geneAnnot, ## data.labels=labels, flag=0, lplots=2, ## savestr="Norm_CDKN2A", titlestr="Normalized coverage", ## pval=0)