CUT&Tag | Cut&Run | 数据分析 | deeptools为核心 | plotFingerprint | profile | heatmap

 

2023年05月04日

在搞懂DiffBind后就几乎没有用过deeptools，因为R可以用jupyter notebook，代码和结果方便撰写和保存。

但是，最近被怼了，质疑了我的peak calling，一个核心的问题就是macs call出来的peak太多了，有6万个左右，我觉得这个比较正常，但如何说服别人呢？【没啥毛病，cistrome也是这么多，想要严谨可以画个overlap图 score-public peak%】

参考：

 

plotFingerprint \
 -b *bam \
--labels H3K27me3 H3K4me1 H3K4me3 H3K9me3 input \
--minMappingQuality 30 --skipZeros \
--region 19 
--numberOfSamples 50000 \
-T "Fingerprints of different samples"  \
--plotFile fingerprints.pdf \
--outRawCounts fingerprints.tab

 

##################
computeMatrix scale-regions -S DMSO_H3K27Ac_1.normalized.bw M60a_H3K27Ac_1.normalized.bw \
                            -R /home/zz950/projects/cut_run/HDAC-b2/pipeline/6.peak/HDAC1_2.merged.narrowPeak \
                            -b 3000 \
                            --regionBodyLength 6000 \
                            -a 3000 \
                            --skipZeros -o H3K27Ac.matrix.mat.gz

plotProfile -m H3K27Ac.matrix.mat.gz  -o H3K27Ac.profile.pdf
# plotHeatmap -m matrix.mat.gz -o sample.heatmap.pdf

　　

 

 

 

最新的可以替代传统ChIP-seq的新技术。

 

传统ChIP-seq流程，了解一下：https://github.com/hbctraining/Intro-to-ChIPseq/blob/master/schedule/2-day.md

 

Cut&Run与ChIP-seq的核心区别

  • pA-Tn5的引入，增强了捕获的特异性，使得起始量需求变少，信噪比变低，需要的数据量变少。改善的信噪比意味着绘制染色质特征所需的 测序量减少了一个数量级

基于 Antibody-tethered Tn5-based 的方法由于 pA-Tn5 系结后严格的样品清洗和接头连接效率高，因此具有很高的灵敏度。相对于 ChIP-seq 而言，改善的信噪比意味着绘制染色质特征所需的 测序量减少了一个数量级，允许样本池 (通常高达 90 个样本) 通过文库的 条形码 barcoded PCR 在 Illumina NGS 测序仪上进行配对末端测序。

 

推荐使用的工具：

小工具，格式转化：

 

基本逻辑：

  • 比对bam file
  • call peak
  • bw-normalization（RPKM）

 

step1：scale bam file 【downsampling】

不确定效果好不好，相关性很差。需要继续评估。

【结论】经过实验，根本就不需要downsampling，macs3会自动矫正，只需要在normalization的时候考虑scale factor即可。详见HDAC的测试。

 

step2：SEACR peak calling【macs3足够好用】

有两种策略，一种是直接用control来call，另一种是对case和control分别call，下游再用DiffBind来处理。

 

step3：generate bigwig - smooth=30 & bin=10

 

step4：select best replicates for downstream analysis

【这个才是最重要的，原始的测序样本才是核心关键，这次处理的SOX9就是shit数据】 

 

step5: visualization (Heatmaps + Fold changes) 

【核心经验】DiffBind跑出流程不难，难就难在很难控制，也很难理解结果。

deeptools的结果高度可控，使用的normalization的方法也是很容易理解，所以我觉得下游所有的差异分析都可以基于deeptools的matrix数据。

就算是后面的单细胞的ATAC-seq数据，我觉得用这套方法来处理也是比较靠谱的。

 

 

一些小技巧，bam转bigwig去IGV可视化。【normalization的方法很重要】

bamCoverage --binSize 25 -p 4 --normalizeUsing RPKM --extendReads $fragLength -b $file -of bigwig -o qc/bw/$base.bw


file1=/home/da528/ATAC_Analysis/cut_run/cut_run1/cutoff25/AD3_H3K27Ac_2_CKDL220009611-1a_HN37TDSX3_L4_1_val_1.bowtie2.bam
bamCoverage --binSize 10 --smoothLength 30 -p 4 --normalizeUsing RPGC --effectiveGenomeSize 2150570000 --extendReads $fragLength -b $file1 -of bigwig -o raw.10x30.bw

file2=/home/da528/ATAC_Analysis/cut_run/zhixin_analysis/downsampling_bam/Exp1/AD3_H3K27Ac_2_CKDL220009611-1a_HN37TDSX3_L4_1_val_1.bowtie2.bam
bamCoverage --binSize 10 --smoothLength 30 -p 4 --normalizeUsing RPGC --effectiveGenomeSize 2150570000 --extendReads $fragLength -b $file2 -of bigwig -o downsampling.10x30.bw

　　

SEACR的两种peak，一种是没有对照的stringent call，另一个种是有对照的call。有spike-in就选non，没有就选择norm。

# first one
bash $SEACR_app $bed \
AD3_IgG_1_CKDL220009612-1a_HN37TDSX3_L1_1_val_1.bowtie2.bam.fragments.normalized.bedgraph \
     norm stringent seacr_control.peaks
# second one
bash /home/da528/ATAC_Analysis/cut_run/zhixin_analysis/software/SEACR/SEACR_1.3.sh AD3_SOX9_1_CKDL220019357-1A_H7MJHDSX5_L1_1_val_1.bowtie2.bam.fragments.normalized.bedgraph AD3_IgG_1_CKDL220019353-1A_H7MJHDSX5_L1_1_val_1.bowtie2.bam.fragments.normalized.bedgraph norm stringent Exp2_AD_SOX9_seacr_control.peaks

　　

指定geneset的peak可视化

cd /mnt/storage/home/da528/ATAC_Analysis/cut_run/cut_run1/cutoff25/qc/bw
out_dir=/mnt/storage/home/da528/ATAC_Analysis/cut_run/zhixin_analysis/deeptools

name1=Absc
peak1=/home/da528/ATAC_Analysis/cut_run/zhixin_analysis/DiffBind/sox9_paper_bed_zx/AbSC3.bed
color="#F8766D"

computeMatrix reference-point -S N3_IgG_1_CKDL220009620-1a_HN2WFDSX3_L1_1_val_1.bowtie2.10x30.bw AD3_IgG_1_CKDL220009612-1a_HN37TDSX3_L1_1_val_1.bowtie2.10x30.bw N3_SOX9_1_CKDL220009618-1a_HN37TDSX3_L4_1_val_1.bowtie2.10x30.bw AD3_SOX9_1_CKDL220009614-1a_HN37TDSX3_L4_1_val_1.bowtie2.10x30.bw  -R $peak1 --referencePoint center -b 500 -a 500 -bs=1 -p=max --outFileName $out_dir/$name1.matrix.gz
plotHeatmap -m  $out_dir/$name1.matrix.gz -o $out_dir/$name1.heatmap.png --legendLocation best --colorList "white,${color}" --samplesLabel N_IgG AD_IgG N_SOX9 AD_SOX9

　　

macs3 peak calling

将同类型【同抗体、同细胞】的数据全部merge起来call peak，以input control作为对照。【不同是指组织不同、TF或者histone不同】

MACS: Model-based Analysis for ChIP-Seq

指定--outdir，以及-n。

Example for regular peak calling on TF ChIP-seq:

macs3 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01

Example for broad peak calling on Histone Mark ChIP-seq:

macs3 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1

Example for peak calling on ATAC-seq (paired-end mode):

macs3 callpeak -f BAMPE -t ATAC.bam -g hs -n test -B -q 0.01

　　

ATAC-seq有自己的指定的call peak程序，hmmratac。

实操

macs2 callpeak -t ${projPath}/alignment/bam/${histName}_rep1_bowtie2.mapped.bam \
      -c ${projPath}/alignment/bam/${controlName}_rep1_bowtie2.mapped.bam \
      -g hs -f BAMPE -n macs2_peak_q0.1 --outdir $projPath/peakCalling/MACS2 -q 0.1 --keep-dup all 2>${projPath}/peakCalling/MACS2/macs2Peak_summary.txt

　　

 

 

 

 

官方建库和分析教程：Efficient low-cost chromatin profiling with CUT&Tag - Nature Protocols - 2022

参考：

 

posted @ 2022-11-24 04:51  Life·Intelligence  阅读(1716)  评论(0)  编辑  收藏  举报
TOP