ChIP-seq | ATAC-seq | Cut&Run | Cut&Tag | DiffBind使用教程
核心可视化:
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里选择合适的样本。
结果一定要符合常识。
核心的可视化:
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就可以画图,非常高大上。
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方向是否正确。
保存所有中间文件,以备后续分析。