Entering edit mode
21 months ago
diqixiaoyaoer
▴
20
Here is my code:
####
library(bsseq)
library(bsseqData)
## ----showData-----------------------------------------------------------------
data(BS.cancer.ex)
BS.cancer.ex <- updateObject(BS.cancer.ex)
BS.cancer.ex
pData(BS.cancer.ex)
## ----smooth,eval=FALSE--------------------------------------------------------
# BS.cancer.ex.fit <- BSmooth(
# BSseq = BS.cancer.ex,
# BPPARAM = MulticoreParam(workers = 1),
# verbose = TRUE)
## ----showDataFit--------------------------------------------------------------
data(BS.cancer.ex.fit)
BS.cancer.ex.fit <- updateObject(BS.cancer.ex.fit)
BS.cancer.ex.fit
## ----cpgNumbers---------------------------------------------------------------
## The average coverage of CpGs on the two chromosomes
round(colMeans(getCoverage(BS.cancer.ex)), 1)
## Number of CpGs in two chromosomes
length(BS.cancer.ex)
## Number of CpGs which are covered by at least 1 read in all 6 samples
sum(rowSums(getCoverage(BS.cancer.ex) >= 1) == 6)
## Number of CpGs with 0 coverage in all samples
sum(rowSums(getCoverage(BS.cancer.ex)) == 0)
## ----poisson------------------------------------------------------------------
logp <- ppois(0, lambda = 4, lower.tail = FALSE, log.p = TRUE)
round(1 - exp(6 * logp), 3)
## ----smoothSplit,eval=FALSE---------------------------------------------------
# ## Split datag
# BS1 <- BS.cancer.ex[, 1]
# save(BS1, file = "BS1.rda")
# BS2 <- BS.cancer.ex[, 2]
# save(BS1, file = "BS1.rda")
# ## done splitting
#
# ## Do the following on each node
#
# ## node 1
# load("BS1.rda")
# BS1.fit <- BSmooth(BS1)
# save(BS1.fit)
# save(BS1.fit, file = "BS1.fit.rda")
# ## done node 1
#
# ## node 2
# load("BS2.rda")
# BS2.fit <- BSmooth(BS2)
# save(BS2.fit, file = "BS2.fit.rda")
# ## done node 2
#
# ## join; in a new R session
# load("BS1.fit.rda")
# load("BS2.fit.rda")
# BS.fit <- combine(BS1.fit, BS2.fit)
## ----keepLoci-----------------------------------------------------------------
BS.cov <- getCoverage(BS.cancer.ex.fit)
keepLoci.ex <- which(rowSums(BS.cov[, BS.cancer.ex$Type == "cancer"] >= 2) >= 2 &
rowSums(BS.cov[, BS.cancer.ex$Type == "normal"] >= 2) >= 2)
length(keepLoci.ex)
BS.cancer.ex.fit <- BS.cancer.ex.fit[keepLoci.ex,]
## ----BSmooth.tstat------------------------------------------------------------
BS.cancer.ex.tstat <- BSmooth.tstat(BS.cancer.ex.fit,
group1 = c("C1", "C2", "C3"),
group2 = c("N1", "N2", "N3"),
estimate.var = "group2",
local.correct = TRUE,
verbose = TRUE)
BS.cancer.ex.tstat
## ----plotTstat,fig.width=4,fig.height=4---------------------------------------
plot(BS.cancer.ex.tstat)
## ----dmrs---------------------------------------------------------------------
dmrs0 <- dmrFinder(BS.cancer.ex.tstat, cutoff = c(-4.6, 4.6))
dmrs <- subset(dmrs0, n >= 3 & abs(meanDiff) >= 0.1)
nrow(dmrs)
head(dmrs, n = 3)
## ----plotSetup----------------------------------------------------------------
pData <- pData(BS.cancer.ex.fit)
pData$col <- rep(c("red", "blue"), each = 3)
pData(BS.cancer.ex.fit) <- pData
# rm(plotRegion)
## ----plotDmr,fig.width=8,fig.height=4-----------------------------------------
# ?plotRegion
plotRegion(BS.cancer.ex.fit, dmrs[1,], extend = 5000, addRegions = dmrs)
plotRegion(BS.cancer.ex.fit, dmrs[1,], extend = 5000)
# names(annot.chr21)[2]
# geneTrack<-data.frame(annot.chr21[2])
plotRegion(BS.cancer.ex.fit, dmrs[1,], extend = 5000,addRegions = dmrs,addTicks = T,annoTrack=annot.chr21)
# source("H:\\R_bsseq\\bsseq-master\\R\\plotting.R")
head(geneTrack)
class(geneTrack)
geneTrack<-data.frame(annot.chr21[2])
unique(geneTrack$group)
names(geneTrack)<-c("exon_number","group_name","chr","start","end","width","strand","id","tx_id","gene_ID","gene_name","type")
# ?plotRegion
plotRegion(BS.cancer.ex.fit, dmrs[1,], extend = 5000,addRegions = dmrs,
addTicks = T,
# col = "black",
# annoTrack=annot.chr21[2],
#annoTrack=annoTrack,
stat = "tstat.corrected",
stat.col = "black", stat.lwd = 1,
# regionCol = "blue",
cex.anno = 0.8,
# byt = "o",
#addPoints = T, pointsMinCov = 5, highlightMain = FALSE,
geneTrack = geneTrack,
cex.gene = 1.5,
lwd = 4)
I have a small quesiton about the bsseq offical code in last part:
plotRegion(BS.cancer.ex.fit, dmrs[1,], extend = 5000,addRegions = dmrs, addTicks = T,
# col = "black", # annoTrack=annot.chr21[2], #annoTrack=annoTrack, stat = "tstat.corrected", stat.col = "black", stat.lwd = 1, # regionCol = "blue", cex.anno = 0.8, # byt = "o", #addPoints = T, pointsMinCov = 5, highlightMain = FALSE, geneTrack = geneTrack, cex.gene = 1.5, lwd = 4)
I want to add gene information, including Exon, promoter , in the figure.
So here I used annot.chr21[2] as the geneTrack.
However, it appeared with wrong information: Error in if (direction == "+") { : argument is of length zero
I don't know why.
And I also checked the bsseq plotting source code: https://github.com/hansenlab/bsseq/blob/master/R/plotting.R
plotAnnoTrack <- function(gr, annoTrack, cex) {
## check may need to be modified
if(!all(sapply(annoTrack, function(xx) is(xx, "GRanges"))))
stop("all elements in 'annoTrack' needs to be 'GRanges'")
plot(start(gr), 1, type = "n", xaxt = "n", yaxt = "n", bty = "n",
ylim = c(0.5, length(annoTrack) + 0.5), xlim = c(start(gr), end(gr)), xlab = "", ylab = "")
lapply(seq(along = annoTrack), function(ii) {
jj <- length(annoTrack) + 1- ii
ir <- subsetByOverlaps(annoTrack[[ii]], gr)
if(length(ir) > 0)
rect(c(start(ir)-0.5,start(ir)-0.5),
c(jj+0.33 ,jj+0.33 ),
c(end(ir),end(ir)),
c(jj +0.43,jj+0.53 ),
col = c("black","#148bda"),
border = NA)
#mtext(names(annoTrack)[1], side = 2, at = jj, las = 1, line = 1,
# cex = cex, adj = 0.9,
# col = "red")
})
}
# Based on plotAnnoTrack() and
# https://jhu-genomics.slack.com/archives/hansen_gtex/p1461760583000142
# TODO: Use standard Bioconductor object for storing gene information rather
# than this ad hoc data.frame (e.g., based on TxDb.* packages)
plotGeneTrack <- function(gr, geneTrack, cex) {
geneTrack_gr <- makeGRangesFromDataFrame(geneTrack)
ol <- findOverlaps(geneTrack_gr, gr)
genes <- geneTrack[queryHits(ol), ]
plot(start(gr), 1, type = "n", xaxt = "n", yaxt = "n", bty = "n",
ylim = c(-1.5, 1.5), xlim = c(start(gr), end(gr)),
xlab = "", ylab = "", cex.lab = 4, lheight = 2, cex.axis = 1)
if (nrow(genes) > 0) {
for (g in 1:nrow(genes)) {
geneind2 = which(geneTrack$gene_name == genes$gene_name[g])
geneind2 = geneind2[which(geneTrack$isoforms[geneind2] == 1)]
direction = unique(geneTrack$strand[geneind2])
ES = geneTrack$start[geneind2]
EE = geneTrack$end[geneind2]
Exons = cbind(ES, EE)
if (direction == "+") {
lines(x = c(min(ES), max(EE)),
y = c(0.65, 0.65))
apply(Exons, 1, function(x) {
polygon(c(x[1], x[2], x[2], x[1]),
c(0.45, 0.45, 0.85, 0.85),
col = "darkgrey")
})
text((max(start(gr), min(ES)) +
min(end(gr), max(EE))) / 2, 1.2,
genes$gene_name[g], cex = cex)
} else {
lines(x = c(min(ES), max(EE)),
y = c(-0.65, -0.65))
apply(Exons, 1, function(x)
polygon(c(x[1], x[2], x[2], x[1]),
c(-0.45, -0.45, -0.85, -0.85),
col = "darkgrey"))
text((max(start(gr), min(ES)
) + min(end(gr), max(EE)
)) / 2, -1.2, genes$gene_name[g], cex = cex)
}
}
}
}
.......
But I still don't know how to solve this problem.
Error in if (direction == "+") { : argument is of length zero
I hope someone can give me some advice. I will be very thankful.
