inst/doc/SigFuge.R

In SigFuge: SigFuge

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