参考文章：
ChIP-seq实践(H3K27Ac,enhancer的筛选和enhancer相关基因的GO分析)
一文讲明白ChIP-seq：高分文章里为什么做ChIP-seq?
ChIP-seq实践（非转录因子，非组蛋白）
一文读懂 ChIPseq
BED文件格式

一、测序数据质量控制

二、序列比对

比对的目的就是“推本溯源”，把我们的reads比对到参考基因组上，利用Bowtie2或这BWA看看我们过滤后的reads能匹配到基因组的什么位置。测序reads和基因组之间并非完全match上，中间会存在几个mismatch，有可能是因为测序错误，也有可能是存在变异位点。

Bowtie2

BWA

#比对
bwa mem -M /public/slst/home/BIO2083/reference/hg38/bwa_index/hg38 /public/slst/home/BIO2083/Class/Class7_ChIP-seq/rawdata/K562_H3K27ac.fastq.gz > K562_H3K27ac.sam
#转换为bam文件
samtools view -bS K562_H3K27ac.sam > K562_H3K27ac.bam 
#bam文件排序
samtools sort K562_H3K27ac.bam -o K562_H3K27ac_sorted.bam 
samtools sort K562_INPUT.bam -o K562_INPUT_sorted.bam
samtools index K562_H3K27ac_sorted.bam

#分别去除PCR duplicates得到K562_H3K27ac.deduped.bam和K562_INPUT.deduped.bam文件
/public/slst/home/BIO2083/software/bin/java -Xmx16g -jar /public/slst/home/BIO2083/software/bin/picard.jar MarkDuplicates -I K562_H3K27ac_sorted.bam -O K562_H3K27ac.deduped.bam -REMOVE_DUPLICATES true -ASSUME_SORTED true -VALIDATION_STRINGENCY LENIENT -VERBOSITY ERROR -M K562_H3K27ac_metrics.txt
/public/slst/home/BIO2083/software/bin/java -Xmx16g -jar /public/slst/home/BIO2083/software/bin/picard.jar MarkDuplicates -I K562_INPUT_sorted.bam -O K562_INPUT.deduped.bam -REMOVE_DUPLICATES true -ASSUME_SORTED true -VALIDATION_STRINGENCY LENIENT -VERBOSITY ERROR -M K562_INPUT_metrics.txt

/public/slst/home/BIO2083/software/bin/java -Xmx40g
#这部分指定了Java可执行文件的路径。Java是运行Picard工具的必要环境。-Xmx40g 是一个Java虚拟机（JVM）的选项，用于设置最大堆内存大小为40GB。这意味着Picard将能够使用最多40GB的内存来处理数据，这对于处理大型BAM文件非常重要。
-jar /public/slst/home/BIO2083/software/bin/picard.jar
#-jar 选项告诉Java虚拟机要运行一个打包成JAR（Java ARchive）文件的应用程序。
/public/slst/home/BIO2083/software/bin/picard.jar 
#是Picard工具集的JAR文件的路径。
MarkDuplicates
#这是Picard工具集中用于标记重复序列的命令。
-I K562_H3K27ac_sorted.bam
#-I 选项指定输入文件，这里是已经排序的BAM文件K562_H3K27ac_sorted.bam。BAM文件是存储测序读段（reads）的二进制格式文件。
-O K562_H3K27ac.deduped.bam
#-O 选项指定输出文件的名称，这里是处理后的去重BAM文件K562_H3K27ac.deduped.bam。
-REMOVE_DUPLICATES true
#这个选项设置为true，意味着在输出文件中将实际移除重复序列，而不仅仅是标记它们。
-ASSUME_SORTED true
#由于输入文件已经被排序，这个选项设置为true可以告诉Picard跳过排序步骤，从而提高处理速度。
-VALIDATION_STRINGENCY LENIENT
#这个选项设置Picard在验证输入文件时的严格程度。LENIENT模式意味着Picard会对输入文件中的小问题采取宽容态度，不会因这些问题而失败。
-VERBOSITY ERROR
#这个选项控制Picard在运行时输出日志信息的详细程度。ERROR级别意味着只有错误信息会被输出，这有助于减少日志文件的噪音。
-M K562_H3K27ac_metrics.txt
#-M 选项指定一个文件来存储关于重复序列的度量信息，这里是K562_H3K27ac_metrics.txt。这个文件将包含有关输入文件中重复序列的统计信息。

三、Peak calling

Peak calling用MACS寻找基因组中大量短读片段富集的区域。实际上表观组学的数据都会用到Peak calling这个概念，都是抓取特定区域的DNA片段，通过测序定量地看这些区域的reads数量，得到Peak 在基因组上的位置信息、peak 富集信息等等。
测序得到的 DNA 片段匹配映射到参考基因组，这些DNA片段其实是随机的，靶蛋白结合的片段越多，测序获得的数据就越多，那么在该位置检测到 DNA 片段堆叠就会越高，反之如果没有蛋白结合，在该位置就会几乎没有DNA 片段堆叠，将这些DNA片段堆叠用柱状图画出来，就会得到文章里出现的峰图 (Peak)。

横向代表基因组坐标，纵向代表ChIP-seq的信号强度，在文章中可能看到过有的峰图不只是向上的，还有向下的，水平线上方的峰代表正方向的,下方的代表互补链，有的一上一下有些错位，是测序造成的。转录因子的结合和组蛋白修饰，二者的峰形差异很明显：
转录因子结合的特征峰，峰型高，而且窄；而组蛋白修饰结合的特征峰，峰型起伏，而且宽。

macs2 callpeak --nomodel -B --SPMR -g hs -t /public/slst/home/BIO2083/student/hemiaomiao2024222204/Class7_test1/K562_H3K27ac.deduped.bam -c /public/slst/home/BIO2083/student/hemiaomiao2024222204/Class7_test1/K562_INPUT.deduped.bam --outdir /public/slst/home/BIO2083/student/hemiaomiao2024222204/Class7_test1/K562_H3K27ac -n K562_H3K27ac --broad

-t 指定实验组 BAM 文件。
-c 指定对照组 BAM 文件。
-f 指定输入文件格式为 BAM。
-g 指定基因组大小（这里使用 hs 代表人类基因组，但你可能需要根据实际情况调整）。
–qvalue-cutoff 设置 FDR 阈值（这里使用 0.00001，与图片中的 <0.00001 略有不同，但意思相近；你可能需要根据实际需求调整）。
–nomodel 表示不使用模型来估计背景 lambda。
–extsize 设置片段扩展大小（这里使用 147，但你可能需要根据实际情况调整）。
–call-summits 表示只调用峰顶位置。
-o 指定输出文件。
–SPMR：这个选项启用单峰模型（Single Peak Model with Random shift），它通过在每个读段（read）的范围内随机移动读段的起始位置来减少局部聚集效应，从而更准确地检测峰。
-n：指定输出文件的前缀名，所有输出文件将以此前缀命名。
–broad：这个选项表示输出的峰可能是宽泛的（broad peaks），这通常用于检测如H3K27ac这样的组蛋白修饰，这些修饰往往覆盖较宽的基因组区域。

```bash
#去除与blacklist overlap的region并且使⽤FDR < 0.00001进⾏peak筛选
bedtools subtract -a K562_H3K27ac_peaks.broadPeak -b /public/slst/home/BIO2083/reference/hg38/hg38-blacklist.v2.bed -A | perl -alne 'print if $F[8]>=5' > K562_H3K27ac_subtracted.bed

四、Peak annotation与可视化

所谓Peak注释，就是得到了靶蛋白在基因组区域的结合峰位置后，对峰位置进行注释。注释有两类，genomic annotation和nearest gene annotation：

• genomic annotation是看peak在基因组的位置，在各种基因组区域（基因上下游，5,3端UTR，启动子，内含子区）分布情况；
• nearest gene annotation是peak相对于转录起始位点（TSS）的距离，不管这个peak是落在内含子或者别的什么位置上，都能够找到一个离它最近的基因（即使它可能非常远），这种主要是应用于基因表达调控，因为启动子区域是重点，所以离TSS最近的基因更有可能被调控，所以这些peak区域附近的基因就作为其候选的调控基因。

目前启动子区域没有明确定义，在基因内部或距离TSS 2.5kb的peak被认为是靶基因。