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方向是否正确。
保存所有中间文件，以备后续分析。
