1. 读入文件

library(DiffBind)
dbObj <- dba(sampleSheet="SampleSheet.csv")

2. 计算亲和结合矩阵

dbObj <- dba.count(dbObj, bUseSummarizeOverlaps=TRUE)
dba.plotPCA(dbObj,  attributes=DBA_FACTOR, label=DBA_ID)
plot(dbObj)

3.计算差异分析

# Establishing a contrast 
dbObj <- dba.contrast(dbObj, categories=DBA_FACTOR,minMembers = 2)
dbObj <- dba.analyze(dbObj, method=DBA_ALL_METHODS)
#  summary of results
dba.show(dbObj, bContrasts=T)
#  overlapping peaks identified by the two different tools (DESeq2 and edgeR)
dba.plotVenn(dbObj,contrast=1,method=DBA_ALL_METHODS)

4. 提取结果

comp1.edgeR <- dba.report(dbObj, method=DBA_EDGER, contrast = 1, th=1)
comp1.deseq <- dba.report(dbObj, method=DBA_DESEQ2, contrast = 1, th=1)

5-1. 保存文件--DEseq2 结果

# DESeq2--所有差异
out <- as.data.frame(comp1.deseq)
write.table(out, file="ACL_vs_WT_deseq2.txt", sep="\t", quote=F, col.names = NA)

ALL.bed <- out[ c("seqnames", "start", "end", "strand", "Fold")]

write.table(ALL.bed, file="ACL_vs_WT_ALL.bed", sep="\t", quote=F, row.names=F, col.names=F)

# 以bed格式保存显著性的差异结果

out <- as.data.frame(comp1.deseq)
#Down_regualted
deseq_down.bed <- out[ which(out$p.value < 0.05 & out$Fold <= -0.58), c("seqnames", "start", "end", "strand", "Fold")]
                 
write.table(deseq_down.bed, file="ACL_vs_WT_Down_deseq2.bed", sep="\t", quote=F, row.names=F, col.names=F)
#Up_regualted    
deseq_up.bed <- out[ which(out$p.value < 0.05 & out$Fold >= 0.58), c("seqnames", "start", "end", "strand", "Fold")]
                                                 
write.table(deseq_up.bed, file="ACL_vs_WT_Up.bed", sep="\t", quote=F, row.names=F, col.names=F)

5-2 保存文件--EdgeR 结果

# EdgeR---所有差异
out <- as.data.frame(comp1.edgeR)
write.table(out, file="HAG_vs_MIX_edgeR.txt", sep="\t", quote=F, col.names = NA)

# 以bed格式保存显著性的差异结果
# Create bed files for each keeping only significant peaks (p < 0.05)
#Down_regulated

out <- as.data.frame(comp1.edgeR)
edge_down.bed <- out[ which(out$p.value <= 0.05 & out$Fold <= -0.58), c("seqnames", "start", "end", "strand", "Fold")]
write.table(edge_down.bed, file="HAG_vs_WT_Down_edgeR_sig.bed", sep="\t", quote=F, row.names=F, col.names=F)

#Up_regulated
edge_up.bed <- out[ which(out$p.value <= 0.05 & out$Fold >= 0.58), c("seqnames", "start", "end", "strand", "Fold")]
write.table(edge_up.bed, file="HAG_vs_WT_Up_edgeR_sig.bed", sep="\t", quote=F, row.names=F, col.names=F)