Chipseq上下游分析流程
步骤解释
Move fastq：管道的第一步是将原始fastq文件复制到指定的分析目录中。这样可以保持原始数据不变，以避免损坏或修改原始数据文件。
Concatenating：如果原始测序数据包含多个通道，此步骤是必需的，以便在分析之前合并这些通道。默认情况下，管道将所有fastq文件视为单个样本。
Trimming：可选的数据清理步骤。这允许使用trimmomatic15修剪低质量的读取或接头序列。
Aligner：默认情况下使用Bowtie216进行比对；也可以指定其他比对工具，如bwa-mem217。Bowtie2比对工具被选为默认工具，因为它特别擅长将相对较短的读取比对到相对较大的基因组，因此非常适合将ChIP-seq和ATAC-seq数据比对到哺乳动物基因组。为了避免任何中间文件，比对器直接传输到samtools view，以将bam文件保存到输出中。对于此规则，用户必须指定用于映射读取的首选基因组版本，例如hg19/hg38（人类），mm10/mm39（小鼠）。
Filtering：保留正确映射的读取，过滤掉低质量的读取。默认：samtools view，参数：-bShuF 4 -f 3 -q 30。
Sort：按最左边的坐标顺序对已比对的读取进行排序。默认：samtools sort（snakemake包装器），参数：-m 4G。
Mark duplicates：识别并标记所有重复读取。用户可以通过更改配置文件参数来决定是否删除它们。默认：Picard MarkDuplicates（snakemake包装器），参数：--REMOVE_DUPLICATES False，用于标记并保留重复项。
Merge bam：如果测序数据由重复或样本组成，可能希望合并为一个bam（即相同的样本在相同条件下多次测序）。在这种情况下，用户可以选择合并bams或在整个分析过程中保持bam文件分开。如果用户选择合并bams（使用samtools merge），必须为合并的bams指定一个公共前缀。
Index：此步骤对排序的坐标进行索引。默认：samtools index（snakemake包装器），使用samtools指定的默认参数。
BamCoverage：此规则从已比对的读取创建一个bigwig覆盖轨迹。使用deepTools中的bamCoverage工具，并计算每个窗口的读取数，其中窗口表示指定大小的范围。在此管道中，bamCoverage的默认参数设置为：-bs 1 -normalizeUsing RPKM -extendReads。
Peak calling：LanceOtron8被选为该管道的默认峰值调用器。与传统峰值调用器主要基于统计测试不同，LanceOtron是基于深度学习的峰值调用器，它结合了基因组富集测量和统计测试，并被证明在性能上优于行业标准峰值调用器MACS29。为了让bigwig与LanceOtron兼容，必须按碱基对计算覆盖度，并进行RPKM标准化；这在BamCoverage步骤的默认设置中有所体现。MACS2可以选择作为备用峰值调用器。将监测新峰值调用器的发布，并在适用时纳入，以保持和优化该分析管道的性能。可选MACS3
TrackDb：这创建了bigwig文件的键值对关联，以便在MLV6或UCSC基因组浏览器18等工具中加载和可视化它们。
创建环境
下载工具
conda create -n chipseq -y && conda activate chipseqmamba install -c bioconda -y bowtie2 fastqc multiqc trimmomatic Bowtie2 samtools Picard deepTools macs3pip install LanceOtron
下载数据
第一步下载fastq
mkdir -p fastq && >down && grep "HAoEC" ena_sra-experiment_20240917-0709.tsv | awk '{gsub("\"","",$1);print$1}'|grep -F -f - filereport_read_run_PRJNA532996_tsv.txt | awk -v file="down" '{gsub("fasp.sra.ebi.ac.uk:","",$3)}{print$3>file}' && ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -l 300M -T -P33001 -k1 --mode recv --host fasp.sra.ebi.ac.uk --user era-fasp --file-list down ./fastq下载GSE131953中的HAoEC内皮细胞的fastq数据
验证md5
>md5_resd && md5sum fastq/* > md5_rescomm -23 <(awk '{match($2,/SRR[0-9]+/)}{print substr($2, RSTART, RLENGTH)}' md5_res | grep -F -f - filereport_read_run_PRJNA532996_tsv.txt | awk '{print$2}' | sort) <(awk '{print$1}' md5_res|sort)
没有任何输出
或者
diff -s <(awk '{match($2,/SRR[0-9]+/)}{print substr($2, RSTART, RLENGTH)}' md5_res | grep -F -f - filereport_read_run_PRJNA532996_tsv.txt | awk '{print$2}' | sort) <(awk '{print$1}' md5_res|sort)
Files /dev/fd/63 and /dev/fd/62 are identical
表示fastq下载完整
判断单端测序还是双端
通过GEO的run SRA Run Selector 进入界面发现是单端测序
fastq质量评估
fastqc和multiqc
mkdir -p 1.qc 2.multiqc && nohub fastqc fastq/* -o 1.qc && multiqc 1.qc/* -o 2.multiqc &查看fastq是否已经过滤
查看一个样本的fastqc的质量报告
对应GEO数据查看原始数据报告
如果文章上有描述更好
质量报告重点
“每个碱基序列质量”图表是FastQC在ChIP-seq分析中最重要的模块；它提供了在每个读取位置上所有碱基质量分数的分布。这些信息可以帮助判断在数据测序过程中，测序是否出现了任何问题。一般来说，我们会发现reads末尾的质量会有所下降，但在reads的开头或中间部分不应该看到质量下降。若出现则需要咨询测序机构。
根据序列质量图，我们可以看到大多数读取的质量都很高，但有些部分却掉入了低质量区域，这表明相当数量的读取中存在低质量的碱基。我们可以对序列的两端进行修剪，或者使用一种对读取末端的低质量碱基可以忽略的对齐工具。
异常的fastqc结果可以看
http://bioinfo-core.org/index.php/9th_Discussion-28_October_2010
如何选择reads长短
由于开始测序时不准确，因此前15bp左右的ATGC含量有波动，需要在后面截掉。而最长reads为76bp，则暗示着我最后至少要保留多长的序列。
一般我的计算思路是：保留的序列长度 =（最长的的reads长度-前面波动的碱基长度）× 40%
去除接头，过滤低质量数据
trimmomatic
如果是双端测序
find ~/anaconda3/pkgs/trimmomatic-0.39-hdfd78af_2/ -regex ".*TruSeq3-PE.*" -exec cp {} ./ \;如果是单端测序
find ~/anaconda3/pkgs/trimmomatic-0.39-hdfd78af_2/ -regex ".*TruSeq3-SE.*" -exec cp {} ./ \;执行代码
mkdir -p 3.clean && for fq in fastq/*;do fq_f=$(basename $fq) trimmomatic SE -phred33 fq 3.clean/${fq_f} ILLUMINACLIP:TruSeq3-SE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 -threads 16;done
参数
Paired End:
对于大多数新数据集，可以使用适当的质量修剪和接头剪切。
通常不需要进行前导和尾随剪切。此外，通常将 keepBothReads 设置为 True 在处理配对末端数据时会很有用，您甚至会保留冗余信息，但这可能会让您的管道管理起来更简单。请注意，True 前面的额外 :2（用于 keepBothReads） - 这是回文模式下接头的最小长度，您甚至可以将其设置为 1。（默认值是非常保守的 8）
trimmomatic PE -phred33 input_forward.fq.gz input_reverse.fq.gz output_forward_paired.fq.gz output_forward_unpaired.fq.gz output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36HEADCROP:从 reads 的开头切掉指定数量的碱基，即从fastqc中看的（本次没有使用）
ILLUMINACLIP：去除接头（ILLUMINACLIP:TruSeq3-PE.fa:2:30:10）
LEADING: 从 reads 的开头切除质量值低于阈值的碱基（低于质量 3）（LEADING:3）
TRAILING: 从 reads 的末尾开始切除质量值低于阈值的碱基（TRAILING:3）
SLIDINGWINDOW: 从 reads 的 5' 端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗，使用宽度为 4 碱基的滑动窗口扫描读取，当每碱基的平均质量低于 15 时进行剪切（SLIDINGWINDOW:4:15）
MINLEN: 如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads，即上面的计算方法：（76-15）×40%≈24 （本次没有使用）
或者
测序片段的长度：MINLEN 应该与测序片段的长度相关。例如，如果片段的原始长度是 150 bp，那么设置 MINLEN 为一个较低值（例如 50-100 bp）可以帮助去除过短的片段，但保留足够长的高质量读段。
下游分析的需求：短的片段可能在比对时不准确，或影响后续的定量分析。因此，如果下游分析要求较长的片段，如比对或组装基因组时，建议选择相对较高的 MINLEN。例如，如果某些分析工具要求最小片段长度为 36 bp 或 50 bp，您应该确保 MINLEN 设置满足该要求。
Single End：
trimmomatic SE -phred33 input.fq.gz output.fq.gz ILLUMINACLIP:TruSeq3-SE:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36这将执行相同的步骤，使用单端接头文件。
比对+去重复PCR+过滤(取消#)
bowtie2
索引下载
bowtie2官网
也可以自己创建索引
mkidr -p reference/{hg38,hg19}cd reference/hg38## hg38的基因组文件wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gzgzip -d hg38.fa.gzbowtie2 build $(pwd)/reference/hg19/hg38.fa hg38 --threads 4
比对(下载的索引)+去重PCR+过滤(线粒体、核糖体等可选)
# 创建比对结果文件夹mkdir -p 4.align 5.filter 6.dedup# 遍历每个压缩的fastq文件for fq in 3.clean/*.gz; do# 提取文件名去掉扩展名fq_f=$(basename ${fq} .fastq.gz)# 1. 使用Bowtie2进行比对并使用samtools进行排序bowtie2 -p 8 -x ref/GRCh38_noalt_as/GRCh38_noalt_as -U $fq | samtools sort -O bam -@ 8 -o 4.align/${fq_f}.bam# 2.1 使用samtools flagstat统计映射比率samtools flagstat 4.align/${fq_f}.bam > 4.align/${fq_f}_flagstat.txt# 2.2. 统计唯一映射读取数（读取深度）并存储到文件中samtools view -q 1 -c 4.align/${fq_f}.bam > 4.align/${fq_f}_read_depth.txt# 3 进行bam质量查看(看下述安装过程)qualimap_v2.3/qualimap_v2.3/qualimap bamqc --java-mem-size=8G -bam 4.align/${fq_f}.bam -outdir 4.align/done#!/bin/bashfor fq in 3.clean/*.gz; dofq_f=$(basename ${fq} .fastq.gz)# 过滤+线粒体和核糖体去除# samtools view -h -F 4 -q 30 4.align/${fq_f}.bam \# | grep -v -e "chrM" -e "MT" -e "rRNA" \# | samtools sort -O bam -@ 10 -o - > 5.filter/${fq_f}_filter.bam# 确保输入BAM文件路径if [ "$(ls -A 5.filter)" ]; theninput_bam="5.filter/${fq_f}_filter.bam"elseinput_bam="4.align/${fq_f}.bam"fi# 去除重复readspicard MarkDuplicates \I=$input_bam \O=6.dedup/${fq_f}_dedup.bam \M=6.dedup/${fq_f}_dedup_metrics.txt \REMOVE_DUPLICATES=true \READ_NAME_REGEX=null# 质控分析qualimap_v2.3/qualimap bamqc --java-mem-size=8G -bam 6.dedup/${fq_f}_dedup.bam -outdir 6.dedup/# 使用samtools flagstat统计映射比率samtools flagstat 6.dedup/${fq_f}_dedup.bam > 6.dedup/${fq_f}_dedup_flagstat.txt# 排序并生成索引samtools index 6.dedup/${fq_f}_dedup.bamdone
注意index路径需要指定到index的前缀
qualimap安装
http://qualimap.conesalab.org/
解压之后测试，若出现如下报错，提示Java虚拟机（JVM）无法识别MaxPermSize选项，因为该选项已在较新的JVM版本中被废弃
解决：
打开文件qualimap
将如下代码
if [ "$JAVA_OPTS" != "" ]thenjava_options=$JAVA_OPTSecho "detected environment java options $java_options"elsejava_options="-Xms32m -Xmx$JAVA_MEM_SIZE -XX:MaxPermSize=1024m"fi
修改为
if [ "$JAVA_OPTS" != "" ]thenjava_options=$JAVA_OPTSecho "detected environment java options $java_options"elsejava_options="-Xms32m -Xmx$JAVA_MEM_SIZE"fi
测试
结果
参数解释
bowtie2参数解释
-p 8: 指定使用8个线程来并行化处理，以加快比对的速度。线程数通常根据可用的CPU核心数来选择。
-x ref/GRCh38_noalt_as/GRCh38_noalt_as: 指定 Bowtie2 使用的索引文件前缀。索引文件通常包含多个 `.bt2` 文件，前缀 `GRCh38_noalt_as` 是使用 `bowtie2-build` 构建索引时生成的。如果索引文件位于 `ref/` 目录下，并且前缀为 `GRCh38_noalt_as`，则这个路径是正确的。
-U $fq: -U参数指定单端测序数据的输入文件。这是单端测序的标志，若是双端测序则会使用 `-1` 和 `-2` 参数分别指定两端文件。
samtools sort 参数解释
-O bam: 指定输出的格式是 BAM 文件（与 SAM 文件不同，BAM 是二进制格式，更紧凑）。
-@ 8: 指定使用8个线程进行排序。这加快了排序的速度，特别是在处理大型 BAM 文件时。
-o 4.align/${fq_f}.bam: 指定输出文件的路径和名称。
samtools flagstat参数解释
- flagstat: Samtools 的 `flagstat` 命令用于生成 BAM 文件的比对统计信息。它会给出比对的总读取数、成功比对的读取数、重复率、质量过滤等信息。
输入文件，即上一步生成的排序后的 BAM 文件。
samtools view 参数解释
- b：输出BAM格式。
- S：输入文件是SAM格式。
- u：输出未压缩的BAM文件。
- F 4：过滤掉未比对的reads。
- f 2：只保留正确成对比对的 reads。
- f 3：只保留比对到正反链上的配对reads（properly paired）。
- q 30：过滤掉比对质量低于30的reads。
picard MarkDuplicates参数解释
I: 输入的BAM文件（即过滤后的BAM文件filtered.bam）。
O: 输出的BAM文件，标记了重复的reads（dedup.bam）。
M: 生成的重复reads标记的统计信息文件（dedup_metrics.txt）。
REMOVE_DUPLICATES=false: 保留重复reads，只进行标记，而不移除。
samtools index参数解释
生成一个.bai文件，与BAM文件配套使用，以加速下游分析如可视化或随机读取位置的检索
PCR重复
背景知识：
测序过程中的随机性：在文库制备过程中，DNA 被超声波或酶切断裂成随机大小的片段。理论上，每个片段的起始位置和长度都是随机的。
高表达量与测序深度：在 RNA-seq 中，高表达基因会产生更多的转录本，因此在测序中会捕获到更多来自这些基因的 reads。
怎么识别 reads 是 PCR 重复而不是真实的转录本？
1. 随机片段化与比对位置的概率
随机性与独立性：在理想情况下，DNA 片段化是完全随机且独立的。因此，不同的 DNA 片段（reads）比对到基因组上完全相同的位置（即起始位点和终止位点都相同）的概率非常低。
-概率计算：考虑到基因组的长度和片段的长度，两个独立片段完全相同的概率可以忽略不计。
2. PCR 重复的产生
PCR 扩增偏好性：在文库构建过程中，为了获得足够的 DNA 量进行测序，需要对片段进行 PCR 扩增。然而，PCR 过程存在偏好性，某些片段可能被过度扩增。
-结果：这导致了来自同一个原始 DNA 分子的多个拷贝（PCR 重复），它们在测序时会产生完全相同的 reads，比对到基因组的相同位置。
3. 高表达基因与 reads 数量
真实的生物学重复：高表达基因会产生更多的 mRNA，因此在 RNA-seq 中会捕获更多来自这些基因的 reads。
随机片段化的影响：即使某个基因被大量转录，产生了大量的 mRNA 分子，但在片段化过程中，这些 mRNA 分子被随机切割成不同的片段。因此，尽管来自同一基因的 reads 数量增加，但它们比对到基因组的位置仍然是随机分布的。
4. 为什么相同位置的重复 reads 被视为 PCR 重复
低概率事件：由于片段化和测序的随机性，不同的分子产生完全相同的片段（即起始位点和终止位点都相同）的概率极低。
判定标准：因此，当发现多个 reads 比对到完全相同的位置时，通常认为它们是由同一个 DNA 分子经过 PCR 扩增产生的重复，即 PCR 重复。
5. 高表达基因的特殊情况
概率增加但仍然低：对于高表达基因，确实有更多的 mRNA 分子被测序。然而，多个独立分子被随机切割并产生完全相同片段的概率仍然很低。
实践观察：实验和数据分析表明，即使在高表达基因中，真正的生物学重复 reads 通常比对到不同的位置，形成覆盖整个基因的 reads 分布。
6. PCR 重复的影响
数据偏差：PCR 重复会导致某些片段被过度代表，从而引入测序偏差。
定量准确性：如果不去除 PCR 重复，可能会夸大某些区域的 reads 数量，影响后续的定量分析。
7. 处理策略
DNA-seq 和 ChIP-seq：
去除 PCR 重复：在这些实验中，去除 PCR 重复是常见的，因为目标是识别基因组上的特定位点，重复 reads 可能导致假阳性。
RNA-seq：
谨慎处理：在 RNA-seq 中，去除 PCR 重复需要谨慎，因为高表达基因的 reads 数量确实更多。
策略选择：
保留重复：一些分析流程选择保留所有 reads，包括 PCR 重复，以更准确地反映基因表达水平。
标记而不去除：也有策略是标记 PCR 重复但不从数据中去除，或者在下游分析中考虑 PCR 重复的影响。
8. 总结
PCR 重复的识别依据：多个 reads 比对到完全相同的位置，极有可能是 PCR 重复。
高表达基因的影响：虽然高表达基因产生更多 reads，但这些 reads 比对的位置应当是随机且分散的，而不是集中在相同的位置。
分析决策：是否去除 PCR 重复应根据实验类型和分析目标决定。在 RNA-seq 中，需要权衡去除 PCR 重复对定量分析的影响。
bam文件查看
samtools view 4.align/SRR8925603.bam | head -n 2 | column -s$'\t' -tQNAME: 查询名称或reads名 - 这是在FASTQ文件头中存在的相同reads名称。
FLAG: 数值，用于提供有关reads mapping的信息以及该reads是否是对的一部分。
注意: FLAG中存储的信息是基于以下信息为真或假的叠加：
标志描述
Flag Description
1 read is mapped
2 read is mapped as part of a pair
4 read is unmapped
8 mate is unmapped
16 read reverse strand
32 mate reverse strand
64 first in pair
128 second in pair
256 not primary alignment
512 read fails platform/vendor quality checks
1024 read is PCR or optical duplicate
对于给定的比对，每个Flag都是开（真）或关（假），指示条件的真实性。FLAG是所有对该比对为真的单个标志（来自上表）的组合。标志值的妙处在于，任何标志的组合只能产生一个总和。
有一些工具可以帮助你翻译位标志，例如Picard提供的这个工具。
RNAME: 是参考序列名称，给出reads map的染色体。示例reads来自染色体1，这解释了为什么我们看到“chr1”。
POS: 指比对的左侧最左侧位置（以1为基准）。
MAPQ: 告诉我们比对的质量，其范围将取决于所使用的比对器。
CIGAR: 是一系列字母和数字，表示将读取匹配到参考所需的编辑或操作。字母是用于指示哪些碱基与参考对齐的操作（即匹配、不匹配、缺失、插入），数字表示每个“操作”的相关碱基长度。
操作描述
M序列匹配或不匹配
I向参考插入
D从参考删除
N跳过参考的区域
接下来的三个字段与配对末端数据更相关。
MRNM：是配对参考名称。
MPOS：是配对位置（基于1，最左侧）。
ISIZE：是推断出的插入大小。
最后，你还有原始FASTQ文件中每个reads的原始序列数据。
SEQ：是原始序列
QUAL：是每个位置的相关质量值。
samtools flagstat结果
14608455 0 in total (QC-passed reads QC-failed reads)
reads总数
37967 0 secondary
出现比对到参考基因组多个位置的reads数
0 0 supplementary
可能存在嵌合的reads数
0 0 duplicates
重复的reads数
14590894 0 mapped (99.88% : N/A)
比对到参考基因组上的reads数
14570488 0 paired in sequencing
属于PE read的reads总数。
7285244 0 read1
PE read中Read 1 的reads 总数。
7285244 0 read2
PE read中Read 2 的reads 总数。
14507068 0 properly paired (99.56% : N/A)
完美比对的reads总数。PE两端reads比对到同一条序列，且根据比对结果推断的插入片段大小符合设置的阈值。
14551500 0 with itself and mate mapped
PE两端reads都比对上参考序列的reads总数。
1427 0 singletons (0.01% : N/A)
PE两端reads，其中一端比上，另一端没比上的reads总数。
26260 0 with mate mapped to a different chr
PE read中，两端分别比对到两条不同的序列的reads总数。
17346 0 with mate mapped to a different chr (mapQ>=5)
PE read中，两端分别比对到两条不同
过滤黑名单区域(可选)
虽然不会执行这一步，但对我们的BAM文件应用额外的过滤通常是一个常见的做法。也就是说，我们会去掉发生在定义的黑名单区域的比对。
黑名单区域代表的是一些伪影区域，这些区域往往会显示出人为的高信号（过量的无结构异常读数映射）。这些区域通常存在于特定类型的重复序列中，如着丝粒、端粒和卫星重复，而且通常看起来是唯一可映射的，因此上面应用的简单可映射性过滤并不能去除它们。ENCODE和modENCODE联盟为多种物种和基因组版本（包括人类、小鼠、线虫和果蝇）编制了黑名单。这些黑名单区域（坐标文件）可以在继续进行峰值调用之前，从我们的比对文件中过滤掉。
bw文件生成
BamCoverage
mkdir -p 7.bamCoverfor i in 6.dedup/SRR89256*.bam;do f=$(echo $i | sed 's/6\.dedup\/\(.*\)\_dedup\.bam/\1/');bamCoverage -b $(pwd)/$i -o $(pwd)/7.bamCover/${f}.bw -bs 1 -p 15 --normalizeUsing RPKM;done
bamCoverage 是 deepTools 工具包中的一个命令行工具，主要用于从对齐的 BAM 文件生成覆盖率 (coverage) 轨迹，通常输出为 BigWig 格式文件。BigWig 文件是压缩的二进制格式，适合在基因组浏览器（如 UCSC Genome Browser）中可视化基因组数据。
bamCoverage的主要用途：
1. 生成覆盖率图 (Coverage tracks)：
- `bamCoverage` 计算每个基因组位置的覆盖率，即测序reads在该位置的覆盖深度。它可以通过指定窗口大小（bins）来平滑覆盖曲线。
- 例如：使用 `-bs 1` 参数将计算逐碱基对的覆盖率；如果设置更大的 bin size，如 `-bs 50`，则会对每50个碱基计算一次平均覆盖。
2. 归一化：
- 它支持多种归一化方式（如 RPKM、CPM、BPM 等），确保不同样本之间的覆盖深度可以进行合理比较。RPKM（每千碱基每百万reads的归一化）是常用的归一化方法。
- 示例命令：`bamCoverage -b sample.bam -o output.bw --normalizeUsing RPKM`
3. 生成 BigWig 文件：
- BigWig 文件比 BAM 文件更小且更易于在基因组浏览器中加载和可视化。
- 它们适合用于远程存储和快速可视化基因组数据，尤其是在 ChIP-seq、RNA-seq 等实验中，通过这些文件可以直观了解感兴趣区域的测序覆盖情况。
4. 扩展 reads：
- 在某些应用中（如 ChIP-seq 分析），可以使用 `--extendReads` 选项来扩展 reads 到估计的片段长度，从而更准确地反映片段的实际覆盖区域。
Peak calling
Peak calling
lanceotron --help
命令及其参数的详细解释：
1. Call Peaks (`callPeaks`)
此命令用于从 BigWig 轨迹中调用峰值。它通过计算每个窗口的滚动平均值来确定候选峰，并根据阈值过滤峰值。
参数解释：
- file:
- 需要分析的 BigWig 文件（即信号覆盖图）。
- t, --threshold:
- 候选峰的阈值。数值越高，表示对信号强度要求越高。默认值为 `4`，表示峰值信号强度需要高于该阈值才能被选为候选峰。
-w, --window:
- 滚动窗口的大小，窗口用于计算滚动平均值。默认值为 `400`，表示窗口大小为 400 bp。
-f, --folder:
- 输出结果保存的目录。默认值为 `./`（当前目录）。
--skipheader:
- 是否跳过输出文件的头信息。默认值为 `False`，表示会保留头信息。
2. Call Peaks with Input (`callPeaks_Input`)
此命令与 `callPeaks` 类似，但允许用户提供一个对照输入轨迹来计算峰值的显著性。这通常用于 ChIP-seq 分析中，控制样本作为对照以提高峰值检测的精度。
参数解释：
- file:
- 需要分析的 BigWig 文件（与目标样本对应）。
-i, --input:
- 控制输入轨迹的 BigWig 文件，用于比较分析以计算峰值的显著性。
-t, --threshold:
- 候选峰的阈值，默认值为 `4`。
- **-w, --window**:
- 滚动窗口的大小，默认值为 `400`。
-f, --folder:
- 输出结果保存的目录，默认值为 `./`。
--skipheader:
- 是否跳过输出文件的头信息，默认值为 `False`。
3. Score a Bed File (`scoreBed`)
此命令用于根据 BigWig 文件中的信号强度对 BED 文件中的现有区域进行评分，计算这些区域的信号覆盖率或其他相关统计数据。
参数解释：
- file:
- 需要分析的 BigWig 文件。
-b, --bed:
- 需要评分的 BED 文件，包含感兴趣的基因组区域。
-f, --folder:
- 输出结果保存的目录，默认值为 `./`。
--skipheader:
- 是否跳过输出文件的头信息，默认值为 `False`。
总结：
- `callPeaks` 和 `callPeaks_Input` 用于峰值调用，后者可以使用控制输入轨迹来提高准确性。
- `scoreBed` 用于对 BED 文件中的现有区域进行信号评分，常用于分析现有注释区域的覆盖情况。
mkdir -p 8.callpeakfor fq in 7.bamCover/*.bw;dolanceotron callPeaks $fq -f 8.callpeakdone
或者
首先整理好bam对应的-t和-c参数，即treat和control的bam文件(根据自己的数据进行整理，也可手动)
ls fastq/ | sed 's/.fastq.gz//g' | grep -F -f - filereport_read_run_PRJNA532996_tsv.txt| awk '{print$1,$4}'| awk 'BEGIN{FS=" "} NR==FNR{info[$1]=$2;next} FNR>1{gsub("\"","",$1);if($1 in info){print "6.dedup/"info[$1]"_dedup.bam","6.dedup/"$6"_"$8".bam"}}' - ena_sra-experiment_20240917-0709.tsv | sed 's/\"//g'| awk 'BEGIN{FS=" "}{system("cp " $1 " " $2)}'使用不同分组的前缀进行分析(根据自己的数据进行整理，也可手动)
mkdir -p masc3 && ls 6.dedup/EC* | sed 's/6\.dedup\/\(.*\)\_HAoEC_.*\.bam/\1/g' | uniq | while read info;do macs3callpeak -t $(pwd)/6.dedup/${info}_HAoEC_H3K27ac.bam -c $(pwd)/6.dedup/${info}_HAoEC_Input.bam -f BAM -g hs -n ${info}_HAoEC_H3K27ac -B -q 0.01 --outdir masc3; macs3 callpeak -t $(pwd)/6.dedup/${info}__HAoEC_H3K4me3.bam -c $(pwd)/6.dedup/${info}_HAoEC_Input.bam -f BAM -g hs -n ${info}__HAoEC_H3K4me3 -B -q 0.01 --outdir masc3;done
与MACS3的区别
1. MACS3：
- MACS3是基于统计学方法的传统峰值调用工具。它通过对目标（ChIP样本）与对照（input样本）的reads信号进行比较，找到目标样本中显著富集的区域。对照样本用于减少背景噪音，特别是在高背景信号时，通过对比降低假阳性峰值。
概述：
`callpeak` 是 MACS3 的主要功能。它可以处理多种格式的比对文件，并在基因组中调用显著富集的区域（即“峰”）。可以通过 `macs3 callpeak` 命令来调用此功能。如果您使用 `-h` 参数，则可以看到完整的命令行选项说明。以下只列出一些基本选项。
关键命令行选项
输入和输出选项：
-t / --treatment：
这是MACS3唯一必需的参数。该文件可以是任何支持的格式（参见 `--format` 选项的详细说明）。如果您有多个比对文件，可以通过 `-t A B C` 来指定。MACS3 将合并所有这些文件。
-c / --control：
控制文件、基因组输入文件或模拟IP（免疫沉淀）数据文件。与 `-t / --treatment` 选项的使用方式相同。
-n / --name：
实验的名称字符串。MACS3 将使用此名称创建输出文件，例如 `NAME_peaks.xls`、`NAME_peaks.bed`、`NAME_summits.bed`、`NAME_model.r` 等。请避免文件名与现有文件冲突。
-f / --format FORMAT：
标签文件的格式可以是 `ELAND`、`BED`、`SAM`、`BAM`、`BOWTIE`、`BAMPE` 或 `BEDPE`。默认值为 `AUTO`，这将允许 MACS3 自动确定格式。如果文件格式为 `BAMPE` 或 `BEDPE`，则必须明确指定。
--outdir：
MACS3 将所有输出文件保存到此选项指定的文件夹中。如果文件夹不存在，MACS3 会自动创建。
-B / --bdg：
如果启用此选项，MACS3 将保存片段堆积信息和控制lambda值为 `bedGraph` 文件。对于处理数据的文件，会保存为 `NAME_treat_pileup.bdg`，对于控制的局部lambda值，则会保存为 `NAME_control_lambda.bdg`。
--trackline：
MACS3 会在输出文件的头部添加轨迹信息，包括 `bedGraph`、`narrowPeak`、`gappedPeak` 和 `BED` 格式文件。如果需要将这些文件上传到UCSC基因组浏览器，可以使用此选项来包含轨迹信息。默认情况下，不会包含轨迹信息。
推荐格式的详细说明：
BED：
BED格式可以在UCSC基因组浏览器网站找到。输入的BED格式必须包含以下必要列：第1列是染色体名称，第2列是起始位置，第3列是结束位置，第6列是链方向信息。请注意，BED格式的坐标是以零为基数的半开区间。
BAM/SAM：
如果使用BAM或SAM格式，请参阅`samtools`的定义。如果BAM文件是为成对末端数据生成的，MACS3 只会保留左端（5'端）的标签。但是，如果指定了BAMPE格式，MACS3 将使用配对reads推断的实际片段来构建片段堆积。
BEDPE或BAMPE：
当指定为BAMPE或BEDPE格式时，会触发特殊模式。在这种情况下，MACS3 会将BAM或BED文件作为成对末端数据处理，而不是通过构建双模分布来预测片段大小，而是使用成对reads的实际插入大小来构建片段堆积。
BAMPE格式是包含成对末端比对信息的BAM格式，例如来自BWA或BOWTIE的比对结果。
BEDPE格式是简化的、更灵活的BED格式，包含前三列定义染色体名称、片段的左侧和右侧位置。请注意，这与bedtools使用的BEDPE格式不同，bedtools的BEDPE版本实际上不符合标准BED格式。
额外工具：
可以使用 MACS3 的子命令 `randsample` 或 `filterdup` 将包含成对末端信息的BAMPE文件转换为BEDPE格式文件。例如：
```bash
macs3 randsample -i the_BAMPE_file.bam -f BAMPE -p 100 -o the_BEDPE_file.bed
```
或者：
```bash
macs3 filterdup -i the_BAMPE_file.bam -f BAMPE --keep-dup all -o the_BEDPE_file.bed
```
Options controling peak calling behaviors
-g / --gsize
这是可测序的基因组大小或有效基因组大小，定义为可以测序的基因组大小。由于染色体上的重复特性，实际可测序的基因组大小通常比原始基因组大小要小，约为总大小的70%-90%。默认情况下，MACS3 为人类基因组推荐使用 `hs ~2.9e9`。
预编译的有效基因组大小参数（来自 deeptools）：
hs: 2,913,022,398（人类 GRCh38）
mm: 2,652,783,500（小鼠 GRCm38）
ce: 100,286,401（线虫 WBcel235）
dm: 142,573,017（果蝇 dm6）
-s / --tsize
序列标签的大小。如果不指定，MACS3 将尝试使用输入文件的前10个序列来自动确定标签大小。手动指定此值可以覆盖自动确定的值。
-q / --qvalue
调用显著区域的q值（最小FDR）阈值。默认值为0.05。对于宽泛的标记，可以尝试0.01。q值是通过Benjamini-Hochberg程序从p值计算得出的。
-p / --pvalue
p值的阈值。如果指定了`-p`，MACS3 将使用p值而不是q值。
--min-length, --max-gap
这两个选项可用于细化峰值调用行为。`--min-length` 指定被调用峰的最小长度，而 `--max-gap` 指定两个相邻区域之间的最大允许间距，如果间距小于此值，这两个区域将被合并为一个。如果未设置，MACS3 将默认使用预测的片段大小 `d` 作为 `min-length`，以及检测到的读长作为 `max-gap`。
--nolambda
启用此标志时，MACS3 将使用背景lambda作为局部lambda值。这意味着MACS3不会考虑峰候选区域的局部偏差。特别建议在没有控制样本的情况下调用峰值时使用此选项。
--slocal, --llocal
这两个参数控制在峰值区域周围检查局部lambda的两个区域大小。默认情况下，MACS3 使用1000bp的小区域(`--slocal`)和10000bp的大区域(`--llocal`)来捕捉来自长距离效应（如开放染色质区域）的偏差。
--nomodel
启用此选项时，MACS3 将跳过构建移位模型。您可以结合使用 `--extsize` 和 `--shift` 来实现所需的效果。
--extsize
当启用了 `--nomodel` 时，MACS3 使用此参数将读取扩展为固定大小的片段。此选项仅在 `--nomodel` 启用时有效，或当 MACS3 无法构建模型且启用了 `--fix-bimodal` 时有效。
--shift
此选项不是 MACS 版本2中的 `--shiftsize`，而是通过 `--extsize` 替换。您可以在这里设置任意的偏移量（以bp为单位）来调整整个库中的reads比对位置。此选项默认为0，推荐在ChIP-Seq数据集中保持为默认值。
--keep-dup
控制 MACS3 对待精确相同位置的重复标签（同一坐标和同一链）的行为。可以设置为 `auto`、`all` 或一个整数值。默认值为1，即最多保留一个重复标签。
--broad
启用此选项可调用宽峰，生成 UCSC `gappedPeak` 格式的结果，封装嵌套的峰结构。`--broad` 选项通过两个不同的阈值区分宽泛的弱峰和窄的强峰，后者随后嵌套在前者中。
--broad-cutoff
宽峰区域的阈值。仅在启用了 `--broad` 时可用。如果设置了 `-p`，这是一个p值阈值；否则，它是一个q值阈值。默认值为0.1。
--scale-to
设置为 `large` 时，将较小的数据集线性缩放到与较大数据集相同的深度。默认行为是将较大的数据集缩放到较小的数据集。
--call-summits
启用此选项后，MACS3 将重新分析信号曲线的形状，以分解大峰中的子峰。建议用于检测相邻的结合事件。
其他选项
--buffer-size
解释：MACS3 使用缓冲区大小来逐步增加内部数组大小，以存储每条染色体或每个contig的读取比对信息。增加缓冲区大小可以加快运行速度，但在某些染色体或contig中读取较少时会浪费更多内存。在大多数情况下，默认值 `100000` 已经足够。然而，如果比对文件中有大量染色体/contig，且每个染色体/contig 的读取量很少，建议指定较小的缓冲区大小以减少内存使用（但读取文件的时间将变长）。读取比对文件所需的最小内存量大约是：染色体数 缓冲区大小 8 字节。
默认值：100000
--cutoff-analysis
解释：启用此选项时，MACS3 将分析在不同截止值下可调用的峰数或峰的总长度，并输出一个总结表，帮助用户决定更合适的截止值。请注意，`minlen` 和 `maxgap` 可能会影响结果。与 `bdgpeakcall` 中的选项不同，`callpeak` 将同时执行调用峰的任务并生成截止值分析报告。
默认值：False
这些选项允许用户根据数据的大小和复杂性进一步优化MACS3的性能，特别是在内存使用和峰值调用截止值选择方面。
输出文件
1. NAME_peaks.xls
这是一个包含call peak信息的表格文件，您可以在Excel中打开并使用排序/筛选功能。信息包括：
染色体名称：显示峰所在的染色体。
峰的起始位置：峰的起点。
峰的结束位置：峰的终点。
峰区域的长度：峰的长度。
峰顶的绝对位置：峰值位置的绝对坐标。
峰顶的堆积高度：峰顶的信号强度。
峰顶的-log10(p值)：例如，p值为1e-10时，该值应为10。
峰顶的富集倍数：基于局部lambda的泊松分布计算的峰顶富集倍数。
峰顶的-log10(q值)：与p值类似，但使用的是q值。
注意，XLS文件中的坐标是1-based（从1开始），而BED格式则是0-based（从0开始）。如果使用宽峰调用模式，XLS文件中的堆积、p值、q值和富集倍数是整个峰区域的平均值，因为在宽峰模式下不会调用峰顶。
2. NAME_peaks.narrowPeak
这是一个`BED6+4`格式的文件，包含峰的位置以及峰顶、p值和q值。如果您计划将其加载到UCSC基因组浏览器，请确保启用了`--trackline`选项。特定列的定义如下：
第5列：显示的整数得分，计算方式为`int(-10 log10(p值))`或`int(-10 log10(q值)`，具体取决于使用的是p值还是q值作为得分的截止值。如果得分超出UCSC ENCODE narrowPeak格式定义的[0-1000]范围，可以使用以下一行代码将其限制在1000以内：
```bash
awk -v OFS="\t" '{$5=$5>1000?1000:$5} {print}' NAME_peaks.narrowPeak
```
第7列：峰顶的富集倍数。
第8列：峰顶的-log10(p值)。
第9列：峰顶的-log10(q值)。
第10列：峰顶相对于峰起始位置的相对位置。
3. NAME_summits.bed
这是一个BED格式文件，包含每个峰的峰顶位置。第5列与narrowPeak文件中的相同。此文件适用于寻找结合位点的motif（序列图案）。可以直接加载到UCSC基因组浏览器。
4. NAME_peaks.broadPeak
这是一个`BED6+3`格式文件，与narrowPeak文件类似，但缺少注释峰顶的第10列。此文件仅在启用了宽峰模式时生成。
5. NAME_peaks.gappedPeak
这是一个`BED12+3`格式文件，包含宽区域和窄峰。可以直接加载到UCSC基因组浏览器。
6. NAME_model.r
这是一个R脚本，您可以使用它根据数据生成模型图像。通过以下命令在R中加载它：
```bash
Rscript NAME_model.r
```
然后，它将在当前目录中生成一个PDF文件。
7. NAME_treat_pileup.bdg 和 NAME_control_lambda.bdg
这些文件是`bedGraph`格式，可以导入到UCSC基因组浏览器，也可以转换为更小的`bigWig`文件。`NAME_treat_pileup.bdg`包含来自ChIP/处理样本的堆积信号，而`NAME_control_lambda.bdg`包含来自控制样本的局部偏差估计值。
Cutoff值分析（Cutoff Analysis）
启用`--cutoff-analysis`选项时，MACS3将分析不同的p值截止值下可以调用的峰的数量或总长度，并生成一个总结表。报告包括以下四列：
score：可能的富集倍数截止值。
npeaks：在该截止值下可以调用的峰数。
lpeaks：所有峰的总长度。
avelpeak：峰的平均长度。
Cutoff值选择建议：
可以使用折线图分析（elbow analysis）找到显著改变趋势的截止值，或者根据期望的峰数或平均峰长来决定截止值。
2. LanceOtron
LanceOtron是基于深度学习的峰值调用工具，利用深度神经网络来预测信号富集区域。它通过训练模型来识别基因组中的峰值区域，而不是依赖传统的统计比较方式。LanceOtron结合了基因组富集测量和机器学习模型，模型经过大量训练，可以有效区分真实的峰值和背景噪音。
由于LanceOtron的深度学习模型已经在训练过程中学习了如何区分信号和背景噪音，因此它不需要通过提供`--control`（对照）样本来进行显著性比较。它的算法本身已经具有足够的鲁棒性来识别目标区域的富集信号。
主要区别：
- MACS3依赖对照样本来减少背景信号，而LanceOtron通过深度学习模型的训练，能够自动学会识别噪音和真实信号的差异，因此可以在没有对照样本的情况下仍然取得高精度的峰值调用结果。
- LanceOtron更侧重于使用数据驱动的模型来进行富集区域的预测，而MACS2则是依赖统计测试。
这使得LanceOtron在某些情况下能够更简化地调用峰值，尤其是在没有合适的对照样本时具有优势。
合并工具差异
合并工具差异
`bedtools intersect`、`bedtools merge`和`samtools merge`的功能及差异：
1.bedtools intersect
`bedtools intersect`是用来寻找两个或多个BED文件、GFF文件或其他基因组区间文件中共同区域（交集）的工具。
-功能：输出不同文件中重叠的区域。它可以精确定位在一个文件中的区域是否与另一个文件中的区域有重叠，并返回重叠的片段。
-用法：常用于寻找两个基因组注释文件中的共同区间，或者验证实验结果（如ChIP-seq、ATAC-seq等）是否在已知基因区域内。
-主要特点：
-支持多种文件格式（如BED、GFF、VCF等）。
-可以指定输出重叠区域的具体格式，支持仅输出A文件的重叠部分或B文件的重叠部分。
2.bedtools merge
`bedtools merge`用于将重叠或相邻的基因组区间合并为一个连续的区间。
-功能：合并同一BED文件中重叠或相邻的区间。特别适用于精简重复的区域或简化峰值文件中的区域。
-用法：如果两个或多个峰值在基因组上彼此相邻或重叠，使用`merge`可以将它们合并为一个连续的区间，避免重复计算。
-主要特点：
-仅用于处理一个文件，将相同染色体上重叠或相邻的区域合并。
-支持指定最大间隔来控制是否合并相邻但不重叠的区域。
3.samtools merge
`samtoolsmerge`是`samtools`中用于合并多个BAM文件的工具。
-功能：合并多个对齐后的BAM文件。这对于处理多个技术或生物学重复的样本时非常有用。合并后的BAM文件可以在后续分析中进行统一处理，如进一步的索引或比对质量检查。
-用法：当您有多个BAM文件时，例如从不同的测序运行中生成的，您可以使用`samtoolsmerge`将它们合并为一个文件以便进行统一的下游分析。
-主要特点：
-用于处理对齐数据，合并后的文件仍保留原始对齐信息。
-通常用于整合多个测序重复的BAM文件。
区别总结：
-bedtools intersect：用于找到两个或多个文件中的重叠区间（交集）。
-bedtools merge：用于合并同一文件中相邻或重叠的区间。
-samtools merge：用于合并多个BAM文件，将测序数据合并为一个文件。
下游分析思路
下游分析思路
第一步：基础分析
1. Peak/Reads在基因组上的分布
方法扩展：
Reads质量控制（QC）：在加载到IGV之前，使用FastQC和MultiQC对测序数据进行质量控制，评估Reads的质量，包括碱基质量分布、GC含量、接头污染等。
比对优化：选择适合的比对工具（如Bowtie2、BWA），根据数据特点调整比对参数，提高比对效率和准确性。
数据处理：使用Samtools、Picard等工具进行BAM文件的排序、标记重复、去除重复、索引等处理，确保数据的规范性和可用性。
多样本比较：在IGV中同时加载多个样本的数据，进行直观的比较分析，观察不同样本之间Reads和Peaks的分布差异。
实现扩展：
整合基因注释和其他功能元素：在IGV中加载基因注释文件（GTF/GFF）、已知的功能元件（如增强子、抑制子）等，提高可视化的丰富性。
识别异常区域：通过观察Reads覆盖度的异常高峰或低谷，可能发现基因组结构变异、重复区域或测序偏差等问题。
动态可视化：利用IGV的功能，制作动画或截图，用于展示特定基因或区域的Reads和Peak分布情况。
2. Peak在元件上的富集
方法扩展：
多层次注释：除了基本的功能元件，还可以考虑非编码RNA、lncRNA、miRNA基因、增强子（Enhancer）和超增强子（Superenhancer）等元素的富集分析。
统计显著性评估：使用超几何检验或Fisher精确检验评估Peak在各功能元件上的富集显著性，计算p值和q值。
背景模型构建：建立随机背景模型，评估富集的可信度，减少偏差影响。
实现扩展：
ChIPseeker高级功能：利用`plotAnnoPie`、`plotAnnoBar`、`plotAnnoDensity`等函数，生成富集结果的多种可视化图形。
功能元件数据库更新：使用最新的功能元件数据库，如ENCODE、FANTOM，确保注释的准确性和全面性。
比较分析：在不同条件或样本间比较Peak在元件上的富集差异，探索可能的调控机制。
3. Peak在基因元件上的分布
方法扩展：
细化基因结构分析：将基因结构细分为启动子、5'UTR、第一内含子、外显子、内含子、3'UTR等，分析Peak在各部分的分布情况。
转录本特异性：考虑可变剪接产生的不同转录本，对Peak在不同转录本上的分布进行分析。
基因功能关联：将Peak分布与基因功能结合，分析在特定功能基因上的Peak富集情况。
实现扩展：
使用Homer的annotatePeaks.pl脚本，生成详细的Peak注释报告。
绘制分布图：使用R的ggplot2或ChIPseeker的可视化函数，绘制Peak在基因元件上的分布柱状图、饼图或密度图。
比较不同条件：分析不同处理条件下，Peak在基因元件分布上的变化，寻找调控模式。
4. Peak的Motif分析
方法扩展：
多工具综合分析：使用多个Motif发现工具（如MEME Suite、Homer、DREME）进行分析，提高结果的可靠性。
已知Motif匹配：将发现的Motif与已知的转录因子结合位点数据库（如JASPAR、TRANSFAC）进行比较，预测可能的调控因子。
Motif位置偏好性：分析Motif在Peak内的位置分布，是否存在位置偏好。
实现扩展：
可视化Motif序列：使用序列Logo图展示Motif的碱基组成和保守性。
功能验证：结合已有文献或数据库，验证预测的转录因子是否与研究主题相关。
高级参数调整：在Motif分析中调整p值阈值、背景模型、序列长度等参数，优化结果。
5. Peak距离TSS位点的距离分析
方法扩展：
距离分布统计：绘制Peak到最近TSS的距离频率直方图或密度曲线，观察Peak的整体分布趋势。
分组分析：将Peak按照距离分组（如01kb、15kb、510kb、>10kb），比较各组的Peak数量和功能特征。
关联基因表达：分析不同距离范围内的Peak关联基因的表达水平，探讨距离与调控强度的关系。
实现扩展：
使用ChIPseeker的`plotDistToTSS`函数，生成距离分布图。
比较不同样本：在不同样本或条件下，比较Peak到TSS距离分布的变化，寻找调控模式的差异。
结合功能富集：对不同距离范围内的Peak关联基因进行功能富集分析，探讨生物学意义。
6. Peak修饰基因的功能分析
方法扩展：
多数据库富集：除了GO和KEGG，还可以使用Reactome、BioCarta、WikiPathways等数据库进行通路分析。
分子功能、细胞组分、和生物过程的深入分析：细化GO分析，分别探讨MF、CC、BP三大类别的富集结果。
疾病关联分析：使用DisGeNET、OMIM等数据库，分析关联基因与疾病的关联性。
实现扩展：
使用clusterProfiler的`enrichGO`、`enrichKEGG`、`enrichPathway`函数，进行多层次的功能富集分析。
可视化富集结果：利用`dotplot`、`emapplot`、`cnetplot`等函数，生成富集结果的可视化图形。
比较分析：比较不同条件下Peak关联基因的功能富集差异，寻找可能的生物学机制。
第二步：差异分析
1. 差异Peak鉴定
方法扩展：
选择合适的统计模型：根据数据特点，选择edgeR、DESeq2、csaw等工具，建立适合的统计模型，提高差异检测的准确性。
归一化处理：对ChIPseq数据进行归一化，如CPM、RPKM，以消除测序深度和样本间差异的影响。
多重检验校正：采用BenjaminiHochberg方法进行FDR校正，控制假阳性率。
实现扩展：
参数优化：调整差异分析的参数，如设定合适的Fold Change阈值、FDR阈值，平衡敏感性和特异性。
质量控制：对差异Peak进行可视化验证，如MA图、火山图，评估分析结果的可靠性。
生物学重复考虑：确保分析中纳入足够的生物学重复，以提高统计分析的可信度。
2. 非时序数据的分析策略
方法扩展：
比较组设计：合理设计比较组，确保处理组与对照组具有可比性，避免混杂因素影响。
批次效应校正：如果存在批次效应，使用limma的`removeBatchEffect`函数进行校正。
聚类分析：对差异Peak进行聚类，识别具有相似模式的Peak集合，可能代表共同的调控机制。
实现扩展：
热图绘制：使用pheatmap或ComplexHeatmap包，绘制差异Peak的信号强度热图，直观展示样本间的差异。
主成分分析（PCA）：评估样本间的全局差异，检查是否存在异常样本或分组不明显的情况。
差异Peak注释：对鉴定出的差异Peak进行功能注释，探索其生物学意义。
3. 时序数据的分析策略
方法扩展：
时间序列模型建立：使用线性模型、非线性模型或混合效应模型，捕捉Peak信号随时间的动态变化。
模式识别：识别不同的表达模式，如瞬时上调、持续下调、周期性变化等。
差异表达趋势分析：比较不同Peak在时间序列中的表达趋势，寻找关键的调控时点。
实现扩展：
使用maSigPro、ImpulseDE2等专门的时序数据分析工具，提高分析的敏感性和特异性。
可视化动态变化：绘制差异Peak在不同时间点的信号强度曲线，直观展示动态变化。
关联分析：将Peak信号的时间变化与基因表达的时间变化进行关联，探索调控关系。
4. 差异Peak关联基因的功能分析
方法扩展：
富集分析比较：比较处理组与对照组差异Peak关联基因的功能富集差异，寻找受影响的通路或生物过程。
上游调控因子预测：使用Ingenuity Pathway Analysis (IPA)等工具，预测可能的上游调控因子。
转录因子结合位点分析：分析差异Peak区域的Motif，预测可能的转录因子。
实现扩展：
使用clusterProfiler的`compareCluster`函数，比较不同条件下的富集结果。
网络分析：构建差异Peak关联基因的相互作用网络，识别关键节点和模块。
结合疾病信息：将差异Peak关联基因与已知的疾病相关基因库比较，探讨与疾病的潜在关联。
5. 差异Peak关联基因的PPI分析
方法扩展：
网络拓扑学分析：计算网络的拓扑参数，如节点度、介数中心性，识别网络中的关键基因。
模块检测：使用MCL、MCODE等算法，识别PPI网络中的功能模块。
网络功能富集：对模块中的基因进行功能富集，探索特定生物过程或通路的调控。
实现扩展：
使用STRING数据库，获取更全面的PPI信息，包括实验验证和预测的相互作用。
Cytoscape插件：利用CytoHubba、ClueGO等插件，进行网络分析和可视化。
动态网络分析：在时序数据中，构建不同时间点的PPI网络，观察网络结构的变化。
6. 感兴趣目标区域的可视化展示
方法扩展：
多层次信息整合：在IGV中叠加多个数据层，如ChIPseq信号、Peak位置、基因表达水平、变异信息等。
区域放大：对感兴趣的基因或区域进行放大观察，详细分析Reads分布和Peak形状。
定制化图形：使用Gviz、ggbio等R包，绘制高质量、可定制的基因组可视化图形。
实现扩展：
标注重要位点：在可视化图形中标注关键的Motif位置、突变位点等信息。
输出高分辨率图片：调整图形参数，导出适合发表的高分辨率图片。
交互式展示：利用JBrowse或UCSC Genome Browser创建交互式浏览器。
第三步：综合分析
1. Meta Genes整体关联
方法扩展：
信号强度聚合：对所有基因的上游、下游和基因体区域进行信号强度的聚合分析，观察全局的信号分布趋势。
基因分类分析：根据基因功能、表达水平、家族等，将基因分类，分别进行Meta Gene分析。
结合表观修饰：将ChIPseq信号与其他表观遗传修饰数据（如甲基化、乙酰化）进行整合分析。
实现扩展：
使用deepTools的`computeMatrix`和`plotHeatmap`、`plotProfile`等函数，生成热图和曲线图。
参数调整：优化窗口大小、步长、归一化方法等参数，获得更清晰的信号模式。
比较分析：在不同条件下，比较Meta Gene的信号分布，寻找全局性的差异。
2. Peak关联基因与DEG对应关联
方法扩展：
关联分析：计算Peak信号强度与基因表达水平之间的相关性，评估Peak对基因表达的影响。
共表达网络构建：基于关联分析，构建Peak关联基因与DEG的共表达网络。
多组学数据整合：将ChIPseq数据与RNAseq、ATACseq等数据进行整合，全面了解调控机制。
实现扩展：
Venn图展示：使用VennDiagram包，绘制Peak关联基因与DEG的交集和并集。
统计检验：使用超几何检验评估共有基因数量的显著性。
功能富集比较：对共有基因和特有基因分别进行功能富集，探讨功能差异。
3. 目标区域和靶基因的筛选
方法扩展：
多维度评分：根据差异倍数、显著性、网络中心性、功能重要性等因素，为候选基因打分。
机器学习方法：利用随机森林、SVM等机器学习方法，综合多种特征筛选关键基因。
实验验证优先级：根据评分和生物学意义，确定优先验证的靶基因和区域。
实现扩展：
文献调研：查阅候选基因的相关研究，了解其在特定生物过程或疾病中的作用。
设计验证实验：根据筛选结果，设计qPCR、Western blot、功能实验等验证方案。
下游分析工具
分析工具选择先看这篇文章
Comprehensive assessment of differential ChIP-seq tools guides optimal algorithm selection
下游分析工具
重复样本处理1.重复样本的交集peaks
和任何高通量实验一样，单个检测往往会受到相当大的变异影响。因此，强烈建议在实验设计中至少设置2-3个生物重复。理论上来说，两个测量相同基础生物学的重复应该具有很高的一致性，但实际上并不总是如此。为了评估重复之间的一致性，我们需要一些客观的指标来评估高通量检测的可重复性。
例如 (是相同状态下多次测序的才进行取交集, 是基于MACS3的narrowPeak结果)
$ bedtools intersect \-a macs2/Nanog-rep1_peaks.narrowPeak \-b macs2/Nanog-rep2_peaks.narrowPeak \-wo > bedtools/Nanog-overlaps.bed
找出两个重复实验（Nanog-rep1 和 Nanog-rep2）的峰值（peaks）之间的重叠区域。接下来，可以对这些重叠区域进行进一步分析，以评估数据质量、探索生物学意义。
1. 统计重叠峰的数量：评估两个重复实验之间的一致性。
重叠峰数=$(wc -l < bedtools/Nanog-overlaps.bed)
2. 计算重叠比例：了解重叠峰占各自实验总峰数的比例。
总峰数_rep1=$(wc -l < macs2/Nanog-rep1_peaks.narrowPeak)总峰数_rep2=$(wc -l < macs2/Nanog-rep2_peaks.narrowPeak)重叠比例_rep1=$(echo "scale=2; $重叠峰数 / $总峰数_rep1 * 100" | bc)重叠比例_rep2=$(echo "scale=2; $重叠峰数 / $总峰数_rep2 * 100" | bc)echo "重叠峰占 rep1 的比例：$重叠比例_rep1%"echo "重叠峰占 rep2 的比例：$重叠比例_rep2%"
3. 可视化重叠峰：通过基因组浏览器查看重叠峰在基因组上的分布。将重叠峰转换为适合浏览器的格式（如 BED、bigBed）。使用 IGV、UCSC Genome Browser 等工具加载并查看。
4. 功能注释：了解重叠峰对应的基因及其功能。
使用 ChIPseeker 进行注释：
install.packages("BiocManager")BiocManager::install(c("ChIPseeker", "TxDb.Hsapiens.UCSC.hg38.knownGene", "org.Hs.eg.db"))# 加载库library(ChIPseeker)library(TxDb.Hsapiens.UCSC.hg38.knownGene)library(org.Hs.eg.db)# 读取峰文件peaks <- readPeakFile("bedtools/Nanog-overlaps.bed")# 注释峰txdb <- TxDb.Hsapiens.UCSC.hg38.knownGenepeakAnno <- annotatePeak(peaks, TxDb=txdb, tssRegion=c(-1000, 1000), annoDb="org.Hs.eg.db")# 查看注释结果head(as.data.frame(peakAnno))
5. 富集分析：探究重叠峰相关基因的生物学功能和通路。
进行 GO 和 KEGG 分析：
# 提取基因列表gene <- as.data.frame(peakAnno)$geneId# 进行 GO 富集分析ego <- enrichGO(gene = gene,OrgDb = org.Hs.eg.db,keyType = 'ENTREZID',ont = "BP",pAdjustMethod = "BH",pvalueCutoff = 0.01,qvalueCutoff = 0.05,readable = TRUE)# 绘制结果library(enrichplot)dotplot(ego)
6. Motif 分析：识别重叠峰中的共有序列，预测潜在的转录因子结合位点。
# 创建输出目录mkdir motif_results# 运行 HOMERfindMotifsGenome.pl bedtools/Nanog-overlaps.bed hg38 motif_results/ -size given
7. 比较信号强度：评估重叠峰在两个重复实验中的信号一致性。
从原始的 narrowPeak 文件中提取信号值（如 fold change、p-value）。合并信号值，生成包含两个实验信号的表格。 使用 R 绘制散点图，分析相关性。
printf("hello wor# 读取原始峰文件rep1 <- read.table("macs2/Nanog-rep1_peaks.narrowPeak")rep2 <- read.table("macs2/Nanog-rep2_peaks.narrowPeak")# 提取重叠峰的信息overlaps <- read.table("bedtools/Nanog-overlaps.bed")# 合并信号值merged_data <- data.frame(rep1_signal=overlaps$V7, rep2_signal=overlaps$V17)# 绘制散点图plot(merged_data$rep1_signal, merged_data$rep2_signal,xlab="Rep1 信号强度", ylab="Rep2 信号强度",main="重叠峰信号强度比较")abline(lm(rep2_signal ~ rep1_signal, data=merged_data), col="red")ld!");
8.绘制Venn图展示重叠情况
library(VennDiagram)# 获取峰列表peaks_rep1 <- read.table("macs2/Nanog-rep1_peaks.narrowPeak")[,1:3]peaks_rep2 <- read.table("macs2/Nanog-rep2_peaks.narrowPeak")[,1:3]# 创建 Venn 图venn.diagram(x = list(Rep1 = peaks_rep1$V1, Rep2 = peaks_rep2$V1),filename = "venn_diagram.png",fill = c("red", "blue"),alpha = 0.5,cex = 2,cat.cex = 2,main = "峰值重叠情况")
可选：过滤重叠峰：
你可以根据重叠长度筛选峰，比如只保留重叠长度较大的区域，以确保峰的可靠性。例如，使用 awk 过滤掉重叠长度小于某个阈值的区域：
awk '$NF >= 100' bedtools/Nanog-overlaps.bed > bedtools/Nanog-overlaps-filtered.bed可选. 差异结合分析：
如果你有多个条件组的 ChIP-seq 数据，可以通过 DiffBind 包进行差异结合位点分析。通过结合这两个样本的峰信息，DiffBind 可以找出在不同条件下显著差异的结合位点。可以利用 .bed 文件和 BAM 文件进行分析
重复样本处理2.IDR
评估重复实验中峰值调用一致性的另一种方法是实施统计方法。一种流行的方法是由Qunhua Li和Peter Bickel的团队开发的IDR框架。它比较一对排名的区域/峰值列表，并分配反映其重复性的值。
基本的思想是，如果两个重复实验测量的是相同的基础生物学，那么最显著的峰值——很可能是真正的信号——在重复实验之间应该具有高度一致性，而显著性较低的峰值，更可能是噪音，在重复实验之间的则应该是一致性较低。如果将包含显著和不显著发现的一对排序列表（峰值）之间的一致性进行绘图，预计会看到一致性的变化（图1C）。这种一致性变化提供了从信号到噪音变化的内部指标，并暗示有多少个峰值被可靠检测到。
Why IDR?
IDR 避免了初始截止点的选择，这对于不同的调用者来说是不可比较的
IDR 不依赖于任意阈值，因此所有区域/峰值都会被考虑。
它是基于排名的，所以不需要将输入信号进行校准或使用特定的固定尺度（只需考虑顺序即可）。
IDR 方法创建了一条曲线，然后通过这条曲线定量评估在何时发现的结果在重复实验中不再一致。主要有三个组成部分：
对应曲线：在排名列表中，匹配峰值的图形表示。定性，不足以用来选择信号。
一种推断程序：总结可重复和不可重复信号的比例。采用定量方法，使用联合混合模型。
不可重现发现率（IDR）：通过推断过程（#2）以类似于FDR的方式导出显著性值，可用于控制选择信号时的不可重现率水平。也就是说，0.05 IDR表示该峰值有5%的概率是一个不可重现的发现。
IDR流程主要有三个步骤：
评估真实重复样本之间的峰值一致性
评估合并伪重复样本之间的峰值一致性
评估每个单独重复样本的自我一致性
要运行IDR，应该使用完整的数据集，完整的BAM文件。
运行IDR时，建议以较宽松的标准运行MACS3，以便为每个重复样本识别更多的峰值。此外，narrowPeak文件必须按-log10(p-value)列进行排序。
例如
$ idr --samples Nanog_Rep1_sorted_peaks.narrowPeak Nanog_Rep2_sorted_peaks.narrowPeak \--input-file-type narrowPeak \--rank p.value \--output-file Nanog-idr \--plot \--log-output-file nanog.idr.log
输出文件
输出文件格式类似输入文件类型，增加了一些额外字段。前10列是标准的narrowPeak文件，涉及两个重复样本的合并峰值。
第5列包含缩放的IDR值，即min(int(log2(-125IDR), 1000)。例如，IDR为0的峰值得分为1000，IDR为0.05的峰值得分为int(-125log2(0.05)) = 540，而IDR为1.0的峰值得分为0。
第11列和第12列分别对应局部和全局IDR值。
全局IDR是用于计算第5列中的缩放IDR值的，它类似于对p值进行多重假设校正以计算FDR。局部IDR则类似于峰值属于不可再现噪音成分的后验概率。
第13至16列对应重复样本1的峰值数据，第17至20列对应重复样本2的峰值数据。
后续功能查看IDR手册
差异peak
DiffBind 是一个 R 包，用于识别在两个或多个样本组之间不同结合位点。它主要处理call peak集（“peak集”），这些peak集是代表每个样本候选蛋白结合位点的基因组区间集合。它包含支持peak集处理的功能，包括在整个数据集中重叠和合并peak集、计算重叠区间中的测序reads，以及基于结合亲和力证据（通过reads密度差异测量）识别统计上显著的不同结合位点。
具体看DiffBind手册
数据可视化
ChIP-seq分析的第一部分使用常见的处理流程，包括将原始读数与基因组匹配、数据过滤和识别富集信号区域（峰值识别）。在第二阶段，各种程序允许对这些峰值进行详细分析、生物学解释和ChIP-seq结果的可视化。
有多种策略用于可视化富集模式，创建bigWig文件，这是一种常用于ChIP-seq数据可视化的标准文件格式。
将我们的比对文件（BAM）转换成 bigWig 文件。bigWig 格式是一种经过索引的二进制格式，适用于密集的连续数据，这些数据会以图形/轨道的形式在基因组浏览器中显示，同时也作为我们将在 deepTools 中运行的一些可视化命令的输入。
deepTools 是一套为高通量测序数据（如 ChIP-seq、RNA-seq 或 MNase-seq）进行高效分析而开发的 Python 工具。deepTools 提供了多种命令。
需要为BAM文件创建索引，以便快速检索与指定区域重叠的比对信息。将使用SAMtools工具套件中的samtools index命令来实现这一点，并通过for循环为每个BAM文件生成索引，以避免重复操作。
for file in ~/chipseq/results/bowtie2/*aln.bamdosamtools index $filedone
现在，为了创建 bigWig 文件，有两个工具非常有用：bamCoverage 和 bamCompare。前者将接受一个 BAM 文件并返回一个 bigWig 文件。后者允许将两个文件进行标准化（即 ChIP 样本相对于输入），并返回一个单一的 bigWig 文件。
使用 bamCoverage 命令为 Nanog 重复样本 2 创建一个 bigWig 文件。除了输入和输出文件外，还有一些我们添加的额外参数。
normalizeTo1x：报告标准化到 1x 测序深度的读取覆盖度（也称为每基因组内容的读取数 (RPGC)）。测序深度定义为：（映射读取的总数 * 片段长度）/ 有效基因组大小。因此，这里提供的数值代表有效基因组大小
https://deeptools.readthedocs.io/en/latest/content/feature/effectiveGenomeSize.html
可以找到一些常用生物的值的示例。
binSize：基因组片段的大小，以碱基为单位
smoothLength：定义一个超过 binSize 的窗口，用于平均读取数量。这有助于生成更连续的图。
centerReads：根据 extendReads 指定的片段长度将读取居中。这个选项有助于在富集区域周围获得更清晰的信号。
$ bamCoverage -b bowtie2/H1hesc_Nanog_Rep2_aln.bam \-o visualization/bigWig/H1hesc_Nanog_Rep2.bw \--binSize 20 \--normalizeTo1x 130000000 \--smoothLength 60 \--extendReads 150 \--centerReads \-p 6 2> ../logs/Nanog_rep2_bamCoverage.log
如果想创建一个大Wig文件，用输入文件进行ChIP归一化，会使用bamCompare。这个命令和bamCoverage非常类似，唯一的区别是你需要两个输入文件（b1和b2）。
## DO NOT RUN THIS$ bamCompare -b1 bowtie2/H1hesc_Pou5f1_Rep1_aln.bam \-b2 bowtie2/H1hesc_Input_Rep1_chr12_aln.bam \-o visualization/bigWig/H1hesc_Pou5f1_Rep1_bgNorm.bw \--binSize 20 \--normalizeTo1x 130000000 \--smoothLength 60 \--extendReads 150 \--centerReads \-p 6 2> ../logs/Pou5f1_rep1_bamCompare.log
轮廓图和热图
一旦有了bigWig文件，就可以用它们来查看指定区域的数据富集模式。在例子中，将评估TSS周围的富集情况，并分别为Nanog和Pou5f1样本（每个图中有两个重复）进行绘制。
为了节省时间，只查看12号染色体上的基因，而不是所有已知基因的TSS。将包含12号染色体上所有基因坐标的BED文件复制到可视化文件夹中。
computeMatrix命令可以接受多个bigWig文件和多个区域文件（BED格式），用于创建一个计数矩阵，作为中间文件。它还可以根据得分来过滤和排序区域。将指定一个在基因转录起始位点（TSS）周围的±1000bp窗口（-b和-a）。对于每个窗口，computeMatrix将根据bigWig文件中的读数密度值计算得分。
$ computeMatrix reference-point --referencePoint TSS \-b 1000 -a 1000 \-R ~/chipseq/results/visualization/chr12_genes.bed \-S /n/groups/hbctraining/chip-seq/full-dataset/bigWig/Encode_Nanog*.bw \--skipZeros \-o ~/chipseq/results/visualization/matrixNanog_TSS_chr12.gz \--outFileSortedRegions ~/chipseq/results/visualization/regions_TSS_chr12.bed$ computeMatrix reference-point --referencePoint TSS \-b 1000 -a 1000 \-R ~/chipseq/results/visualization/chr12_genes.bed \-S /n/groups/hbctraining/chip-seq/full-dataset/bigWig/Encode_Pou5f1*.bw \--skipZeros -o ~/chipseq/results/visualization/matrixPou5f1_TSS_chr12.gz \--outFileSortedRegions ~/chipseq/results/visualization/regionsPou5f1_TSS_chr12.bed
注意：通常来说，基因组区域指的是基因，可以通过UCSC表浏览器获取。或者，你可以查看其他感兴趣的区域，这些区域与基因组特征无关（例如，某个感兴趣蛋白的结合区域）。
使用该矩阵，可以创建一个轮廓图，基本上是一个密度图，用于评估所有转录起始位点的读数密度。对于Nanog，可以看到，与重复试验1相比，重复试验2在转录起始位点的信号特别强。
$ plotProfile -m visualization/matrixNanog_TSS_chr12.gz \-out visualization/figures/TSS_Nanog_profile.png \--perGroup \--colors green purple \--plotTitle "" --samplesLabel "Rep1" "Rep2" \--refPointLabel "TSS" \-T "Nanog read density" \-z ""$ plotHeatmap -m visualization/matrixNanog_TSS_chr12.gz \-out visualization/figures/TSS_Nanog_heatmap.png \--colorMap RdBu \--whatToShow 'heatmap and colorbar' \--zMin -4 --zMax 4
在差异富集区域可视化富集情况
之前，评估了研究中两个因素之间的差异富集。我们发现几乎所有识别出的峰值都是特定于Nanog的，只有一个区域在Pou5f1中显著富集。我们可以使用用DiffBind生成的BED文件作为deepTools的输入，来可视化这些区域的富集情况，以评估读取密度的差异。
$ computeMatrix scale-regions \-R ~/chipseq/results/visualization/Nanog_enriched.bed \-S /n/groups/hbctraining/chip-seq/full-dataset/bigWig/Encode_Pou5f1*.bw /n/groups/hbctraining/chip-seq/full-dataset/bigWig/Encode_Nanog*.bw \--skipZeros -p 6 \-a 500 -b 500 \-o ~/chipseq/results/visualization/matrixAll_Nanog_binding_sites.gz$ plotProfile -m visualization/matrixAll_Nanog_binding_sites.gz \-out visualization/figures/Allsamples_NanogSites_profile.png \--perGroup --plotTitle "" \--samplesLabel "Pou5f1-Rep1" "Pou5f1-Rep2" "Nanog-Rep1" "Nanog-Rep2" \-T "Nanog only binding sites" -z "" \--startLabel "" --endLabel "" \--colors red red darkblue darkblue
$ computeMatrix scale-regions \-R ~/chipseq/results/visualization/Pou5f1_enriched.bed \-S /n/groups/hbctraining/chip-seq/full-dataset/bigWig/Encode_Pou5f1*.bw /n/groups/hbctraining/chip-seq/full-dataset/bigWig/Encode_Nanog*.bw \--skipZeros -p 6 \-a 500 -b 500 \-o ~/chipseq/results/visualization/matrixAll_Pou5f1_binding_sites.gz$ plotProfile -m visualization/matrixAll_Pou5f1_binding_sites.gz \-out visualization/figures/Allsamples_Pou5f1Sites_profile.png \--perGroup --plotTitle "" \--samplesLabel "Pou5f1-Rep1" "Pou5f1-Rep2" "Nanog-Rep1" "Nanog-Rep2" \-T "Pou5f1 only binding sites" -z "" \--startLabel "" --endLabel "" \--colors red red darkblue darkblue
IGV进行质量评估
评估alignment质量的另一种方法是通过基因组浏览器进行可视化。将使用https://www.broadinstitute.org/igv/（IGV）。你应该已经在你的笔记本电脑上下载了这个工具。
IGV是一个互动工具，允许探索大型的综合基因组数据集。能够可视化不同的基因组数据表示方式，包括BAM（alignment）、BED（call peaks）、BigWig以及更多其他格式。此外，它还有一个大型的公共数据存储库可以使用。
Peak 注释
为样本找到了高度可信的peak，下一步就是对这些peak进行注释，以识别查询peak与基因/基因组特征之间的相对位置关系，从而获得一些生物学背景。
ChIPseeker的输入是一个或多个call peak 具体用法详见手册
annotatePeak函数默认使用TSS方法，并提供参数来指定最大距离截止值。还有一个选项可以报告在这个距离内的所有基因，无论是否与TSS重叠。对于注释基因组区域，annotatePeak函数在基因组区域是外显子或内含子时，会报告详细信息。例如，“Exon (uc002sbe.3/9736, 外显子80中的第69个外显子)”意味着该峰值与转录本uc002sbe.3所拥有的80个外显子的第69个外显子重叠，对应的Entrez基因ID是9736。
基因组区域注释
转录因子结合位点相对于转录起始位点的分布
这个图展示了每个峰值相对于转录起始位点（TSS）的位置。x 轴显示的是位点的百分比，而颜色则代表与 TSS 的距离。由于观察到几乎 50% 的峰值位于远端的基因间区域，因此大部分位点远离 TSS 也就不足为奇了。
Motif Analysis
MEME序列分析工具套件提供了大量的功能，用于基序分析和发现。DREME是一种基序发现算法，旨在寻找真核转录因子的短核心DNA结合基序，并优化以处理大型ChIP-seq数据集。
DREME专门针对真核数据，聚焦于短基序（4到8个核苷酸），涵盖了大多数真核单体转录因子的DNA结合区域。因此，它可能会因为大型转录因子复合体的结合而错过更宽的基序。
Tomtom判断识别出的基序是否与已知转录因子的结合基序相似，我们可以将这些基序提交给Tomtom，它会搜索已知基序的数据库，以找到潜在的匹配，并提供基序之间相似性的统计度量。
MEME-ChIP
MEME-ChIP 是 MEME 套件中的一个工具，专门用于 ChIP-seq 分析。MEME-ChIP 除了进行 DREME 和 Tomtom 分析外，还使用工具评估哪些基序最集中富集（基序应该位于峰中心），并将相关基序合并成相似性聚类。它能够识别更长的基序，长度小于 30bp，但运行时间会更长。
与RNAseq整合
由于基因表达调控机制的复杂性，多种组学数据的整合分析，从不同的层面探究生物问题越来越重要。从RNA-Seq层面，我们可以探究哪些基因具有显著差异，上调或下调；从ChIP-Seq层面，我们可以研究某个特定转录因子的调控作用；从ATAC-Seq我们可以了解到染色质可及性的动态变化，由于染色质的可及性与调控元件或转录因子的结合密切相关，在转录调控中发挥着重要的作用。因此，整合分析可以进一步探究调控某一生物学过程的关键因子（包括顺式调控元件和转录因子），以及哪个转录因子调控了感兴趣的基因，以及感兴趣的转录因子的靶基因等。
一是直接比较
首先得到差异基因与ChIP-Seq靶基因的overlap，然后选择一些关键基因比较一下谱图。
二是使用BETA工具
BETA (Binding and Expression Target Analysis)是 Shirley Liu实验室开发的工具，通过整合转录因子或染色质调控因子的ChIP-Seq与差异基因的表达直接预测靶基因，而且有可能发现增强子区的蛋白质的靶基因。
BETA有三个功能：
预测转录因子的功能是激活还是抑制
预测转录因子的靶基因
鉴定转录因子的motif以及调控转录因子激活或抑制的其他因子
BETA包括三个命令取决于输入数据格式和想要的输出数据
BETA basic: 预测转录因子的功能是激活还是抑制，直接检测靶标
BETA plus：BETA basic 靶标区域的motif分析
BETA minus：只基于结合数据的调控潜能的值预测TF靶基因
Chipseq学习
Chipseq学习
ChIP-seq分析（初学者的最佳链接之一）
https://github.com/crazyhottommy/ChIP-seq-analysis
使用R/Bioconductor进行NGS数据分析：ChIP-Seq工作流程
http://biocluster.ucr.edu/~rkaundal/workshops/R_feb2016/ChIPseq/ChIPseq.html
ChIP-seq预处理的管道
https://github.com/shenlab-sinai/chip-seq_preprocess
ChIP-seq数据分析
https://galaxyproject.org/tutorials/chip/
使用SPP和MACS2进行峰值调用，使用IDR进行一致性分析的ENCODE TF ChIP管道
https://github.com/mforde84/ENCODE_TF_ChIP_pipeline
ChIP-seq分析
https://www.ebi.ac.uk/training/online/course/ebi-next-generation-sequencing-practical-course/gene-regulation/chip-seq-analysis
ChIP-seq分析的实践性介绍 - VIB培训
http://www.biologie.ens.fr/~mthomas/other/chip-seq-training/
ChIP-seq数据综合分析的实用指南
http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003326
同学们赶快学起来吧
chipseq分析工具很多一定要结合文献和自己研究目的进行选择。
