diffbind

samples <- read.csv(file.path(system.file(“extra”, package=“DiffBind”),
“tamoxifen.csv”))
names(samples)
samples

setwd(system.file(‘extra’,package=‘DiffBind’))

tamoxifen=dba(sampleSheet=“tamoxifen.csv”,dir=system.file(“extra”, package=“DiffBind”))

tamoxifen

plot(tamoxifen)

data(tamoxifen_peaks)
tamoxifen

看看现在的peaks数量

看看现在的peaks数量

peakdata <- dba.show(tamoxifen)$Intervals

去除以后

tamoxifen <- dba.blacklist(tamoxifen, blacklist=DBA_BLACKLIST_HG19,
greylist=FALSE)

去除以后的数量

peakdata.BL <- dba.show(tamoxifen)$Intervals

看看每个样本去除了多少的blacklist区域peaks

peakdata - peakdata.BL

[1] 0 0 1 1 1 0 0 0 0 0 0

可以看到，这里只有中间的三个样本（MCF7）去掉了1的peaks

tamoxifen <- dba.count(tamoxifen)

data(tamoxifen_counts)
#load(“C:/Users/Administrator/Documents/R/win-library/4.0/DiffBind/data/tamoxifen_counts.rda”)
tamoxifen_counts

tamoxifen <- dba.normalize(tamoxifen)
norm <- dba.normalize(tamoxifen, bRetrieve=TRUE)
norm

info <- dba.show(tamoxifen)
libsizes <- cbind(LibReads=info R e a d s , F R i P = i n f o Reads, FRiP=info Reads,FRiP=infoFRiP,
PeakReads=round(info R e a d s ∗ i n f o Reads * info Reads∗infoFRiP))
rownames(libsizes) <- info$ID
libsizes

这里默认也是基于测序深度进行normalize

tamoxifen <- dba.normalize(tamoxifen)

可以看一下具体的细节

norm <- dba.normalize(tamoxifen, bRetrieve=TRUE)
norm
normlibs <- cbind(FullLibSize=norm l i b . s i z e s , N o r m F a c s = n o r m lib.sizes, NormFacs=norm lib.sizes,NormFacs=normnorm.factors,NormLibSize=round(norm l i b . s i z e s / n o r m lib.sizes/norm lib.sizes/normnorm.factors))

rownames(normlibs) <- info$ID

tamoxifen <- dba.contrast(tamoxifen,
reorderMeta=list(Condition=“Responsive”))

tamoxifen <- dba.analyze(tamoxifen)

#> dbObj <- dba.contrast(dbObj,

+ reorderMeta=list(Condition=“NC”))

#Computing results names…
#> dbObj <- dba.analyze(dbObj)
#analyze失败采用下面的

#dbObj <- dba.analyze(dbObj, method=DBA_ALL_METHODS,bBlacklist=FALSE,bGreylist=FALSE)
#If no normalization has been specified, the reads will be normalized via a call to dba.normalize using defaults.

dba.show(tamoxifen, bContrasts=TRUE)
tamoxifen.DB <- dba.report(tamoxifen)

sum(tamoxifen.DB$Fold>0)

sum(tamoxifen.DB$Fold<0)

dba.plotVenn(tamoxifen, contrast=1, bDB=TRUE,
bGain=TRUE, bLoss=TRUE, bAll=FALSE)

如果使用全部的2845 sites

dba.plotPCA(tamoxifen,DBA_TISSUE,label=DBA_CONDITION)

如果使用差异的246sites

dba.plotPCA(tamoxifen, contrast=1, label=DBA_TISSUE)

dba.plotVolcano(tamoxifen)

使用全部的sites

corvals <- dba.plotHeatmap(tamoxifen)

默认使用DBA_SCORE_NORMALIZED，还可以用另一种计算normalization score的方法：RPKM fold (RPKM of the ChIP reads divided by RPKM of the control reads)

dba.plotHeatmap(tamoxifen, score=DBA_SCORE_RPKM_FOLD)

如果只用差异的246sites

hmap <- colorRampPalette(c(“red”, “black”, “green”))(n = 13)

readscores <- dba.plotHeatmap(tamoxifen, contrast=1, correlations=FALSE,scale=“row”, colScheme = hmap)

setwd(“C:/Users/Administrator/Documents/R/win-library/4.0/DiffBind/data”)
data(tamoxifen_peaks)
tamoxifen_peaks <- load(“C:/Users/Administrator/Documents/R/win-library/4.0/DiffBind/data/tamoxifen_peaks.rda”)
olap.rate <- dba.overlap(tamoxifen,mode=DBA_OLAP_RATE)

plot(olap.rate,type = ‘b’,ylab=’# peaks’,xlab = ‘Overlap at least this many peaksets’)

names(tamoxifen$masks)

dba.overlap(tamoxifen,tamoxifen m a s k s masks masksMCF7 & tamoxifen m a s k s masks masksResponsive,
mode=DBA_OLAP_RATE)

#dba.overlap(dbObj,dbObj m a s k s masks masksNC & dbObj m a s k s masks masksHD,

mode=DBA_OLAP_RATE)

dba.plotVenn(tamoxifen, tamoxifen m a s k s masks masksMCF7 & tamoxifen m a s k s masks masksResponsive)

tamoxifen_consensus <- dba.peakset(tamoxifen,
consensus=c(DBA_TISSUE,DBA_CONDITION),
minOverlap=0.66)

tamoxifen_consensus <- dba(tamoxifen_consensus,
mask=tamoxifen_consensus m a s k s masks masksConsensus,
minOverlap=1)

consensus_peaks <- dba.peakset(tamoxifen_consensus, bRetrieve=TRUE)

tamoxifen <- dba.count(tamoxifen, peaks=consensus_peaks)

data(tamoxifen_peaks)
tamoxifen <- dba.peakset(tamoxifen, consensus=DBA_TISSUE, minOverlap=0.66)

cons.ol <- dba.plotVenn(tamoxifen, tamoxifen m a s k s masks masksConsensus)

data(tamoxifen_peaks)
dba.overlap(tamoxifen,tamoxifen m a s k s masks masksResistant,mode=DBA_OLAP_RATE)
dba.overlap(tamoxifen,tamoxifen m a s k s masks masksResponsive,mode=DBA_OLAP_RATE)

tamoxifen <- dba.peakset(tamoxifen, consensus=DBA_CONDITION, minOverlap=0.33)
dba.plotVenn(tamoxifen,tamoxifen m a s k s masks masksConsensus)

tamoxifen.OL <- dba.overlap(tamoxifen, tamoxifen m a s k s masks masksConsensus)

tamoxifen.OL o n l y A t a m o x i f e n . O L onlyA tamoxifen.OL onlyAtamoxifen.OLonlyB

tamoxifen <- dba.peakset(tamoxifen,tamoxifen m a s k s masks masksConsensus,
minOverlap=1,sampID=“OL Consensus”)

tamoxifen <- dba.peakset(tamoxifen,!tamoxifen m a s k s masks masksConsensus,
minOverlap=3,sampID=“Consensus_3”)
dba.plotVenn(tamoxifen,14:15)

data(tamoxifen_analysis)
tamoxifen.rep <- dba.report(tamoxifen,bCalled=TRUE,th=1)
onlyResistant <- tamoxifen.repKaTeX parse error: Expected 'EOF', got '&' at position 12: Called1>=2 &̲ tamoxifen.repCalled2<3
sum(onlyResistant )
onlyResponsive <- tamoxifen.repKaTeX parse error: Expected 'EOF', got '&' at position 12: Called2>=3 &̲ tamoxifen.repCalled1<2
sum(onlyResponsive)

bothGroups <- tamoxifen.repKaTeX parse error: Expected 'EOF', got '&' at position 13: Called1>= 2 &̲ tamoxifen.repCalled2>=3
sum(bothGroups)

tamoxifen.DB <- dba.report(tamoxifen,bCalled=TRUE)
onlyResistant.DB <- (tamoxifen.DBKaTeX parse error: Expected 'EOF', got '&' at position 15: Called1 >= 2) &̲ (tamoxifen.DBCalled2 < 3)
sum(onlyResistant.DB)

onlyResponsive.DB <- (tamoxifen.DBKaTeX parse error: Expected 'EOF', got '&' at position 15: Called2 >= 3) &̲ (tamoxifen.DBCalled1 < 2)
sum(onlyResponsive.DB)

bothGroups.DB <- (tamoxifen.DBKaTeX parse error: Expected 'EOF', got '&' at position 15: Called1 >= 2) &̲ (tamoxifen.DBCalled2 >= 3)
sum(bothGroups.DB)

neitherGroup.DB <- (tamoxifen.DBKaTeX parse error: Expected 'EOF', got '&' at position 14: Called1 < 2) &̲ (tamoxifen.DBCalled2 < 3)
sum(neitherGroup.DB)

标签：peaks,diffbind,DBA,tamoxifen,DB,masks,dba
来源： https://blog.csdn.net/qq_39306047/article/details/109990452