转录组数据对比之后是进行定量，而在chip-seq中后续的步骤是去除pcr重复，peaks calling，以及可视化。
本次分析步骤包括：环境部署——数据下载——查看数据(非过滤)——数据质控清洗——数据比对——去除pcr重复——peaks calling
分析步骤
1.延续上一个推文。数据比对，查看flagstat的比对结果
cat ./Z_Projects/chipseq/align/SRR30273120_flagstat.txt
45844006 + 0 in total (QC-passed reads + QC-failed reads)：表示总共有 45,844,006 个读段（包含通过质量控制的读段和未通过的读段），这里没有未通过质量控制的读段。 45844006 + 0 primary：这里的 45,844,006 个读段都是主要比对（primary alignment）。主要比对表示每个读段的主要比对位置。 0 + 0 secondary：没有次级比对（secondary alignment），次级比对指的是同一读段在基因组上其他位置的次要比对。 0 + 0 supplementary：没有补充比对（supplementary alignment），补充比对通常用于跨多个位置的读段（如嵌合读段）。 0 + 0 duplicates：没有重复读段，重复读段通常是指 PCR 扩增时产生的冗余片段。 0 + 0 primary duplicates：没有主要重复读段，这意味着在比对过程中没有检测到重复读段。 44185647 + 0 mapped (96.38% : N/A)：44,185,647 个读段（约占 96.38%）成功比对到参考基因组。这里的 N/A 表示没有适用的失败比对读段的比例。 44185647 + 0 primary mapped (96.38% : N/A)：同样地，有 44,185,647 个主要读段比对成功（约 96.38%），表示比对效果非常好。 0 + 0 paired in sequencing：没有成对的读段，表示这是 单端测序的数据。 0 + 0 read1：没有 read1，这与单端测序的特征一致。如果是双端测序，read1 表示配对中的第一个读段。 0 + 0 read2：没有 read2，也表明数据来自 单端测序。 0 + 0 properly paired (N/A : N/A)：没有正确配对的读段，因为这里使用的是单端测序数据，因此不适用。 0 + 0 with itself and mate mapped：没有读段和它的配对同时比对到基因组上，因为这是单端测序数据。 0 + 0 singletons (N/A : N/A)：没有 单独比对的读段（即只有一端比对成功），因为这里不是双端测序。 0 + 0 with mate mapped to a different chr：没有读段的配对比对到不同的染色体上，因为是单端测序数据，所以不适用。 0 + 0 with mate mapped to a different chr (mapQ>=5)：同样地，没有配对比对到不同染色体且比对质量分数大于等于 5 的读段，这是针对双端数据的统计。
2.去除PCR重复
# 创建文件
mkdir -p ${path}/rmdup
# 设定环境变量
export path="/home/lm/Z_Projects/chipseq"
nohup bash -c "
cat SRR_Acc_List.txt | while read id;
do
echo \"Processing \${id}...\"
samtools markdup -r \${path}/align/\${id}.bam \${path}/rmdup/\${id}_rmdup.bam
samtools flagstat \${path}/rmdup/\${id}_rmdup.bam > \${path}/rmdup/\${id}_rmdup_flagstat.txt
done
" > run_log.txt 2>&1 &
查看flagstat的比对结果
cat ./rmdup/SRR30273120_rmdup_flagstat.txt
3.peaks calling(查找DNA结合位点)
在peaks calling分析之前可以先重命名样本，这样做的目的是1、用生物学特征去替换SRA数据库ID号，有助于我们了解每个样本是什么并做过什么处理。2、也可以发现有多少样本存在生物学重复，为后续合并合并做准备。
cd ~/chipseq/rename
vim rename.tsv
# 把如下内容输入进去，中间的空格是tab键
SRR30273118 control
SRR30273119 p300
SRR30273120 BRD4
# 然后按下ESC键并输入:wq保存且退出rename.tsv文件
cat以下文件，查看内容使用上述内容修改样本文件名称
export path="/home/lm/Z_Projects/chipseq"
for i in {1..3}
do
id=$(cat rename.tsv | awk '{print $1}'| head -n $i | tail -1)
filename=$(cat rename.tsv | awk '{print $2}'| head -n $i | tail -1)
mv ${path}/rmdup/${id}_rmdup.bam ${path}/rename/${filename}.bam
samtools sort ${path}/rename/${filename}.bam -o ${path}/rename/${filename}_sort.bam
done
在这里会按照不同的命名以及sort每组分别生成两个文件如果一个样品分成了多个lane进行测序，那么进行peaks calling的时候，需要把bam进行合并，使用samtools merge, 知道即可，这个数据集不需要,以下代码在正式使用时需要修改！！！merge之后的文件需要重新构建索引。
B站视频中示例代码-merge数据-这里不采用这个
mkdir mergeBam
cd align
# 不用循环
# samtools merge control.merge.bam control_1_trimmed.bam control_2_trimmed.bam
# 通常的循环处理方式
# 这里使用“_”进行分割，但是针对不同的数据需要采用不同的分割方式！
ls *.bam | cut -d "_" -f 1 |sort -u |while read id do samtools merge
../mergeBam/$id.merge.bam $id.bam; done
推文中示例代码—笔者所分析的数据虽然不需要合并，但还是按照这个步骤去重新生成文件并构建索引
mkdir intersect
cd intersect
cat id_norep.tsv # 这个文件电脑上手动处理会更快
不需要合并数据代码如下
cp -r ./rename/*_sort.bam ./intersect/
path="/home/lm/Z_Projects/chipseq"
cat ${path}/intersect/id_norep.tsv | while read id;
do
index_file="${path}/intersect/${id}_sort.bam" # Target index file path
# Index the BAM file and specify the index file's location
samtools index $index_file
done > ${path}/indexing_output.log 2>&1 &
如果需要合并代码如下，但要进行修改
path="/home/lm/Z_Projects/chipseq"
cat ${path}/intersect/id_norep.tsv | while read id;
do
bedtools intersect -a ${path}/rename/${id}_rep1_sort.bam -b ${path}/rename/${id}_rep2_sort.bam > ${path}/intersect/${id}_intersect.bam
samtools index ${path}/intersect/${id}_intersect.bam
done
# bedtools intersect参数：
# -a 需要合并的文件1
# -b 需要合并的文件2
使用macs2正式进行peak calling分析
Peak Calling 分析 通常用于识别染色质免疫共沉淀测序（ChIP-Seq）、开放染色质测序（ATAC-Seq）或其他类似实验中的基因组区域（称为峰或peaks），这些区域通常代表了蛋白质-基因组相互作用位点或活跃的染色质区域。
这些峰是指在基因组特定区域中测序读数显著富集的区域，通过Peak Calling 方法，可以从背景噪声中提取这些生物学信号。
path="/home/lm/Z_Projects/chipseq"
cd ${path}
mkdir macs2
for i in {1..3}
do
id_p300=$(cat ${path}/intersect/id_norep.tsv | grep "p300" | head -n ${i} | tail -1)
id_BRD4=$(cat ${path}/intersect/id_norep.tsv | grep "BRD4" | head -n ${i} | tail -1)
macs2 callpeak -t ${path}/intersect/${id_p300}_sort.bam \
--control ${path}/intersect/control_sort.bam \
-g hs --outdir ${path}/macs2 \
-n ${id_p300}_intersect -q 0.01
macs2 callpeak -t ${path}/intersect/${id_BRD4}_sort.bam \
--control ${path}/intersect/control_sort.bam \
-g hs --outdir ${path}/macs2 \
-n ${id_BRD4}_intersect -q 0.01
done > ${path}/macs2_output.log 2>&1 &
# macs2 callpeak参数：
# -t 输入文件
# --control 对照组文件。
# -g 基因组大小。输入具体的数字（如1.0e+9或1000000000）。
# 对于人、小鼠、线虫或果蝇，可以分别用hs、mm、ce或dm。
# --outdir 输出文件的路径
# -n 文件名前缀
# -q Q值
bamCoverage
mkdir bamCoverage
path="/home/lm/Z_Projects/chipseq"
for i in {1..3}
do
id_p300=$(cat ${path}/intersect/id_norep.tsv | grep "p300" | head -n ${i} | tail -1)
id_BRD4=$(cat ${path}/intersect/id_norep.tsv | grep "BRD4" | head -n ${i} | tail -1)
bamCoverage --binSize 10 -p 15 --normalizeUsing RPKM -b ${path}/intersect/${id_p300}_sort.bam -o ${path}/bamCoverage/${id_p300}.bw
bamCoverage --binSize 10 -p 15 --normalizeUsing RPKM -b ${path}/intersect/${id_BRD4}_sort.bam -o ${path}/bamCoverage/${id_BRD4}.bw
done > ${path}/bamCoverage.log 2>&1 &
--binSize 10：设置 bin 的大小为 10 bp，这是计算覆盖度时使用的数据区块大小。 -p 15：指定并行处理的线程数为 15，提高处理速度。 --normalizeUsing RPKM：一共有RPKM,CPM,BPM,RPGC,None这几种方式，选择使用 RPKM 方法标准化数据。 -b {id_p300}_sort.bam 和 -b {id_BRD4}_sort.bam：指定输入的 BAM 文件，根据提取的 ID 动态指定。 -o {id_p300}.bw 和 -o {id_BRD4}.bw：指定输出的 BigWig 文件路径和名称。
既往推文
转录组上游分析流程(一）：转录组上游分析流程(一） 转录组上游分析流程(二）：转录组上游分析流程(二） 转录组上游分析流程(三）：转录组上游分析流程(三） 转录组上游分析流程(四）：转录组上游分析流程(四）
参考资料
生信技能树：https://mp.weixin.qq.com/s/KOFd60USQdKY2C4Sc4OIvA 生信学习日记：https://mp.weixin.qq.com/s/z-qbVZqo-cKT3Tfv4LT6OQ Tobin笔记：https://mp.weixin.qq.com/s/Uf8GgK4kZAWbRKKoh6IiyA 老俊俊的生信笔记：https://mp.weixin.qq.com/s/pfRiLhP8ONMNGcUI06IPog 生信媛：https://mp.weixin.qq.com/s/iHnogzM6QSu8tlImE26_Wg
致谢：感谢曾老师以及生信技能树团队全体成员。
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟
- END -