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

2025年01月15日

一个更直接的Heatmap出图方式,直接从bigwig读取数据。

教程:https://crazyhottommy.github.io/reproduce_genomics_paper_figures/07_figure1_i_j_k.html

https://crazyhottommy.github.io/reproduce_genomics_paper_figures/ 

 

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的问题

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
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;

 

参考:

 


 

2022年12月03日

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

  1. 样本的数量
  2. 送样的细胞数量;
  3. 下机的测序深度;

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

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

结果一定要符合常识。

 

核心的可视化:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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。

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
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 @   Life·Intelligence  阅读(2484)  评论(0编辑  收藏  举报
(评论功能已被禁用)
相关博文:
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
· spring官宣接入deepseek,真的太香了~
历史上的今天:
2021-12-03 医保目录谈判 | 国家带量采购 | DRG/DIP支付方式改革 | 腾笼换鸟
TOP
点击右上角即可分享
微信分享提示