inst/doc/groHMM.R

In groHMM: GRO-seq Analysis Pipeline

### R code from vignette source 'groHMM.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex(use.unsrturl=FALSE)
savewd = getwd()


###################################################
### code chunk number 2: install groHMM (eval = FALSE)
###################################################
## if (!requireNamespace("BiocManager", quietly=TRUE))
##     install.packages("BiocManager")
## BiocManager::install("groHMM")


###################################################
### code chunk number 3: install packages (eval = FALSE)
###################################################
## BiocManager::install("GenomicFeatures")
## BiocManager::install("org.Hs.eg.db")
## BiocManager::install("edgeR")
## BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")


###################################################
### code chunk number 4: use groHMM
###################################################
library(groHMM)
options(mc.cores=getCores(4))


###################################################
### code chunk number 5: load data
###################################################
S0mR1 <- as(readGAlignments(system.file("extdata", "S0mR1.bam",
        package="groHMM")), "GRanges")
S0mR2 <- as(readGAlignments(system.file("extdata", "S0mR1.bam",
        package="groHMM")), "GRanges") # Use R1 as R2
S40mR1 <- as(readGAlignments(system.file("extdata", "S40mR1.bam",
        package="groHMM")), "GRanges")
S40mR2 <- as(readGAlignments(system.file("extdata", "S40mR1.bam",
        package="groHMM")), "GRanges") # Use R1 as R2


###################################################
### code chunk number 6: combine replicates
###################################################
# Combine replicates
S0m <- c(S0mR1, S0mR2)
S40m <- c(S40mR1, S40mR2)


###################################################
### code chunk number 7: write wiggle files (eval = FALSE)
###################################################
## writeWiggle(reads=S0m, file="S0m_Plus.wig", fileType="wig", strand="+",
##         reverse=FALSE)
## writeWiggle(reads=S0m, file="S0m_Minus.wig", fileType="wig", strand="-",
##         reverse=TRUE)
## 
## # For BigWig file:
## # library(BSgenome.Hsapiens.UCSC.hg19)
## # si <- seqinfo(BSgenome.Hsapiens.UCSC.hg19)
## # writeWiggle(reads=S0m, file="S0m_Plus.wig", fileType="BigWig", strand="+",#
## #               reverse=FALSE, seqinfo=si)
## 
## # Normalized wiggle files
## expCounts <- mean(c(NROW(S0m), NROW(S40m)))
## writeWiggle(reads=S0m, file="S0m_Plus_Norm.wig", fileType="wig", strand="+",
##            normCounts=expCounts/NROW(S0m), reverse=FALSE)


###################################################
### code chunk number 8: transcript calling
###################################################
Sall <- sort(c(S0m, S40m))
# hmmResult <- detectTranscripts(Sall, LtProbB=-200, UTS=5,
#                   threshold=1)
# Load hmmResult from the saved previous run 
load(system.file("extdata", "Vignette.RData", package="groHMM"))
txHMM <- hmmResult$transcripts


###################################################
### code chunk number 9: called transcripts
###################################################
head(txHMM)


###################################################
### code chunk number 10: use ucsc genes
###################################################
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
kgdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

library(GenomicFeatures)
# For refseq annotations:
# rgdb <- makeTxDbFromUCSC(genome="hg19", tablename="refGene")
# saveDb(hg19RGdb, file="hg19RefGene.sqlite")
# rgdb <- loadDb("hg19RefGene.sqlite")
kgChr7 <- transcripts(kgdb, filter=list(tx_chrom = "chr7"),
                      columns=c("gene_id", "tx_id", "tx_name"))
seqlevels(kgChr7) <- seqlevelsInUse(kgChr7)


###################################################
### code chunk number 11: make consensus annotations
###################################################
# Collapse overlapping annotations
kgConsensus <- makeConsensusAnnotations(kgChr7, keytype="gene_id",
        mc.cores=getOption("mc.cores"))
library(org.Hs.eg.db)
map <- select(org.Hs.eg.db,
    keys=unlist(mcols(kgConsensus)$gene_id),
        columns=c("SYMBOL"), keytype=c("ENTREZID"))
mcols(kgConsensus)$symbol <- map$SYMBOL
mcols(kgConsensus)$type <- "gene"


###################################################
### code chunk number 12: evaluation of transcript calling
###################################################
e <- evaluateHMMInAnnotations(txHMM, kgConsensus)
e$eval


###################################################
### code chunk number 13: tuning HMM
###################################################
tune <- data.frame(LtProbB=c(rep(-100,3), rep(-200,3), rep(-300,3)),
           UTS=rep(c(5,10,15), 3))

Fp <- windowAnalysis(Sall, strand="+", windowSize=50)
Fm <- windowAnalysis(Sall, strand="-", windowSize=50)

# evals <- mclapply(seq_len(9), function(x) {
#        hmm <- detectTranscripts(Fp=Fp, Fm=Fm, LtProbB=tune$LtProbB[x],
#                  UTS=tune$UTS[x])
#        e <- evaluateHMMInAnnotations(hmm$transcripts, kgConsensus)
#        e$eval
#        }, mc.cores=getOption("mc.cores"),  mc.silent=TRUE)
tune <- cbind(tune, do.call(rbind, evals))  # evals from the previous run


###################################################
### code chunk number 14: check tuning
###################################################
tune
which.min(tune$total)
which.min(tune$errorRate)


###################################################
### code chunk number 15: expressed genes
###################################################
getExpressedAnnotations <- function(features, reads) {
     fLimit <- limitToXkb(features)
     count <- countOverlaps(fLimit, reads)
     features <- features[count!=0,]
     return(features[(quantile(width(features), .05) < width(features))
    & (width(features) < quantile(width(features), .95)),])}
conExpressed <- getExpressedAnnotations(features=kgConsensus,reads=Sall)


###################################################
### code chunk number 16: fig2plot
###################################################
td <- getTxDensity(txHMM, conExpressed, mc.cores=getOption("mc.cores"))
u <- par("usr")
lines(c(u[1], 0, 0, 1000, 1000, u[2]), c(0,0,u[4]-.04,u[4]-.04,0,0),
      col="red")
legend("topright", lty=c(1,1), col=c(2,1), c("ideal", "groHMM"))
text(c(-500,500), c(.05,.5), c("FivePrimeFP", "TP"))
td


###################################################
### code chunk number 17: fig2
###################################################
td <- getTxDensity(txHMM, conExpressed, mc.cores=getOption("mc.cores"))
u <- par("usr")
lines(c(u[1], 0, 0, 1000, 1000, u[2]), c(0,0,u[4]-.04,u[4]-.04,0,0),
      col="red")
legend("topright", lty=c(1,1), col=c(2,1), c("ideal", "groHMM"))
text(c(-500,500), c(.05,.5), c("FivePrimeFP", "TP"))
td


###################################################
### code chunk number 18: Non-mammalian (eval = FALSE)
###################################################
## # mysql --user=genome --host=genome-mysql.cse.ucsc.edu ce10 -e \
## #  "select chrom, txStart, txEnd, name, exonCount, strand, name2 from refGene \
## #   where chrom not like chrom!='chrM' and cdsStart != cdsEnd" | tail -n +1 > refgene.bed
## 
## # G <- read.table("refgene.bed", header=TRUE, stringsAsFactors=FALSE, sep="\t")
## # ce <- GRanges(G$chrom, IRanges(G$txStart, G$txEnd), strand=G$strand, \
## #               access=G$name, gene_id=G$name2)
## # ceConsensus <- makeConsensusAnnotations(ce, keytype="gene_id", \
## #                mc.cores=getOption("mc.cores"))


###################################################
### code chunk number 19: repairing called transcripts
###################################################
bPlus <- breakTranscriptsOnGenes(txHMM, kgConsensus, strand="+")
bMinus <- breakTranscriptsOnGenes(txHMM, kgConsensus, strand="-")
txBroken <- c(bPlus, bMinus)
txFinal <- combineTranscripts(txBroken, kgConsensus)


###################################################
### code chunk number 20: plotTxDensity_final (eval = FALSE)
###################################################
## tdFinal <- getTxDensity(txFinal, conExpressed, mc.cores=getOption("mc.cores"))


###################################################
### code chunk number 21: DE for trascripts
###################################################
# For called transcripts
library(edgeR)
txLimit <- limitToXkb(txFinal)
ctS0mR1 <- countOverlaps(txLimit, S0mR1)
ctS0mR2 <- countOverlaps(txLimit, S0mR2)
ctS40mR1 <- countOverlaps(txLimit, S40mR1)
ctS40mR2 <- countOverlaps(txLimit, S40mR2)

pcounts <- as.matrix(data.frame(ctS0mR1, ctS0mR2, ctS40mR1, ctS40mR2))
group <- factor(c("S0m", "S0m", "S40m", "S40m"))
lib.size <- c(NROW(S0mR1), NROW(S0mR2), NROW(S40mR1), NROW(S40mR2))
d <- DGEList(counts=pcounts, lib.size=lib.size, group=group)
d <- estimateCommonDisp(d)
et <- exactTest(d)
de <- decideTestsDGE(et, p=0.001, adjust="fdr")
detags <- seq_len(NROW(d))[as.logical(de)]
# Number of transcripts regulated at 40m
cat("up: ",sum(de==1), " down: ", sum(de==-1), "\n")


###################################################
### code chunk number 22: fig3plot
###################################################
plotSmear(et, de.tags=detags)
# 2 fold up or down
abline(h = c(-1,1), col="blue")


###################################################
### code chunk number 23: fig3
###################################################
plotSmear(et, de.tags=detags)
# 2 fold up or down
abline(h = c(-1,1), col="blue")


###################################################
### code chunk number 24: DE for annotated genes
###################################################
# For ucsc knownGenes
kgChr7 <- transcripts(kgdb, filter <- list(tx_chrom = "chr7"),
                      columns=c("gene_id", "tx_id", "tx_name"))
map <- select(org.Hs.eg.db,
        keys=unique(unlist(mcols(kgChr7)$gene_id)),
        columns=c("SYMBOL"), keytype=c("ENTREZID"))
missing <- elementNROWS(mcols(kgChr7)[,"gene_id"]) == 0
kgChr7 <- kgChr7[!missing,]

inx <- match(unlist(mcols(kgChr7)$gene_id), map$ENTREZID)
mcols(kgChr7)$symbol <- map[inx,"SYMBOL"]

kgLimit <- limitToXkb(kgChr7)
ctS0mR1 <- countOverlaps(kgLimit, S0mR1)
ctS0mR2 <- countOverlaps(kgLimit, S0mR2)
ctS40mR1 <- countOverlaps(kgLimit, S40mR1)
ctS40mR2 <- countOverlaps(kgLimit, S40mR2)

counts <- as.matrix(data.frame(ctS0mR1, ctS0mR2, ctS40mR1, ctS40mR2))
group <- factor(c("S0m", "S0m", "S40m", "S40m"))
lib.size <- c(NROW(S0mR1), NROW(S0mR2), NROW(S40mR1), NROW(S40mR2))
d <- DGEList(counts=counts, lib.size=lib.size, group=group)

d <- estimateCommonDisp(d)
et <- exactTest(d)
de <- decideTestsDGE(et, p=0.001, adjust="fdr")
detags <- seq_len(NROW(d))[as.logical(de)]
symbols <- mcols(kgChr7)$symbol
# Number of unique genes regulated at 40m
cat("up: ", NROW(unique(symbols[de==1])), "\n")
cat("down: ", NROW(unique(symbols[de==-1])), "\n")
plotSmear(et, de.tags=detags)
abline(h = c(-1,1), col="blue")


###################################################
### code chunk number 25: metagene analysis
###################################################
upGenes <- kgChr7[de==1,]
expReads <- mean(c(NROW(S0m), NROW(S40m)))
# Metagene around TSS
mg0m <- runMetaGene(features=upGenes, reads=S0m, size=100,
                       normCounts=expReads/NROW(S0m), sampling=FALSE,
               mc.cores=getOption("mc.cores"))
mg40m <- runMetaGene(features=upGenes, reads=S40m, size=100,
                        normCounts=expReads/NROW(S40m), sampling=FALSE,
                mc.cores=getOption("mc.cores"))


###################################################
### code chunk number 26: plot metagene
###################################################
plotMetaGene <- function(POS=c(-10000:+9999), mg, MIN, MAX){
    plot(POS, mg$sense, col="red", type="h", xlim=c(-5000, 5000),
        ylim=c(floor(MIN),ceiling(MAX)), ylab="Read Density",
        xlab="Position (relative to TSS)")
     points(POS, (-1*rev(mg$antisense)), col="blue", type="h")
     abline(mean(mg$sense[5000:8000]), 0, lty="dotted")
}

MAX <- max(c(mg0m$sense, mg40m$sense))
MIN <- -1*max(c(mg0m$antisense, mg40m$antisense))
plotMetaGene(mg=mg0m, MIN=MIN, MAX=MAX)
plotMetaGene(mg=mg40m, MIN=MIN, MAX=MAX)


###################################################
### code chunk number 27: metagene figures
###################################################
png(filename="metagene0m.png")
plotMetaGene(mg=mg0m, MIN=MIN, MAX=MAX)
dev.off()
png(filename="metagene40m.png", width=480, height=480)
plotMetaGene(mg=mg40m, MIN=MIN, MAX=MAX)
dev.off()


###################################################
### code chunk number 28: save data for vignette (eval = FALSE)
###################################################
## save(hmmResult, evals, file="Vignette.RData")


###################################################
### code chunk number 29: sessionInfo
###################################################
toLatex(sessionInfo())


###################################################
### code chunk number 30: groHMM.Rnw:599-600
###################################################
setwd(savewd)