Cut&Tag数据分析流程（一）

写在开头：在表观遗传学领域，染色质免疫沉淀（ChIP）是研究蛋白质（转录因子、组蛋白修饰等）与染色质相互作用的传统方法，可以分析特定蛋白在整个基因组上的分布情况。但是由于ChIP-seq需要较高的细胞起始量，这对于某些特殊的样品类型很难获得，比如哺乳动物原始生殖细胞。近年来，一类以融合蛋白proteinA-微球菌核酸酶(pA-MNase)和proteinA-Tn5转座酶(pA-Tn5)为核心的新型染色质靶向酶切测序技术兴起，并逐渐被用于低起始细胞量和单细胞的组蛋白修饰和转录因子调控网络研究。核酸酶的靶向切割和释放策略(cl-eavage under target and release using nuclease,CUT&RUN)以及靶向切割和标记策略(cleavage und-er target and release using nuclease,CUT&Tag)。
CUT&RUN的原理为利用特异抗体结合目标靶点蛋白，随后用pA-MNase结合抗体进行靶向切割，切割后产生的片段可通过自由扩散进入溶液上清，回收后即可用于建库测序。后续pAG-MNase融合蛋白，高钙/低盐酶切条件和残留大肠杆菌DNA归一化的方法的改进，使得CUT&RUN流程进一步成熟。与传统的X-ChIP-seq相比，CUT&RUN具有多种优点：①CUT&RUN可以在非交联的天然活细胞中原位进行，将细胞与伴刀豆球蛋白A固定以便于洗涤。②降低了起始细胞量要求。③CUT&RUN方法避免了超声处理和快速离心操作，更适合自动化。
CUT&Tag是CUT&RUN和Tn5转座酶技术的结合。它可以通过Tn5转座酶和特定测序接头序列的结合可以同时完成基因组的片段化和测序接头的引入，从而避免了极少量样品在后续二代测序文库构建时的困难。由于Tn5转座酶与MNase的性质不同，CUT&Tag选择与CUT&RUN不同的高盐反应缓冲液来抑制Tn5的非特异性切割活性，并且需要添加二抗来放大Tn5的信号。CUT&Tag在文库构建上更加容易操作，并且相比于CUT&RUN，CUT&Tag有更好的重复性，更低规模的细胞数量下限。

接下来我将介绍cut&tag数据上游分析流程
首先拿到公司测序得到的fastq文件（Raw Data)，上传到服务器的特定文件夹中，例如/home/zyn/cut&tag/fastq
质控

cd /home/zyn/
#注意mkdir -p ./cut&tag/fastq报错，是因为&符号在Linux shell中有特殊意义，它用于将一个命令放入后台执行。
mkdir ./cut-tag/fastq

# --------------------------->
# QC and pre-processing
# --------------------------->
mkdir -p /home/zyn/zbw_cut-tag/1.QC
cd /home/zyn/zbw_cut-tag/1.QC 
fastqc -t 10 ./fastq/*.fq.gz -o /home/zyn/zbw_cut-tag/1.QC
multiqc /home/zyn/zbw_cut-tag/1.QC

#预处理，过滤
#Trim Galore 工具要求在同一次运行中不能同时指定自定义序列和 Nextera 转座酶序列
#--quality参数用于指定截短低质量碱基的阈值。Trim Galore主要基于单个碱基的质量评分进行截短。它并不会直接使用FastQC图中的蓝线、黄色盒子或黑色胡须线作为判断依据，而是逐个碱基检查其质量评分。--length 75参数是去除掉小于75bp的reads，--adapter参数过滤掉overrepresented sequences，--nextera参数去除接头，这些参数会过滤掉很多reads，对于测序平台常见的adapter，trim_galore都会检测并去除的。对于大部分命令，只要不是自己特殊需要，都可以用默认参数。
#首先，对于fastqc得到的结果中存在大量的overrepresented sequences，先不要着急去掉，可能是具有生物学意义的，后面的命令再去除也可以，直接用trim_galore命令的默认参数过滤即可
vi samples.txt  #写入
C
C2
E
E2
cat samples.txt |while read a ;do echo "trim_galore --paired --gzip /home/zyn/zbw_cut-tag/fastq_last/${a}_1.fq.gz /home/zyn/zbw_cut-tag/fastq_last/${a}_2.fq.gz --fastqc -o ./";done  >trim.sh
cat trim.sh
#Mamba下载ParaFly，该工具可以并行运行脚本中的多个命令
nohup ParaFly -c trim.sh -CPU 4 &
zcat /home/zhangyina/zbw_cut-tag/2.trim_last/E_1_val_1.fq.gz |wc -l

比对

# --------------------------->
# mapping with bowtie2
# --------------------------->
mkdir -p /home/zyn/zbw_cut-tag/3.map/bowtie2_result
#--very-sensitive会产生更多比对结果，在比对过程中进行更细致的搜索，允许更多错配、插入或缺失，从而捕获那些在默认设置下可能被忽略的比对。
#样本DNA比对
cat ./sample.txt |while read id; do echo "bowtie2 --end-to-end --no-mixed --no-discordant --phred33 -I 10 -X 700 -p 10 -x /home/zyn/zbw_cut-tag/3.map_last/GRCm38_bt2ref/GRC38_PRI_index -1 /home/zyn/zbw_cut-tag/2.trim_last/${id}_1_val_1.fq.gz -2 /home/zyn/zbw_cut-tag/2.trim_last/${id}_2_val_2.fq.gz -S ./bowtie2_result/${id}_bowtie2.sam"; done > map.sh
nohup ParaFly -c map.sh -CPU 10 &

#用建库试剂盒给的spike in序列构建索引，spike in DNA比对
bowtie2-build Ecoli_904.fasta ./E.coli_index
mkdir -p /home/zyn/zbw_cut-tag/3.map/bowtie2_result
spikeinref="/home/zyn/zbw_cut-tag/3.map/E.coli_index/E.coli_index"

cat ./sample.txt |while read id; do echo "bowtie2 --end-to-end --no-mixed --no-discordant --phred33 -I 10 -X 700 -p 10 -x /home/zyn/zbw_cut-tag/3.map_last/E.coli_index/E.coli_index -1 /home/zyn/zbw_cut-tag/2.trim_last/${id}_1_val_1.fq.gz -2 /home/zyn/zbw_cut-tag/2.trim_last/${id}_2_val_2.fq.gz -S ./bowtie2_result/${id}_bowtie2_spikeIn.sam"; done > map_Ecoli.sh
nohup ParaFly -c map_Ecoli.sh -CPU 10 &

grep "rate" 3.map_last/nohup.out #搜索日志文件包含 "rate" 的行，并将这些行显示出来。这通常用于快速查看各个样本的比对成功率
#如果 spike-in DNA 和实验生物体基因组有相似的序列，可能会导致 spike-in DNA 的 reads 错误地比对到实验基因组上，或者实验基因组的 reads 错误地比对到 spike-in DNA 上。因此比对 spike-in DNA 时，特别添加 --no-overlap 和 --no-dovetail 参数，主要是为了确保 spike-in 序列与实验基因组序列之间的比对准确性，避免交叉比对
#--no-overlap：禁止成对的reads在比对时有重叠部分。也就是说，这个参数确保read1和read2在比对到参考序列时，不会有任何重叠的碱基对。避免那些实际上是相邻或部分重叠的reads被错误地比对到同一个区域，从而减少错配的可能性。
#--no-dovetail：禁止成对的reads形成“燕尾”状比对，即read1的5'端和read2的5'端相对，这意味着两个read的末端都比对到同一个位置，但方向相反。确保成对的reads比对时方向一致，不会导致一个read的5'端和另一个read的5'端比对到同一位置，从而减少错配的可能性。
#发现./spikein.out:0.00% overall alignment rate， 检查索引文件和输入文件似乎都有效，尝试去掉一些严格的参数，如--no-overlap 和 --no-dovetail 参数

#一般spike in 的比对率在1%左右，不要过多。0.09%也正常

比对结果过滤

# --------------------------------------->
# filtering and file format conversion
# --------------------------------------->
# sam转成bam，并且留下Q10以上的
#samtools view命令中，-F INT 参数用于过滤掉具有特定FLAG的reads。FLAG是SAM/BAM文件中的一个字段，用于描述read的各种状态，比如是否比对、是否是多重比对的一部分、是否是反向互补序列等。每种状态都有一个对应的值（用二进制表示），FLAG是一个整数，但实际上是多个二进制位的组合，每个位代表不同的状态。
#FLAG的常用值包括：0x4：该read未比对。0x8：该read的配对未比对。0x10：该read比对到参考序列的反向互补链。
#-bS: 输出BAM格式，-S表示输入文件是SAM格式

cat ./sample.txt |while read a;do echo "samtools view -bS -F 0x04 -q 10 /home/zyn/zbw_cut-tag/3.map_last/bowtie2_result/${a}_bowtie2.sam > ./${a}_bowtie2.q10.bam";done > filter.sh
nohup ParaFly -c filter.sh -CPU 10 &

cat ./sample.txt |while read a;do echo "samtools sort -@ 10 -o ${a}_bowtie2.q10.sorted.bam ${a}_bowtie2.q10.bam";done > sort.sh
nohup ParaFly -c sort.sh  -CPU 10 &

如果没加入spike in DNA，不需要计算scale factor进行标准化了
使用 Picard 标记去除重复

cat ./sample.txt | while read a;do echo "picard MarkDuplicates I=${a}_bowtie2.q10.sorted.bam O=${a}_bowtie2_redup.q10.bam M=${a}_duplication_metric.txt REMOVE_DUPLICATES=true"; done >remove_dup.sh
nohup ParaFly -c remove_dup.sh -CPU 10 &

cat ./sample.txt |while read a;do echo "samtools sort -@ 10 ${a}_bowtie2_redup.q10.bam -o ${a}_bowtie2_redup.q10.sort.bam";done > sort_redup.sh
nohup ParaFly -c sort_redup.sh -CPU 10 &

cat ./sample.txt | while read a;do echo "samtools index ${a}_bowtie2_redup.q10.sort.bam";done > bai.sh
nohup ParaFly -c bai.sh -CPU 10 &
#前面已按照序列在染色体上的位置进行排序
#MarkDuplicates 命令，标记重复的reads。输入排序后的文件，METRICS_FILE生成的统计信息文件。REMOVE_DUPLICATES=true，即可以去除重复

评估样品的再现性[相似性]，转换为bw文件方便导入IGV中进行可视化

###################评估样品的再现性[相似性]#############################
#当需要评估ChIP-seq类测序数据的相关性时，deeptools 是一个可行且方便的工具
#需要先进行RPKM标准化（无scale factor情况）
cat /home/zhangyina/zbw_cut-tag/4.filter_last/sample.txt |while read a;do echo "bamCoverage -b /home/zhangyina/zbw_cut-tag/4.filter_last/${a}_bowtie2_redup.q10.bam --normalizeUsing RPKM --binSize 1000 --extendReads 300 -o ${a}.1kb.RPKM.bw";done > narmalize.sh
nohup ParaFly -c narmalize.sh -CPU 12 &
##multiBamSummary / multiBigwigSummary 两个命令可以分别计算 bam/bw 文件在基因组或特定区域测序覆盖度。
#计算bw文件的coverage。对于转录组数据，通常通过样本的表达谱来计算样本间的相关性，对于chip-seq等没有明确定量结果的数据，通常的策略是将基因组划分为等长的区间，称之为bin, 计算每个区间内的覆盖度，然后通过比较不同样本间的覆盖度来计算样本相关性。发现bin为10kb时，再将spearman换为pearson，组内的相似性最好。
#-b输入bam或bw文件，-o 参数指定输出文件的路径和名称。.npz 是 NumPy 的压缩文件格式，存储的是 bw文件中每个窗口的读数，这个文件可以用于进一步分析或绘图。
#-bs 参数指定 bin 的大小，也就是每个窗口的大小（以碱基对为单位）。将基因组划分为的等大小区域，工具会计算每个区域中有多少 reads 覆盖。较小的 bin 尺寸会提供更高的分辨率，较大的 bin 尺寸则有助于更快的计算和更广泛的覆盖统计。
#--outRawCounts:用于将每个 bin 的原始读数输出到指定的文本文件中。输出文件是 readCounts.tab。该文件以表格格式保存，包含每个 bin 的基因组坐标及其在每个 BAM 文件中的读取覆盖数，这样可以进一步检查和处理这些原始数据。
nohup multiBigwigSummary bins -b C.10kb.RPKM.bw C2.10kb.RPKM.bw E.10kb.RPKM.bw E2.10kb.RPKM.bw -o results10kb.npz -bs 10000 --outRawCounts readCounts.10kb.tab &

##--corData，-in:由multiBigwigSummary或multiBamSummary生成的值压缩矩阵。--corMethod, -c： 计算相关性方法，spearman, pearson, Correlation method。
#--whatToPlot, -p： 图类型，选项有:heatmap, scatterplot。--plotFile, -o： 保存热图的文件格式是:.png， .eps， .pdf和.svg。--skipZeros： 排除所有样本中有零或缺失(nan)值的基因组区域(在所有样本中都没有信号的区域)
#--labels, -l: 用户定义的标签，而不是文件名称中的默认标签，多个标签必须用空格分隔，例如-labels sample1 sample2 sample3。--plotTitle, -T: 图标题
#--removeOutliers： 去除离群值。--outFileCorMatrix指定输出相关性矩阵的文件名称,这个文件包含样本之间的相关性值。--plotHeight： 图高 (cm)。(默认值:9.5)，--plotWidth： 图宽。最小值为1cm。(默认值:11)
#--zMin, -min： 热图强度的最小值。如果未指定，则自动设置该值。--zMax, -max：热图强度的最大值。如果未指定，则自动设置该值。
#--colorMap 参数指定热图的颜色方案。RdYlBu 是一个红-黄-蓝的渐变色彩映射方案。--plotNumbers:在热图的每个单元格中显示实际的相关性数值。此选项仅在绘制热图时有效。
#--xRange： X轴范围。默认的缩放会显示所有点的范围。--yRange：Y轴范围。默认的缩放会显示所有点的范围。--log1p：对数据执行 log(1 + x) 转换。注意，这仅用于绘图，相关性不受影响。
plotCorrelation -in results10kb.npz --corMethod pearson --skipZeros --plotTitle "pearson Correlation of Read Counts" --whatToPlot heatmap --colorMap RdYlBu_r  --zMin 0.75 --zMax 1.0 -o heatmap_pearsonCorr10kb_2.pdf --outFileCorMatrix pearsonCorr_readCounts10kb_2.tab

添加spike-in DNA的情况

计算有多少对reads比对到E.coli上，再写到一个文件里，方便后面计算scale factor

######计算有多少对reads比对到E.coli上，再写到一个文件里，方便后面计算scale factor##########
#遍历以_bowtie2_spikeIn.sam为结尾的sam文件。使用 samtools 查看每个 SAM 文件，-F 0x04过滤掉未成功比对的读段（即在 SAM 文件中标志位为 0x04 的读段）。wc -l：计算通过 samtools 输出的行数，即比对成功的reads数。
cat ./sample.txt |while read id; do echo "echo $(($(samtools view -F 0x04 ./3.map/bowtie2_result/${a}_bowtie2_spikein.sam|wc -l) /2)) > ./3.map/bowtie2_result/${a}_bowtie2_spikein.seqDepth";done > scale_factor.sh
nohup ParaFly -c scale_factor.sh -CPU 10 &
#由于在双端测序中，每对reads会有两行，将步骤 2 中计算出的reads数除以 2，以得到read对数。$((expression))计算算术表达式的值
# F 是 SAMtools 的选项，用于过滤出符合特定标志的记录。unmap 标志表示过滤出未比对到参考基因组上的读数。
#$(...) 是命令替换，用于在命令行中执行嵌套的命令并返回其输出。$((...)) 是算术扩展，用于执行数学运算。它将括号内的表达式计算并返回结果

对比对后的片段进行scale factor标准化，生成的bw文件以进行可视化，生成的bedgraph文件输入到seacr中识别峰值。同时标准化以后可以使用前面的deeptools工具评估样本重复性
基本假设是，对于使用相同细胞数量的一系列样本，映射到primary基因组和大肠杆菌基因组的片段比例是相同的。基于这一假设，不在各组之间或不同批次纯化pATn5之间进行标准化，因为这些批次可能包含非常不同量的大肠杆菌DNA。

# bc：一个命令行计算器工具。-l：选项让 bc 使用数学库，以支持浮点运算。
mkdir /home/zyn/zbw_cut-tag/4.filter/bedgraph
#自己构建chromsize文件
#利用samtools进行提取。samtools的faidx命令可以获取fasta文件中的序列长度信息，从其生成的后缀为fai的文件中可以获得chrom.sizes文件
samtools faidx GRCm39.primary_genecode.genome.fa
cut GRCm39.primary_genecode.genome.fa.fai -f1,2 > GRCm39.chrom.sizes
chromSize=/home/zyn/zbw_cut-tag/6.normalize/GRCm39.chrom.sizes

cat ./sample.txt | while read a
do
    echo "seqDepth=\$(cat ${a}_bowtie2_spikeIn.seqDepth)" >> bw.sh
    echo "if [[ \"\$seqDepth\" -gt \"1\" ]]" >> bw.sh
    echo "then" >> bw.sh
    echo "scale_factor=\$(echo \"10000 / \$seqDepth\" | bc -l)" >> bw.sh
    echo "echo \"Scaling factor for ${a}_bowtie2_spikeIn.seqDepth is: \$scale_factor\"" >> bw.sh
    echo "bamCoverage -b /home/zyn/zbw_cut-tag/4.filter/${a}_bowtie2_redup.q10.sort.bam -o /home/zyn/zbw_cut-tag/4.filter/${a}_bowtie2.fragments.normalized.bw --scaleFactor \$scale_factor --binSize 1000" >> bw.sh
    echo "fi" >> bw.sh
done
#使用 \$seqDepth 和 \$scale_factor 来确保变量引用在写入脚本时是有效的。使用\$、\"确保 bash 可以正确解析符号。

#seqDepth=cat $name``：读取当前文件的测序深度。
#在双方括号内，[[ ... ]] 表示进行条件判断，而 -gt 表示 "greater than"，用于比较两个数值的大小。"$seqDepth" 是要比较的变量，"1" 是用于比较的值。如果 $seqDepth 的值大于 1，条件判断返回真（true）。then：如果条件判断为真，则执行 then 下面的语句。fi：表示条件语句的结束
#将标准化系数 scaling factor:定义为S = C / (fragments mapped to E. coli genome) 这个常数是一个任意乘数，通常为 10,000。在 -l 参数下，bc 会加载一个标准数学库，这样就可以使用一些额外的函数和常量，例如三角函数等。这在进行更加复杂的数学计算时更加方便和灵活。
#然后计算Normalized coverage = (primary_genome_coverage) * S
#如果要生成bedgraph文件，设置-o output.bedgraph --outFileFormat bedgraph即可
#归一化基因组覆盖率（normalized genome coverage）是将基因组测序数据的覆盖率进行标准化处理，以便在不同样本之间进行公平比较
#文件内容包括染色体名称 起始位置 结束位置 归一化的覆盖率

识别peak。
相对于ChIP-seq数据来说，Cut&Tag数据的峰相对较窄，信号更集中，且背景噪音低。所以一些为 ChIP-seq 数据设计的适用于识别宽峰的工具，如sicer2，可能不适合处理cut&tag数据。SEACR是专为 CUT&Tag 和 CUT&RUN 设计的工具，能够更好地处理这些技术生成的高分辨率数据。MACS2是目前最常用的峰识别工具之一，适用于窄峰和宽峰的识别。
SEACR需要成对末端测序（paired-end sequencing）的bedGraph文件作为输入。SEACR在识别来自因子结合位点的窄峰值和某些组蛋白修饰特征的宽域方面都非常有效。

# ------------->
# peak calling
# ------------->
# SEACR,调用 SEACR 工具进行峰值识别（peak calling）
# SEACR包用于从非常低的背景（即没有reads覆盖的区域）的染色质中富集区域并call peaks。需要来自双端测序的bedGraph文件作为输入
#有或没有scale factor的情况下，都可以用这个命令生成seacr的输入文件.bedgraph，有scale factor的情况使用 --scaleFactor参数进行标准化。没有scale factor的情况下，使用 --normalizeUsing进行标准化
cat /home/zyn/zbw_cut-tag/4.filter_last/sample.txt | while read a; do
    echo "bamCoverage -b /home/zhangyina/zbw_cut-tag/4.filter_last/${a}_bowtie2_redup.q10.bam \
    --normalizeUsing RPKM \
    --binSize 1000 \
    --extendReads 300 \
    --outFileFormat bedgraph\
    -o ${a}.1kb.RPKM.bedgraph";done> bedgraph.sh
nohup ParaFly -c bedgraph.sh -CPU 12 &
#首先要先下载seacr.sh和Rscript,然后修改脚本，将目录换为自己的文件目录。还有保证安装了R和bedtools
# bash SEACR_1.3.sh <experimental bedgraph>.bg [<control bedgraph>.bg | <FDR threshold>] [norm | non] [relaxed | stringent] output prefix
# 加不加0.01差很多，把top1%作为peak
mkdir -p /home/zyn/zbw_cut-tag/7.peakcalling/
cat ./sample.txt |while read a;do echo "bash /home/zyn/zbw_cut-tag/7.peakcalling/SEACR_1.3.sh ${a}_1kb.RPKM.bedgraph 0.01 non stringent /home/zyn/zbw_cut-tag/7.peakCalling/${a}_seacr_top0.01.peaks ";done >seacr.sh
nohup ParaFly -c seacr.sh -CPU 12 &
#由于我们已经使用大肠杆菌read count对reads计数进行了标准化，或者RPKM进行标准化，所以我们将SEACR的标准化选项设置为“non”，non：表示不需要提供IgG对照组数据，适用于没有IgG对照组的情况，不会使用 IgG 对照组数据来排除背景信号。在选择“norm”选项时，需要提供实验数据和 IgG 对照组数据的 bedGraph 文件，SEACR 将使用 IgG 对照组数据来区分背景信号和实际富集区域。
#stringent：表示使用严格的峰值识别模式。FDR阈值为0.01
#文献中特定基因的峰图都是igv画出来的
#对于多个实验组和对照组应该如何设计，从而使得peak calling更加可靠，后续我会继续分享的。

########alternative method with MACS2（可选的另一种方法），组蛋白修饰一般是宽峰#################
# 单末端测序：给所有的读段延长到200/300bp的fragment
# 双末端测序：直接基于bam文件，用真正的长度。（前面sam转成bam，并且留下Q10以上的，过滤掉未成功比对的reads）。需要-f BAMPE 
mkdir -p /home/zyn/zbw_cut-tag/7.peakCalling/MACS2
cat /home/zhangyina/zbw_cut-tag/4.filter_last/sample.txt |while read a;do echo "macs2 callpeak -t /home/zhangyina/zbw_cut-tag/4.filter_last/${a}_bowtie2_redup.q10.sort.bam -g mm -B -f BAMPE -n ${a}_macs2_peak0.01 --outdir /home/zyn/zbw_cut-tag/7.peakcalling_last/MACS2 --broad -p 0.01 --keep-dup all";done > macs2.sh
nohup ParaFly -c macs2.sh -CPU 12 &
###macs2 callpeak：调用 MACS2 的命令进行峰值识别。-t $pair：指定实验组的 经过过滤的bam 文件作为输入，如果有多个可以用空格隔开输入-t A B C，如果没有对照组，也可以单独使用-t进行peak calling，但假阳性率会有所提升，通过spike in标准化可以不加对照组。-c：对照组的文件。-g mm：指定基因组的大小，通常会根据所研究的物种指定相应的基因组大小，MACS2给出了几个内置的基因组大小。-f BAMPE：指定输入文件的格式为 BAMPE，表示双末端测序的 bam 文件。
#-n：指定输出文件的前缀。--outdir：指定输出结果文件的目录。--broad：指定进行广峰值识别。-p 1e-5：指定阈值为 1e-5。--keep-dup all：保留所有的重复reads。在重复率低的情况下，--keep-dup all是合适的，确保所有数据都被利用。如果测序数据的重复率很高，这可能表明PCR扩增过程产生了大量的冗余reads，这些冗余reads可能不代表真实的生物学信号。保留这些重复的reads可能会导致峰值调用中的假阳性增加。

#当需要进行宽峰的peak calling时，可以使用以下参数 --broad：使用该参数后，MACS进行宽峰的callpeak。--broad-cutoff：宽峰的统计值cutoff。
#宽峰（Broad Peaks）：宽度通常超过1kb，有时可以延伸到几kb甚至更宽的区域，信号相对分散。如H3K27me3、H3K9me2、H3K36me3。窄峰（Narrow Peaks）：宽度通常在几十到几百个碱基对之间，信号集中在很窄的区域。

#--min-length：peak最短的长度，MACS2没有为--min-length设置默认值，它会自动识别并处理所有检测到的peak。如果不指定，MACS2不会强制要求peak的最小长度。
#--max-gap：定义了两个相邻peaks之间的最大距离，MACS2的--max-gap参数默认是根据数据自动调整的，具体取决于你使用的数据类型和模型。如果不设置，MACS2会根据数据的特性和统计模型自行决定相邻peaks合并的最大距离。
#--nomodel：禁用shifting model。在默认情况下，MACS2会使用shifting model，即它会自动估算并调整峰的延伸长度。这种模式是基于信号分布和背景噪声来推测峰的实际宽度，并进行调整，以更好地反映实际的生物学现象。使用--nomodel时，MACS2不会调整峰的宽度，而是使用用户指定的--extsize参数来定义peak的宽度。这对于已经有明确峰宽知识的实验，可能非常有用。
#--extsize：reads的延伸长度，在默认模式下（即没有使用--nomodel），必须手动指定--extsize。--extsize定义了从read的5'端向3'端延伸的碱基对数。如果你了解目标结合区域的宽度，可以手动设置这个值。
#-B/--bdg：生成bedGraph文件,如果不使用-B或--bdg参数，MACS2不会生成bedGraph格式的信号强度文件和对照文件。默认情况下，MACS2只输出peak文件，而不输出中间的信号文件。使用这个参数时，MACS2会输出bedGraph格式的文件，这些文件可以用于进一步的信号可视化和分析，如在IGV中查看信号分布。

用于识别基因组上的motif

#------------->
# motif
# ------------->
# MEME输入需要fasta格式
# 其他功能：想知道基因组上那些位置有某个motif，去基因组上scan
#用于motif分析的参考数据库和文件路径，包括motif数据库、参考基因组FASTA文件以及黑名单区域BED文件。黑名单区域通常是指在基因组中可能存在系统性误差或低质量比对的区域，需要在分析中排除。
motifdb=/public/workspace/shaojf/Course/NGS/Reference/MEME/JASPAR2018_CORE_vertebrates_non-redundant.meme
reffasta=/public/workspace/shaojf/Course/NGS/Reference/bowtie2_Index/GRCh38.fa
blacklist=/public/workspace/shaojf/Course/NGS/Reference/ENCODE/blacklist.v2/hg38-blacklist.v2.bed
bedtools intersect -u \
   -a /home/zyn/zbw_cut-tag/7.peakCalling/CTCF2_seacr_top0.01.peaks.stringent.bed \
   -b /home/zyn/zbw_cut-tag/7.peakCalling/CTCF1_seacr_top0.01.peaks.stringent.bed | \
   bedtools intersect -v -a - -b <(sed 's?chr??' $blacklist) > /home/zyn/zbw_cut-tag/7.peakCalling/CTCF.overlapped.bed
##bedtools intersect: 用于比较两个BED文件中的区域，进行交集或差集操作。是前面peak calling得到的BED文件
#第一个 bedtools intersect 命令使用 -u 选项输出第一个BED文件与第二个BED文件的交集，不包括重叠的具体信息，只保留A文件（第一个文件）中的原始行。-v: 这个参数表示只输出第一个输入文件（-a参数指定的文件）中不与第二个输入文件（-b参数指定的文件）中的区域相交的区域。
#-a 参数表示第一个输入文件，而 - 表示|前输出的数据，因为在 -a 后面没有指定具体的文件路径。-b <(sed 's?chr??' $blacklist) 指定了第二个输入文件（$blacklist所代表的 BED 文件，将其染色体名中的 chr 字符串删除，并将处理后的内容作为输入提供给 bedtools intersect 命令中的 -b 选项。
#使用了 <() 来处理第二个文件。sed 's?chr??' $blacklist 这个命令使用sed工具去掉了黑名单文件中的染色体名中的 chr 字符串，这样就可以和第一个文件中的染色体名匹配。<(sed 's?chr??' $blacklist)，s：表示替换命令。?chr??：表示将$blacklist文件中的 chr 替换为空，即删除 chr。<(...) 是一个 shell 特性，它允许将命令的输出作为一个文件进行处理

# 给坐标提取fasta
 bedtools getfasta -fi $reffasta -fo /home/zyn/zbw_cut-tag/5.motif/CTCF.overlapped.fa -bed /home/zyn/zbw_cut-tag/7.peakCalling/CTCF.overlapped.bed
#bedtools getfasta: 用于从参考基因组中提取给定的BED文件所对应的序列。-fi指定参考基因组的fasta文件，-bed 指定峰的BED文件，-fo指定输出的fasta文件路径。

meme-chip -o /home/zyn/zbw_cut-tag/5.motif/meme.CTCF.overlapped -db $motifdb /home/zyn/zbw_cut-tag/5.motif/CTCF.overlapped.fa
#meme-chip: 用于从ChIP-seq数据中发现motif。-db指定motif数据库的路径。输入前面得到的fasta文件路径。