ChIP-seq | ATAC-seq | Cut&Run | Cut&Tag | DiffBind使用教程 | dba.plotProfile

2024年09月09日

Diffbind其实功能非常强大，我还在探索。可以极速分析成图，也可以精细的个性化定制。工程化做得非常好，代码都被封装起来了，看不到源码。

dba.plotProfile，这个核心的Heatmap成图就很厉害。

今天，我花了快一个多小时来探索，如何自定义Heatmap的颜色，正常的需求就是按列的样本来设定不同的color

我才发现dba.plotProfile函数被用了两次，第一次是生成数据，第二次是绘图。

核心的参数有三个：

sites=tmp.gr，这个决定了row的peak list，如果就是一个Grange，那就是一个颜色，如果是Grangelist，那就是自动分组；
samples = sample.list，这个决定了你用哪些样本，可以是简单的vector，但如果想指定颜色，那就必须是个list；
matrices_color = color.list，这个决定了最核心的Heatmap的color，这里用的circlize::colorRamp2，如果直接用colorRampPalette就会有scale的问题

代码：

tmp <- dmso.tf.dba.obj$samples$SampleID
# for subset samples
tmp <- 1:length(dmso.tf.dba.obj$samples$SampleID)
names(tmp) <- dmso.tf.dba.obj$samples$SampleID
tmp

tmp.gr <- GRanges(norm.count[,1:3])

tmp.order <- c('DMSO_SOX9_2_Paul',#'DMSO_SMARCC1_2_Paul',
'dT29_SC1_DM1_Rag_b4','dT29_ARD1_DM2_Rag_b4',
'dT29_BRD9_DM1_Rag_b4','dT29_PBR1_DM1_Rag_b4')

library(RColorBrewer)
tmp.colors <- brewer.pal(5,"Set1")
names(tmp.colors) <- names(tmp[tmp.order])

sample.list <- list()
color.list <- list()
for (i in names(tmp[tmp.order])) {
    print(i)
    sample.list[[i]] <- tmp[i]
    # color.list[[i]] <- colorRampPalette(c("white", tmp.colors[i]))(10)
    color.list[[i]] <- circlize::colorRamp2(c(0, 6), c("white", tmp.colors[i]))
}

profiles <- dba.plotProfile(dmso.tf.dba.obj, sites=tmp.gr, merge=NULL, samples = sample.list)

options(repr.plot.width=10, repr.plot.height=7)
dba.plotProfile(profiles, use_raster=F, matrices_color = color.list, all_color_scales = F)

参数文档在这里：?profileplyr::generateEnrichedHeatmap

代码地址：http://localhost:17449/lab/tree/projects/BAF_SOX9/diffbind/6.DMSO_only.ipynb

这个需要一点R功底，解决之后真的非常有成就感！！！

核心可视化：

https://content.cruk.cam.ac.uk/bioinformatics/software/DiffBind/plotProfileDemo.html

2023年02月20日

没想到已经两个月整了，感觉项目丝毫没有任何进展，不知道自己每天都在干什么！

目前我对表观的数据已经有了比较深刻的认识了，cut&run的一个核心就是对IgG和spike-in进行校正，这是分析环节的最大门槛。

如果不对IgG进行校正，只任其跟随lib-size校正，那背景噪音将会非常大，尤其是在bigwig里可视化，所以必须按spike-in进行down-sampling。

但spike-in又不是万能的，正常样本里的差异也是很大的，甚至是成倍的，如果无视其直接校正，那一切都是mess，所以对正常样本（IgG）我们可以忽视spike-in。

我个人认为IgG在分析时基本可以忽视，按spike-in得到一个normalized bigwig文件即可，结局几乎是一定的。

RNA-seq的QC比较简单，所以如果Cut&Run有对应的RNA-seq，那就直接看DEG周围的peak是否有差异。

结果是必然的！

Diffbind高级分析法：notebooks/projects/cut_run/HDAC-b2/pipeline/Diffbind.ipynb

DBA_NORM_LIB ("lib") Normalize by library size only. Library sizes can be specified using the library parameter. Normalization factors will be calculated to give each equal weight in a manner appropriate for the analysis method. See also the libFun parameter, which can be used to scale the normalization factors for DESeq2.

这篇教程不错：DiffBind使用探索

2022年12月14日

现在基本掌握了所有技术。

normalization的逻辑和实操，IgG和experiment分开norm，组内用library size。
目前我们只能考虑proximal的peak，distal基本可以忽略，除非有Hi-c的数据。
差异分析也基本搞通，初学者不建议使用DiffBind，用好deeptools足矣。

2022年12月09日

经过一周多的学习训练，已经掌握了数据分析的大致流程，但有一些核心的东西还没搞懂。

如何对不同样本的数据做normalization，因为包含了control，复杂了一些。因为control本身起始DNA就很低，如何按照library size直接做normalization，那必然就会有噪音，这在IGV里可视化尤为明显；
H3K27Ac的Peak应该在center有一个drop，这才是可靠的Peak，有点类似footprint。
抗体的特异性不够，DNA量不够怎么办？比对到hg38的rate太低怎么办？

Alignment frequencies are expected to be >80% for high-quality data.

CUT&Tag data typically has very low backgrounds, so as few as 1 million mapped fragments can give robust profiles for a histone modification in the human genome.

Profiling of less-abundant transcription factors and chromatin proteins may require 10 times as many mapped fragments for downstream analysis.

低于80%的mapping rate就说明抗体特异性有问题；
对于histone modification，1 M的mapped fragments即可；
而对于丰度没那么多的transcription factors，则需要10倍，即10 M的mapped fragments；

参考：

ChIP-Seq数据中包含了spike-in怎么分析

2022年12月03日

基本分析流程都挺容易得，就是设计层面有些东西需要注意。

样本的数量
送样的细胞数量；
下机的测序深度；

根据目前经验，样本数和下机的测序深度对结果影响非常大【差异peak的数量】，一定需要选择和校正。

可以对原始的bam文件进行downsampling，也需要选取合适的样本进行差异分析，在PCA里选择合适的样本。

结果一定要符合常识。

核心的可视化：

pdf("plot.pdf", width=8, height=8)
plot(exp1_2)
dev.off()

pdf("plotPCA.pdf", width=8, height=8)
dba.plotPCA(exp1_2,  attributes=DBA_CONDITION, label=DBA_REPLICATE, labelSize=0.4)
dev.off()

pdf("plotHeatmap", width=10, height=10)
dba.plotHeatmap(exp1_2)
dev.off()

pdf("plotHeatmap", width=10, height=10)
dba.plotProfile(exp1_2)
dev.off()

深度可视化

TSS heatmap - 只要有bam和bed就可以画图，非常高大上。

dba.plotProfile Demonstration

pdf("plot.pdf", width=8, height=8)
plot(exp1_2)
dev.off()

pdf("plotPCA.pdf", width=8, height=8)
dba.plotPCA(exp1_2,  attributes=DBA_CONDITION, label=DBA_REPLICATE, labelSize=0.4)
dev.off()

pdf("plotHeatmap.pdf", width=10, height=10)
hmap <- colorRampPalette(c("red", "black", "green"))(n = 13)
dba.plotHeatmap(exp1_2, correlations=FALSE,scale="row", colScheme = hmap)
dev.off()

N3_SOX9_E1 <- dba.contrast(exp1_2, contrast=c("Condition","N3_SOX9_EXP1","N3_IgG_EXP1"))

exp1_2$config$RunParallel <- TRUE
profiles <- dba.plotProfile(exp1_2)

exp1_2$masks$AD_SOX9 <- exp1_2$masks$AD3_SOX9_EXP1
exp1_2$masks$AD_SOX9[exp1_2$masks$AD_SOX9] <- F
exp1_2$masks$AD_SOX9[c("AD3_SOX9_EXP1_1","AD3_SOX9_EXP1_3","AD3_SOX9_EXP2_1","AD3_SOX9_EXP2_2")] <- T

exp1_2$masks$N_SOX9 <- exp1_2$masks$AD3_SOX9_EXP1
exp1_2$masks$N_SOX9[exp1_2$masks$N_SOX9] <- F
exp1_2$masks$N_SOX9[c("N3_SOX9_EXP1_1","N3_SOX9_EXP1_3","N3_SOX9_EXP2_1","N3_SOX9_EXP2_2")] <- T


AbSC <- read.table("sox9_paper_bed/AbSC3.bed")
ISC <- read.table("sox9_paper_bed/ISC3.bed")
Mustata <- read.table("sox9_paper_bed/Mustata3.bed")
Nusse <- read.table("sox9_paper_bed/Nusse3.bed")
Vallone <- read.table("sox9_paper_bed/Vallone3.bed")

repList <- GRangesList(AbSC=paste(paste(AbSC3$V1,AbSC3$V2,sep=":"), AbSC3$V3, sep="-"),
						ISC=paste(paste(ISC$V1,ISC$V2,sep=":"), ISC$V3, sep="-"),
						Mustata=paste(paste(Mustata$V1,Mustata$V2,sep=":"), Mustata$V3, sep="-"),
						Nusse=paste(paste(Nusse$V1,Nusse$V2,sep=":"), Nusse$V3, sep="-"),
						Vallone=paste(paste(Vallone$V1,Vallone$V2,sep=":"), Vallone$V3, sep="-")
						)
profiles <- dba.plotProfile(exp1_2, sites=repList, score="Fold",
							samples=list(Normal=exp1_2$masks$N_SOX9,
                                         Adenoma=exp1_2$masks$AD_SOX9)
                                         )

pdf("plotProfile.pdf", width=10, height=15)
dba.plotProfile(profiles, use_raster=F)
dev.off()

挺好用，只需要bam和peak就能做差异分析。

准备工作，一个meta info的sample file。主要是bam和bed的位置，bed可以用macs一行命令合并所有bam来做peak calling。

library(DiffBind)
# library(tidyverse)
samples <- read.csv("EXP1.csv")
exp1 <- dba(sampleSheet=samples)
exp1 <- dba.count(exp1) 

raw.count <- dba.peakset(exp1, bRetrieve=TRUE, DataType = DBA_DATA_FRAME)
rownames(raw.count) <- paste(raw.count$CHR, raw.count$START, raw.count$END, sep = ":")

exp1 <- dba.normalize(exp1, normalize=DBA_NORM_TMM)

norm.count <- dba.peakset(exp1, bRetrieve=TRUE, DataType = DBA_DATA_FRAME)
rownames(norm.count) <- paste(norm.count$CHR, norm.count$START, norm.count$END, sep = ":")

# first is case, second is control
AD3_SOX9 <- dba.contrast(exp1, contrast=c("Condition","AD3_SOX9_EXP1","AD3_IgG_EXP1"))
AD3_SOX9 <- dba.analyze(AD3_SOX9)
AD3_SOX9.report <- dba.report(AD3_SOX9, DataType=DBA_DATA_FRAME)
rownames(AD3_SOX9.report) <- paste(AD3_SOX9.report$Chr, AD3_SOX9.report$Start, AD3_SOX9.report$End, sep = ":")

N3_SOX9 <- dba.contrast(exp1, contrast=c("Condition","N3_SOX9_EXP1","N3_IgG_EXP1"))
N3_SOX9 <- dba.analyze(N3_SOX9)
N3_SOX9.report <- dba.report(N3_SOX9, DataType=DBA_DATA_FRAME)
rownames(N3_SOX9.report) <- paste(N3_SOX9.report$Chr, N3_SOX9.report$Start, N3_SOX9.report$End, sep = ":")

# save
save(samples, raw.count, norm.count, AD3_SOX9, AD3_SOX9.report, N3_SOX9, N3_SOX9.report, file = "DiffBind.result.Rdata")

核心就是设置好condition，后期可以指定任意两个condition来比较。

其次就是normalization，这里用的edgeR的TMM方法，注意check是否normalization成功。

最后就是手动check 软件生成的log2FC方向是否正确。

保存所有中间文件，以备后续分析。

posted @ 2022-12-03 07:05 Life·Intelligence 阅读(2199) 评论(0) 编辑收藏举报

刷新页面返回顶部

Digital-LI

ChIP-seq | ATAC-seq | Cut&Run | Cut&Tag | DiffBind使用教程 | dba.plotProfile