多组学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日
基本分析流程都挺容易得,就是设计层面有些东西需要注意。
- 样本的数量
- 送样的细胞数量;
- 下机的测序深度;
根据目前经验,样本数和下机的测序深度对结果影响非常大【差异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就可以画图,非常高大上。
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方向是否正确。
保存所有中间文件,以备后续分析。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
· spring官宣接入deepseek,真的太香了~
2021-12-03 医保目录谈判 | 国家带量采购 | DRG/DIP支付方式改革 | 腾笼换鸟