chip-seq全流程

2024.9.23
染色体名称很重要一般都是chr1…这样
1.方面参考基因组文件中第一列要改，改完之后其索引也要改要不然之后进行获得peak之后提取序列时会出错
2.其次在索引时，如果是从biwtie2官网下载的索引文件，其染色体命名不是chr1…而是NC_050096.1…

我首先查看自己的bam文件中染色体名称是否规范

	#1.bam转sam
	samtools view -h OE_CHIP1.Aligned.sortedByCoord.final.bam > OE_CHIP1.sam
	samtools view OE_CHIP1.sam |head -n 2
	#-h: 表示输出包括头部信息
	#内容如下
	@HD VN:1.0 SO:unsorted
	@SQ SN:NC_050096.1 LN:1000000
	
	#不加-h如下
	samtools view OE_CHIP1.sam |head -n 2
	#内容如下
	E250020332L1C011R01702112136    147     NC_050096.1     259     30      143M    =       259     -143    ATGT.......
	
 	#可以看到染色体都是NC_050096.1这样的
 	#而正确的应该是这样的 chr1.......
 	#@HD VN:1.5 SO:unsorted GO:query
	#@SQ SN:chr1 LN:308452471
	#E250020332L1C011R01702112136 99 chr1 259 255 143M = 259 143 ATGTTTCCTTGT

————————————————————————————————————————————————————————

ChIP-seq 学习笔记总结

1. 数据过滤

目的：在 ChIP-seq 分析中，只保留编码蛋白质的基因，确保转录因子（TFs）只注释到基因的上下游区域，而不是非编码区域如核糖体。

步骤：

处理 GTF 文件：使用以下命令过滤 GTF 文件，保留编码蛋白质的基因信息。

grep "protein_coding" genes.gtf > genes.filter.gtf
#cat genes.gtf |wc
#1302223 30674924 354831452

#cat genes.filter.gtf |wc
#1284140 30281238 350388284

结果：生成新的 GTF 文件 genes.filter.gtf，只包含编码蛋白质的基因信息。
处理fa文件：ensamble上的参考基因组，染色体名称为1，2…，需要变为chr1…

sed 's/>/>chr/' file
 #我个人的真实数据中使用如下
 >chr1 dna:chromosome chromosome:Zm-B73-REFERENCE-NAM-5.0:1:1:308452471:1 REF
 >chrscaf_695 dna:scaffold scaffold:Zm-B73-REFERENCE-NAM-5.0:scaf_695:1:30084:1 REF
 #使用如下命令
 sed -i 's/dna:scaffold scaffold:Zm-B73-REFERENCE-NAM-5.0:scaf_[^:]*:[0-9]*:[0-9]*:[0-9]* REF//g' genome.fa
 sed -i 's/dna:chromosome chromosome:Zm-B73-REFERENCE-NAM-5.0:[^ ]* REF//g' genome.fa

2. 确定基因的可变剪切

目的：通过比较 GTF 文件中的基因数目和转录本数目，确定是否存在可变剪切（一个基因对应多个转录本）。

步骤：

统计基因和转录本数目：

awk '$3=="gene"' gene.gtf | wc -l
awk '$3=="transcript"' gene.gtf | wc -l
#awk '$3=="gene"' genes.gtf |wc
#44303 620242 5369267

#awk '$3=="transcript"' genes.gtf |wc
# 77341 1635426 18991023

结果：如果转录本的数量大于基因的数量，说明存在可变剪切。

3. 提取最长转录本

目的：从存在可变剪切的基因中提取出最长的转录本，简化后续分析。
步骤：
1. 使用 LongestTranscript.py 脚本：
```
srun -A xxxxx -J longesttranscript -N 1 -n 1 -p Debug python ../scripts/LongestTranscript.py genes.gtf genes_longest.bed6 genes_longest.gtf
```
- 读取 GTF 文件：解析文件中的基因和转录本信息。
  - 查找最长转录本：找到每个基因的最长转录本。
  - 保存结果：将最长转录本及相关信息保存为新的 GTF 文件，并转换为 BED6 格式。
1. 结果：生成新的 GTF 文件和 BED6 文件，包含每个基因的最长转录本信息。

4. 测序数据的过滤和质控

目的：去除低质量的 reads、去除有接头的 reads、去除过短的 reads，确保数据的高质量。

步骤：

处理低质量 reads 的两种方式：
- Remove：直接移除低质量的 reads。
- Trim：修剪低质量部分，保留较高质量部分。

使用 fastp 进行质控和过滤：

基本命令结构如下：

fastp -i input_directory/input_file.fastq.gz \
  -o output_directory/output_file.fastq.gz \
  --html output_directory/output_report.html \
  --json output_directory/output_report.json \
  --thread num_threads \
  --adapter_sequence adapter_sequence_r1 \
  --adapter_sequence_r2 \
  --qualified_quality_phred quality_threshold \
  --length_required minimum_length \
  1>log_directory/quality_control.log 2>&1

参数解释：
- -i：输入的 FASTQ 文件路径。
- -o：输出的过滤后 FASTQ 文件路径。
- --html：生成 HTML 格式的质控报告。
- --json：生成 JSON 格式的质控报告。
- --thread：指定用于处理的线程数量。
- --adapter_sequence：指定接头序列。
- --qualified_quality_phred：指定用于判断合格的质量阈值（Phred 质量分数）。
- --length_required：指定 reads 的最小长度，低于此长度的 reads 会被丢弃。
- 1>log_directory/quality_control.log 2>&1：将标准输出和标准错误重定向到日志文件。

   #批量操作
   $ ls ../data/raw/*.gz | xargs -n 2 | while read i j; do \
   base_name=$(basename "$i" | cut -d '.' -f 1) \
   echo fastp --thread 4 -i "$i" -I "$j" -o   "${base_name}_1.fq.gz"  \
     -O "${base_name}_2.fq.gz" --adapter_sequence 	AGATCGGAAG \
     --adapter_error_rate 0.2 --n_base_limit 10 --qualified_quality_phred 10 \
     --qualified_percentage 50 ;done
    ## fastp --thread 4 -i ../data/raw/OE_CHIP1.R1.fq.gz -I ../data/raw/OE_CHIP1.R2.fq.gz -o OE_CHIP1_1.fq.gz -O OE_CHIP1_2.fq.gz --adapter_sequence AGATCGGAAG --adapter_error_rate 0.2 --n_base_limit 10 --qualified_quality_phred 10 --qualified_percentage 50
     ##fastp --thread 4 -i ../data/raw/OE_CHIP2.R1.fq.gz -I ../data/raw/OE_CHIP2.R2.fq.gz -o OE_CHIP2_1.fq.gz -O OE_CHIP2_2.fq.gz --adapter_sequence AGATCGGAAG --adapter_error_rate 0.2 --n_base_limit 10 --qualified_quality_phred 10 --qualified_percentage 50

5. 合并 fastp 的 JSON 报告文件

目的：将多个 fastp 生成的 JSON 质控报告文件合并为一个 TSV 文件，便于统一查看各样本的质控结果。
步骤：
1. 使用 R 脚本合并 JSON 文件：
  - 提取关键质控指标：如 total_reads, total_bases, q30_rate, gc_content。
  - 保存结果：生成合并后的 TSV 文件和 RData 文件。
2. 结果：生成一个汇总的 TSV 文件，包含所有样本在过滤前后的质控指标。

6. 查看 TSV 文件

工具：在 Linux 系统中使用 csvlook 工具查看 TSV 文件，使数据更易于阅读。

7. 比对到参考基因组

目的：将过滤后的 reads 比对到参考基因组，获取其在基因组上的定位信息。
常用工具：
- DNA 数据：通常使用 BWA 或 Bowtie2 进行比对。
- RNA 数据：通常使用 STAR 或 HISAT2 进行比对。当然，HISAT2 和 STAR 也可用于 DNA 数据的比对。

8. 构建参考基因组

步骤：
1. 使用 bowtie2-build 构建参考基因组索引：
```
bowtie2-build --threads num_threads input_reference.fa output_index_prefix \
  1>log_directory/bowtie2_build.log 2>&1
```
2. 下载已构建的索引：可以直接从官网上下载已经构建好的参考基因组索引，省去自行构建的时间。
  使用迅雷会快一些，之后对下载的索引重命名
```
 ls *.bt2 | xargs -n 1 | while read i; do mv ${i} ${i/Zm-B73-REFERENCE-NAM-5.0/genome}; done
```

9. 进行比对

步骤：

使用 bowtie2 进行比对：

bowtie2 -p num_threads \
  -x input_index_prefix \
  -U input_reads.fastq.gz \
  -S output_alignment.sam \
  1>log_directory/alignment.log 2>&1

结果：生成比对后的 SAM 文件，可以使用 less 命令查看内容。

10. SAM 文件处理

步骤：

将 SAM 文件转换为 BAM 文件：

samtools sort \
  -@ num_threads \
  -o output_sorted.bam \
  input_alignment.sam \
  1>log_directory/sort_bam.log 2

&1
```

为 BAM 文件构建索引：
```
samtools index output_sorted.bam
```

11. BAM 文件过滤

目的：由于 ChIP-seq 的起始量较高，可能导致文库质量过高，因此需要对 BAM 文件进行进一步过滤。
常用工具：
- Picard：标记 duplicates。
- Samtools：去除各种低质量比对。
- DeepTools：使用 alignmentSieve 一步完成上述操作，并且速度更快。
过滤内容：
- duplicates：经过 Picard 标记后的重复序列。
- 低质量比对：MAPQ 小于 20 或 30 的 reads。
- 未比对上的 reads：SAM flag 4。
- 补充比对：SAM flag 2048。
- 比对到黑名单区域的 reads：如黑名单、线粒体、叶绿体、大肠杆菌等区域。

步骤：

使用 alignmentSieve 进行过滤：

alignmentSieve \
  --numberOfProcessors num_threads \
  --bam input_sorted.bam \
  --outFile output_filtered.bam \
  --filterMetrics log_directory/sieve_alignment.log \
  --ignoreDuplicates \
  --minMappingQuality min_quality \
  --samFlagExclude flag_value \
  --blackListFileName blacklist_file.bed \
  1>log_directory/sieve_alignment.log 2>&1

为过滤后的 BAM 文件构建索引：
```
samtools index output_filtered.bam
```

12. 估计过滤效果

工具：使用 DeepTools 包中的 estimateReadFiltering 工具估计每个步骤中过滤掉了多少 reads。

步骤：

运行 estimateReadFiltering：

estimateReadFiltering -p 1 \
-b raw_bams/OE_CHIP1.sorted.bam raw_bams/OE_CHIP2.sorted.bam \
--sampleLabels OE_CHIP1 OE_CHIP2 \
-o quality_control/estimate_read_filtering.tsv \
--ignoreDuplicates \
--minMappingQuality 25 \
--samFlagExclude 260
# --blackListFileName blacklist_file.bed \ 有黑名单可以选
  1>log_directory/estimate_read_filtering.log 2>&1

13. IGV 可视化与 BAM 文件转换为 BigWig 格式

目的：通过 IGV 对 ChIP-seq 数据进行可视化，发现对照组的噪声较大。为更好地处理和展示数据，将 BAM 文件转换为 BigWig 格式。

BAM文件是单碱基分辨率的比对文件，而bigwig文件通常是以固定窗口大小来表示信号强度的，这个窗口大小可以是任意的，比如可以是50个碱基对，但这并不是固定的，取决于生成bigwig文件时所使用的窗口大小设置。因此，bigwig文件更适合于查看整体的信号趋势，而不是单个碱基的精确比对信息。

BAM文件是单碱基分辨率的比对文件，而bigwig文件通常是以固定窗口大小来表示信号强度的，这个窗口大小可以是任意的这样有什么好处？
BigWig文件将BAM文件中的单碱基分辨率数据进行汇总，并以固定窗口大小（通常称为bin size或step size）来表示信号强度。这种方式带来了以下几个好处：

数据压缩与存储效率：
• BAM文件包含每个测序读段的精确位置信息，数据量庞大且不便于直接进行快速查询和可视化。而BigWig文件通过将信号强度在特定窗口内进行汇总，大大减少了数据量。这种压缩形式使得BigWig文件占用的存储空间更小，适合于在有限的存储资源下处理大规模数据集。
提高数据处理速度：
• 由于BigWig文件已经将信号强度进行预先汇总和归一化，因此在后续的数据处理和分析（如在基因组浏览器中进行可视化）时，可以避免逐个碱基读取原始数据的繁琐过程，从而显著提高了数据读取和显示的速度。
更灵活的信号展示：
• 固定窗口大小的信号强度表示允许研究者在不同的分辨率下观察数据。例如，在较大的窗口大小下可以快速获取整体趋势，在较小的窗口大小下可以更详细地观察局部变化。这种多尺度的信号展示方式在探索基因组区域时非常有用。
减少噪音影响：
• 将信号汇总到一定大小的窗口中，可以平均掉一些随机的噪音信号，使得整体趋势更为清晰。在一些分析场景中，特别是当数据本身噪音较大时，这种方式可以帮助更好地识别出真实的生物学信号。
增强对大规模数据集的可视化能力：
• 在浏览大规模基因组数据时，逐个碱基显示信息可能非常繁琐且不切实际。通过使用BigWig文件的固定窗口大小汇总，可以更好地适应常见的基因组浏览器（如UCSC Genome Browser）进行数据展示，方便研究者对全基因组的信号分布进行快速浏览和定位。
灵活的可调节性：
• 研究者可以根据分析的具体需求选择合适的窗口大小。例如，在探索全基因组范围内的趋势时，可以使用较大的窗口，而在分析局部热点区域时，可以使用较小的窗口，以便在不同的分辨率下获取更适合的信号强度信息。这种灵活性使得BigWig文件在各种生物信息学应用中都能很好地适应。

步骤：
1. 使用 bamCoverage 命令将 BAM 文件转换为 BigWig 文件：
```
bamCoverage -p num_threads \
  --bam input_bam.bam \
  -o output_bigwig.bw \
  --effectiveGenomeSize genome_size \
  --binSize bin_size \
  --normalizeUsing normalization_method \
  1>log_directory/convert_bam_to_bigwig.log 2>&1
```
2. 参数解释：
  - -p num_threads：指定用于处理的线程数量，提升处理速度。
  - --bam input_bam.bam：输入的 BAM 文件路径和名称。
  - -o output_bigwig.bw：输出的 BigWig 文件路径和名称，用于后续的可视化。
  - --effectiveGenomeSize genome_size：基因组的有效大小（以碱基对为单位），用于计算覆盖度。
  - --binSize bin_size：将 BAM 文件转换为 BigWig 文件时使用的窗口大小（以碱基对为单位），影响数据的分辨率。
  - --normalizeUsing normalization_method：归一化方法，例如 RPKM，用于标准化数据。
  - 1>log_directory/convert_bam_to_bigwig.log 2>&1：将标准输出和标准错误重定向到日志文件，便于调试和记录。

归一化一般都选PRKM
在将BAM文件转化为BigWig（BW）文件时，归一化通常是一个可选步骤，但它非常常见且重要，具体是否进行归一化取决于数据分析的需求和所使用的工具。归一化的目的是为了消除测序深度或其他技术性因素的差异，使得不同样本或实验条件之间的信号更具可比性。以下是一些情况和方法：

无归一化：
• 如果不进行归一化，BigWig文件中的信号强度会直接反映原始BAM文件中每个位置的读段数。这种情况下，不同样本之间的信号强度可能无法直接比较，因为它们的测序深度和覆盖度可能不同。
RPKM/FPKM（Reads/Fragments Per Kilobase of transcript, per Million mapped reads）：
• RPKM或FPKM归一化方法会考虑基因长度和测序深度，将读段数标准化为每千碱基和每百万读段上的读段数。这在基因表达水平的比较中非常常见。
CPM（Counts Per Million）：
• CPM归一化是一种相对简单的方法，将每个位置的读段数标准化为每百万总读段数。这种方法通常用于基因组范围内的比较，而不考虑基因长度。
RPM（Reads Per Million mapped reads）：
• RPM是CPM的一种形式，常用于ChIP-seq数据的归一化。它仅基于测序深度，将信号标准化为每百万测序读段中的读数。
其他归一化方法：
• 其他归一化方法如TPM（Transcripts Per Million）、DESeq2或EdgeR的归一化方法，可以针对特定类型的分析进行调整。

总结

以上内容涵盖了从 GTF 文件处理、检测基因的可变剪切、提取最长转录本，到测序数据的过滤和质控、比对到参考基因组、处理 BAM 文件，以及最终的可视化操作。这些步骤帮助你一步步优化数据质量，为下游分析奠定坚实基础。

如有进一步的学习内容或问题，欢迎继续分享。

14. IGV 使用bigwig格式文件绘制Reference-Point 图和Scaled Region Plot

在得到bw格式之后就可以进行Reference-Point 图的绘制

computeMatrix reference-point \
	-S $workdir/03.mapping/clean_bams/OE_CHIP1.bw \
	$workdir/03.mapping/clean_bams/OE_CHIP2.bw \
	--samplesLabel OE_CHIP1 OE_CHIP2 \
	-R $workdir/ref/genes_longest.bed6 \
	--referencePoint TSS \
	-b 2000 -a 2000 \
	--binSize 50 \
	-o $workdir/03.mapping/deeptools/tss_matrix_test.gz \
	-p 2

  -R <reference_bed1> \   # 参考基因组的bed文件，如最长基因
     <reference_bed2> \   # 参考基因组的bed文件，如打乱顺序的基因
  --referencePoint TSS \  # 指定参考点为转录起始位点（TSS）
  -b 5000 \               # 参考点上游区域距离（单位：bp）
  -a 5000 \               # 参考点下游区域距离（单位：bp）
  --binSize 50 \          # 区域内的bin大小（单位：bp）
  -o <output_file> \      # 输出文件路径，包含计算后的矩阵
  -p 1 > <log_file> 2>&1  # 使用一个线程，日志重定向到指定文件、

#画图
plotHeatmap \
	-m $workdir/03.mapping/deeptools/tss_matrix_test.gz \
 	-o $workdir/03.mapping/deeptools/tss_heatmap.pdf \
	--missingDataColor 0.5 \
	1>>$workdir/03.mapping/deeptools/plot_heatmap_reference_point.log 2>&1

在 ChIP-seq 数据分析中，绘制 Reference-Point 图的原理整理如下：

1. 选择参考点

定义参考点：常见的参考点包括转录起始位点（TSS）、转录因子结合位点、启动子区域或其他功能区域的中心。
目标区域：确定参考点后，通常会指定参考点上下游各 2 kb 的区域进行信号分析。

2. 信号提取

原始数据：ChIP-seq 实验产生的 BAM 文件包含 DNA 片段与参考基因组的比对信息。
提取信号：在参考点周围指定区域内提取信号强度，常用的指标包括覆盖度（coverage）、RPKM、CPM 等。

3. 信号对齐

对齐参考点：将所有目标区域的信号数据对齐到参考点，确保所有参考点在相同的位置。
对齐后的信号矩阵：形成一个矩阵，每一行代表一个参考点附近的信号分布，每一列代表参考点的一个特定位置。

4. 平均信号计算

平均化：计算所有目标区域在参考点附近的平均信号强度，消除单个基因或区域之间的噪音，展示全局趋势。
生成平均曲线：绘制平均信号强度在参考点周围的曲线图。

5. 可视化

曲线图：展示参考点周围平均信号强度的变化，直观地显示信号在参考点上下游的分布趋势。
热图：展示每个参考点在上下游区域的信号分布，通常按信号强度或排序后的基因排列。

6. 归一化和校正

信号归一化：使用 RPKM、CPM 或 TPM 等方法，去除测序深度和片段长度的影响，确保不同样本间的可比性。
背景校正：进行背景信号的校正，减少非特异性结合或实验噪音的影响。

这种图的可视化和分析有助于理解特定基因组位置的染色质状态或蛋白质结合情况，对于理解基因表达调控机制至关重要。

在chip-seq中Scaled Region Plot图的绘制原理

使用 computeMatrix scale-regions 来分析和转换信号数据到矩阵格式

命令如下
bash

computeMatrix scale-regions \
	-S $workdir/03.mapping/clean_bams/OE_CHIP1.bw \
	$workdir/03.mapping/clean_bams/OE_CHIP2.bw \
	--samplesLabel OE_CHIP1 OE_CHIP2 \
	-R $workdir/ref/genes_longest.bed6 \
	--regionBodyLength 4000 \
	-b 2000 -a 2000 \
	--binSize 50 \
	-o $workdir/03.mapping/deeptools/tss_tes_matrix.gz \
	-p 1 \
	1>$workdir/03.mapping/deeptools/plot_heatmap_scale_regions.log 2>&1

plotHeatmap \
  -m $workdir/03.mapping/deeptools/tss_tes_matrix.gz \
  -o $workdir/03.mapping/deeptools/tss_tes_heatmap.pdf \
  --missingDataColor 0.5 \
  1>>$workdir/03.mapping/deeptools/plot_heatmap_scale_regions.log 2>&1

ChIP-seq 数据分析中，Scaled Region Plot（缩放区域图）是为了展示在一组基因组区域（如一组基因、启动子区域或增强子区域）内部的信号分布模式。绘制 Scaled Region Plot 的原理是将不同长度的区域缩放到相同的长度，使得这些区域能够在同一尺度上进行比较和可视化。这种图通常用于分析这些区域的整体特征或寻找共同的信号模式。以下是 Scaled Region Plot 的绘制原理：

选择目标区域
• 定义区域：首先，需要选择一组目标区域，这些区域可以是基因、启动子、增强子、或者任何感兴趣的基因组区域。每个区域的长度可能不同。
• 区域类型：这些目标区域的类型会影响图的解读。例如，选择全基因范围的区域可能用于展示基因体内的修饰分布，而选择启动子区域则可能用于展示转录起始位点附近的信号。
区域缩放
• 统一区域长度：由于这些区域的长度可能不一致，在绘制 Scaled Region Plot 时，必须将每个区域缩放到相同的长度。例如，如果选择的是基因体区域，所有基因将被缩放到相同的基因体长度（通常通过线性插值实现）。
• 区域分段：每个缩放后的区域通常会被划分为若干等长的分段（bins），每个分段表示该区域的一个相对位置。这样，每个区域的信号就可以被比较。
信号提取
• 计算信号强度：在每个缩放后的区域中，计算各个分段的 ChIP-seq 信号强度，通常表示为覆盖度（coverage）或其他归一化的信号值。这样，尽管原始区域的长度不同，信号在每个缩放后区域内的分布都可以被直接比较。
• 平均信号：对所有目标区域的信号强度进行平均，得到每个分段的平均信号值。这些平均信号值将用于生成最终的 Scaled Region Plot。
可视化
• 生成图表：绘制一条展示所有区域平均信号的曲线，横轴表示区域的相对位置（从区域的起点到终点），纵轴表示信号强度。通过这种方式，可以观察到信号在这些区域内部的分布模式。
• 热图展示：除了信号曲线，通常还会生成一个热图，展示每个具体区域在缩放后的信号分布情况。每行表示一个区域，每列表示该区域内的一个分段。热图能够展示不同区域之间的信号分布差异，以及在这些区域内部的信号变化模式。
实例
举个简单的例子，如果我们有以下两个基因：
• 基因 A：长度为 2000 bp
• 基因 B：长度为 1000 bp
我们希望将它们缩放到 10 个相等的分段。
• 对于基因 A，每个分段将覆盖 200 bp（2000 bp / 10 分段）。
• 对于基因 B，每个分段将覆盖 100 bp（1000 bp / 10 分段）。
然后，对每个基因的每个分段，计算该分段内的 ChIP-seq 信号（例如，覆盖度的平均值），将这些信号归一化到相同的分段数。

生成 genes_longest.bed6 文件

目的：genes_longest.bed6 文件包含了可能是该物种中基因长度最长的集合。这些信息用于定位感兴趣的区域（如基因体），并对 ChIP-seq 信号在转录起始位点（TSS）上下游的分布进行分析。
方法：从 GTF 文件中提取最长的基因。使用如下 AWK 命令筛选正链上的“gene”特征，并生成 BED 文件。
bash
awk ‘$3==“gene” && $7==“+”’ genes.gtf > genes_longest.gtf
然后，你需要将 GTF 转换为 BED 格式，可能还需要进一步处理以确保选择每个基因的最长转录本。

生成 genes_shuffle.bed6 文件

目的：genes_shuffle.bed6 文件包含一组随机化的基因或区域，作为对照使用。这有助于评估信号在特定区域的特异性富集，或是否也在基因组的其他随机区域出现。
实现：使用 bedtools shuffle 命令重新定位 genes_longest.bed6 文件中的区域，随机地将其分配到基因组的其他位置，同时保持区域的长度和染色体分布不变，但排除已知的敏感区域。
bash
bedtools shuffle -i genes_longest.bed6 -g genome.fa.fai -o genes_shuffle.bed6
这个工具将输入文件中的区域随机地重新定位，每个区域的染色体位置、起始和终止位置都可能改变，而区域的名称和长度保持不变。
为什么使用 shuffle 而不是 random？
• bedtools shuffle：用于将已存在的基因组区域重新随机定位到基因组的其他位置，适合当需要保持区域长度和染色体分布时。这对于评估特定基因组事件（如结合位点、甲基化位点）的生物学意义是有帮助的。
• bedtools random：生成完全随机的基因组区域，适用于需要完全随机位置和长度的场合，常用于生成模拟数据或作为统计分析的背景数据集。
这两种工具各有用途，选择使用哪一个取决于你的具体需求和研究目标。

``bash

				第四步

                 |

   ------------------------------------------

   |         	   |                           |

 Peak 过滤    多样本 Peak 合并               提取序列
   |            |                           |
cutoff analysis   bedtools merge           narrowPeak
   |            |                           |
intersect         narrowPeak          --> 提取 summit 附近序列
   |            |                           |
IDR --> narrow Peak, 2个生物学重复    broadPeak
                |                           |
       --> 提取 summit 附近序列      --> 提取整个 peak 序列

第15步：峰值检测 call peak

在Scaled Region Plot和Reference-Point生成后，接下来进行峰值检测，通常使用MACS3软件进行。
关键步骤：估计d值
• 双末端测序（PE）：
o 使用参数 -f BAMPE，MACS3会自动统计d值。

• 单末端测序（SE）：
o 使用参数 -f BAM，让MACS3自动估计d值。
o 或者，可以手动指定d值，使用参数 -nomodel -extsize d -shift。

峰的类型选择
• Narrow Peak：
o 通常用于转录因子（TF）分析。
• Broad Peak：
o 当两个窄峰连在一起时，可能被计算为Broad Peak。此时需要设置参数 -broad。
MACS3常用命令

macs3 callpeak \
  -g 2178268108 \ #有效基因组大小
  -t $workdir/03.mapping/OE_CHIP1.sorted.bam \
  -c $workdir/03.mapping/OE_CHIP2.sorted.bam \
  --name $workdir/04.macs3_call_peaks/t_1_c_2/OE_CHIP \
  --format BAMPE \
  --keep-dup all \ #
  --qvalue 0.075 \
  --cutoff-analysis \
  1>$workdir/04.macs3_call_peaks/t_1_c_2/OE_CHIP_callpeak_with_macs3.log 2>&1

生成的文件
• _peaks.narrowPeak 或 _peaks.broadPeak 文件：
o 这些是根据分析类型（窄峰或广峰）生成的，符合BED格式的标准。
o 文件包含峰的位置和额外的统计信息（如信号强度、p-value 和 q-value）。
o 这些文件适用于在基因组浏览器中的可视化或作为后续分析步骤的输入。
后续步骤
• 可以将 callpeak 结果在IGV（Integrative Genomics Viewer）中进行显示，以可视化分析结果。

峰值过滤和筛选

1.对于无生物重复的

对于没有生物学重复的样本，首先需要根据之前的cutoff analysis结果来确定合适的qscore值：

绘制图表：
o X轴：qscore
o Y轴：npeaks（峰的数量）
确定拐点：
o 如果图中有明显的拐点，可以选择拐点处的qscore值作为阈值。
o 如果没有明显的拐点，可以根据保留多少峰（npeaks）的需求来确定qscore值。
选择Fold Enrichment阈值
• 选择一个合适的Fold Enrichment阈值，用于进一步过滤峰值数据。
• Fold Enrichment表示信号与背景噪音的比率，是一个重要的筛选标准。
使用awk命令进行数据过滤
为了根据特定阈值过滤峰值数据，可以使用awk命令。该命令允许选择符合特定条件的行，例如选择某列值大于或小于指定阈值的行。

awk命令示例：
awk '$9>3' \  #第9列为q值，3根据上面qscore图确定
  xxxxxxxxxxxxxx.narrowPeak > ......
#我筛选后的文件
$ awk '{print $1}' OE_CHIP_peaks.narrowPeak |sort -V |uniq
chr1
chr2
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chr10
chrscaf_22
chrscaf_23
chrscaf_27
chrscaf_28
.......
#仅保留chr1:chr10
awk -v OFS='\t' 'BEGIN {FS=OFS} \ $1 ~ /^chr[1-9]|^chr10/ {print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10}' \
	OE_CHIP_peaks.narrowPeak >noscaf_OE_CHIP_peaks.narrowPeak
1.#-v OFS='\t'：设置输出字段分隔符为制表符。
2.#BEGIN {FS=OFS}：在处理任何输入行之前，设置输入字段分隔符与输出字段分隔符相同。
3.#$1 ~ /^chr1-9|^chr10/：
#这是一个正则表达式匹配，用于检查第一列（$1）是否符合以下条件之一：
#以 "chr1" 到 "chr9" 开头，后面紧跟冒号或行尾（表示 "chr1" 到 "chr9" 本身）。
#正好是 "chr10"。
4.#~ 是 awk 中用于模式匹配的操作符。
5.#print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10：如果第一列符合上述条件，则打印所有十个字段。
6.#OE_CHIP_peaks.narrowPeak：指定输入文件。
7.#| sort -k1,1：将 awk 命令的输出通过管道传递给 sort 命令，按照第一列进行排序。
 #k1,1表示从1列到1列，n表示按数值排列
8.#| uniq -w 1：将排序后的输出传递给 uniq 命令，-w 1 选项指定只比较第一列进行去重

#检测一下
awk '{print $1}' noscaf_OE_CHIP_peaks.narrowPeak |sort -V |uniq
chr1
chr2
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chr10
 
 #改一下第4列内容
 awk -v OFS='\t'  \
 'BEGIN {FS=OFS} $4=$1":"$2"_"$3; \ {print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10}' \
 noscaf_OE_CHIP_peaks.narrowPeak >1

2.处理有生物学重复的峰值
步骤1：查看每个重复的峰值数量
• 对于有生物学重复的实验，首先需要查看每个重复样本中的峰值数量。
• 这一步有助于理解每个样本的峰值分布情况，以便进行后续分析。
步骤2：求交集
• 为了确保生物学重复的可靠性，通常会对每个重复样本的峰值进行交集分析。
• 交集分析可以使用bedtools intersect工具来实现。建议在求交集前，先使用cutoff过滤掉不符合条件的峰值，以提高结果的精确度。
bedtools intersect命令示例
bedtools intersect
-f <overlap_fraction> \ # 重叠百分比，指定两个区域至少需要重叠指定的百分比才认为它们相交
-r \ # 相互重叠，要求A和B文件中的区域互相重叠指定的百分比
-a <file_A> \ # A文件路径，指定第一个要比较的文件
-b <file_B> \ # B文件路径，指定第二个要比较的文件
-wa \ # 输出 A 文件的重叠部分

<output_file> # 输出文件路径，将结果保存到新的文件
示例说明：
• <overlap_fraction>：指定两个区域至少需要重叠的百分比。例如，-f 0.5表示区域需要有50%的重叠才算相交。
• -r：要求A和B文件中的区域互相重叠指定的百分比，这有助于提高交集的精确度。
• <file_A>：指定第一个要比较的文件路径。
• <file_B>：指定第二个要比较的文件路径。
• -wa：输出A文件中与B文件重叠的部分。
• <output_file>：指定输出文件的路径，将交集结果保存到这个文件中。
选择A文件的考虑因素
• 峰值数量：通常选择峰值数量较多的文件作为A文件。
• 测序质量：如果某个样本的测序质量更高，则可优先选择该样本作为A文件。
• 比对率：比对率高的样本可能更具有代表性，可以考虑作为A文件。
IDR分析
为什么使用IDR分析？
传统的交集方法只考虑了峰值的位置，而未考虑峰值的强度。IDR分析能够有效地将生物学重复中的峰值强度纳入考量，从而筛选出更加可靠的峰值。
步骤1：对重复样本的峰值进行排序
• 在进行IDR分析之前，需要将每个生物学重复样本中的峰值按照qscore进行排序。
• qscore越高，峰值越显著。
步骤2：使用IDR进行分析
IDR分析结合了多个生物学重复中的峰值强度，并根据预设的阈值筛选出具有高重复性的峰值。
IDR命令示例：
idr
–samples \ # 指定输入文件，即用于IDR分析的样本文件列表
<sample_file_1> \ # 第一个样本文件路径
<sample_file_2> \ # 第二个样本文件路径
–idr-threshold <idr_threshold> \ # 指定IDR阈值，一般为0.05
–output-file <output_file_path> \ # 输出文件路径，将结果保存到指定文件中
–plot \ # 生成IDR分析的图表
–input-file-type <file_type> \ # 指定输入文件的类型
–rank <ranking_metric> \ # 指定用于排序的指标
1><log_file_path> \ # 标准输出重定向到日志文件路径
2>&1 # 标准错误重定向到与标准输出相同的文件
示例说明：
• <idr_threshold>：指定IDR的阈值，通常设置为0.05。
• <ranking_metric>：指定用于排序的指标，通常为qscore。
• <file_type>：输入文件的类型，如narrowPeak格式。
• <sample_file_1>和<sample_file_2>：分别是重复样本的峰值文件路径。
• <output_file_path>：指定输出文件的路径。
步骤3：生成的IDR文件解释
IDR分析的输出文件包括以下重要信息：
• 前1-10列：包含IDR合并后的Peaks信息，以narrowPeak格式表示。
• 第11列及以后：包含原始各个重复样本的Peaks信息。
• IDR score：使用公式min(int(log2(-125*IDR), 1000))计算的得分。
• localIDR：局部IDR值，以-log10表示。
• globalIDR：全局IDR值，以-log10表示。
步骤4：提取通过IDR阈值的峰值并转换为narrowPeak格式
可以使用awk命令提取通过IDR分析阈值的峰值，并将其转换为narrowPeak格式。
awk命令示例：
awk
-v OFS=“\t” \ # 设置输出字段分隔符为制表符
‘BEGIN {FS=OFS} \ # 设置输入和输出的字段分隔符为制表符
{ $4=$1":“$2”_"$3; print $1, $2, $3, $4, $5, $6, $7, $8, $9, $10 }’ \ # 将第4列修改为染色体名、起始位置和结束位置的组合，并打印指定的列
<input_file> | \ # 指定输入文件路径
sort -k1,1 -k2,2n -k3,3n \ # 按照第1列（染色体名）、第2列（起始位置）和第3列（结束位置）进行排序
<output_file> # 输出排序后的结果到指定文件
示例说明：
• <input_file>：IDR分析的输出文件路径。
• <output_file>：转换后的narrowPeak格式文件的保存路径。

16Motif分析

在这里插入图片描述

#对于broad peak
#按q值提取前1000条（最多5000-6000条）进行motif分析
sort -k8,8nr 1 |head -n 1000 \
|awk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6}' >broad_1000_reads_motif
#提取序列
bedtools getfasta 	\
-fi ~/yiyaoran/chip-seq/mychip-seq/ref/genome.fa \
-bed broad_1000.bed6 >broad_1000_reads_motif


#对于narrow peak
sort -k8,8nr made_noscaf_OE_CHIP_peaks.narrowPeak \
|head -n 1000|awk 'BEGIN {OFS="\t"} \ {print $1,$2+$10-100,$2+$10+100,$4,$5,$6}' - >narrow_1000.bed6
#提取序列
bedtools getfasta 	\
-fi ~/yiyaoran/chip-seq/mychip-seq/ref/genome.fa \
-bed narrow_1000.bed6 >narrow_1000_reads_motif

Motif 分析中sequence logo的基本概念
• Motif 长度：横坐标通常在 6-12 bp 之间。
• 信息量 (Bits)：纵坐标表示信息量范围为 0-2 bits。
• Narrow Peak 分析：提取 summit 两侧 100 bp 或 50 bp 进行 motif 分析。
• Broad Peak 分析：提取整个 peak 序列进行 motif 分析。

一、broad peak

IDR 分析后提取前 1000 个 peaks
• 排序并提取前 1000 个 peak：
bash
sort -k<列号>,<列号><排序方式> <输入文件> | head -1000 > top1000_peaks.bed
例如，如果你想根据第五列进行数值降序排序并提取前 1000 条：
bash
sort -k5,5nr <输入文件> | head -1000 > top1000_peaks.bed
• 生成 BED6 文件（可选）：
bash
sort -k5,5nr <输入文件> | head -1000 | awk ‘BEGIN{OFS=“\t”}{print $1,$2,$3,$4,$5,$6}’ > head1000_peaks.bed
提取序列信息
• 使用 bedtools getfasta 从基因组中提取序列：
bash
bedtools getfasta -fi <reference.fa> -bed top1000_peaks.bed -fo output.fa
o -fi <reference.fa>：参考基因组文件。
o -bed <regions.bed>：包含你要提取的坐标的 BED 文件。
o -fo <output.fa>：输出的 FASTA 文件。

二、Narrow Peak 处理步骤
1.注释
2.搜库
3.搜基因组

Motif 注释与分析

Motif Discovery（motif 发现）：
o 可以使用 MEME 或 HOMER 进行 motif 发现分析。
o MEME 适合处理长度为 1-49 bp 的输入文件，而 STREME 适合处理 >50 bp 的序列。
o 在 MEME 的网页端进行操作时，选择 motif discovery，并在 input sequences 中上传你的 FASTA 文件。

操作如下