ChIP-seq数据分析实例（针对长染色体无法建立index导致bw文件生成失败）

写在前面: 本流程以染色体较长的黑麦基因组为例，每条染色体达到1Gb左右。

Step1：对下机数据进行质控，去接头

fastp -i 20P13_CENH3_ChIp_R1.fq.gz -I 20P13_CENH3_ChIp_R2.fq.gz -o 20P13_CENH3_ChIp_1.clean.fq.gz -O 20P13_CENH3_ChIp_2.clean.fq.gz -z 4 -q 20 -u 30 -n 10 -h 20P13_CENH3_ChIp.clean.html

Step2：利用bwa将测序数据比对到参考基因组中

Note：由于黑麦基因组的染色体ID为GWHASIY00000001-8，为方便后续对bam文件排序，因此可以将其染色体ID转换位Chr1-Chr8。

更改黑麦基因组染色体之前的ID：

image

使用脚本对黑麦基因组染色体ID进行批量更改

for i in {1..8}; do sed -i "s/GWHASIY0000000$i/Chr$i/g" GWHASIY00000000.genome.fa;done

image

对参考基因组建立索引文件

bwa index -a bwtsw  GWHASIY00000000.genome.fa

（适合于大于2G的参考基因组建立索引）

利用bwa mem进行比对

bwa mem -M -t 20 GWHASIY00000000.genome.fa 20P13_CENH3_ChIp_1.clean.fq.gz 20P13_CENH3_ChIp_2.clean.fq.gz >20P13_CENH3_ChIp.sam 2> 20P13_CENH3_ChIp_bwa.log

Tip：-t 后接线程数，根据服务器情况自行设定线程数，GWHASIY00000000.genome.fa为基因组文件，且基因组目录下需要放置对于基因组文件的索引文件

Step3：对sam文件进行排序

Note：通常情况下，需要将sam先转化为bam再进行排序，但是samtools对于长染色体基因组排序存在无法正常排序的问题。下面就是用samtools排序后的bam文件，对samtools排序后的bam文件的第三列（染色体ID列）进行查看，发现并没有正确排序；

samtools view PI428373_CENH3_ChIp.sorted.bam |cut -f 3 | uniq

如果排序正常输出的应该是Chr1-8，但是用上述脚本检查，发现并没有正确排序，这也是无法对该bam文件建立索引的原因。

因此我们需要使用perl脚本对sam文件进行排序（先排序再转bam）

perl sort_sam.pl 20P13_CENH3_ChIp.sam
$cat sort_sam.pl

#!/usr/bin/perl
use warnings;
use strict;
open IN,"$ARGV[0]";
my ($name) = $ARGV[0] =~ /(.*).sam/;
open OUT,">$name.perl.sorted.sam";
my @data = ();
while (<IN>) {
    my @a = split/\s+/;
    next unless $a[2] =~ /^Chr/;
    push @data,[$a[2],$a[3],$_];
}
close IN;
foreach  (sort{$a->[0] cmp $b->[0] || $a->[1] <=> $b->[1]}@data) {
    print OUT "$_->[2]";
}
close OUT;

生成20P13_CENH3_ChIp.perl.sorted.sam文件（需要用cat命令为sam文件加头文件）

Step4：将排序后的sam文件转化为bam文件

samtools view -Sb -o 20P13_CENH3_ChIp.perl.sort.bam 20P13_CENH3_ChIp.perl.sort.sam

Step5：去PCR重复并建立索引

samtools rmdup 20P13_CENH3_ChIp.perl.sort.bam 20P13_CENH3_ChIp.perl.sort_rmdup.bam

samtools index -c 20P13_CENH3_ChIp.perl.sort_rmdup.bam

Step6：bam文件转bw（bigwig格式）文件

bw格式详解：BW格式详解

bamCoverage --normalizeUsing RPGC --effectiveGenomeSize 7730188902 --binSize 5 --bam  20P13_CENH3_ChIp.perl.sort_rmdup.bam -o  20P13_CENH3_ChIp.perl.sort_rmdup.bw

其中effectiveGenomeSize后接基因组大小，可以用下面的python脚本进行计算

python GenomeSize.py GWHASIY00000000.genome.fa
$cat GenomeSize.py

#coding=utf-8
import sys
aList=[]
fa_file = sys.argv[1]
with open(fa_file,'r') as f:
    for line in f:
        line = line.strip()
        line = line.upper()
        if not line.startswith(">"):
            baseA = line.count("A")
            baseT = line.count("T")
            baseC = line.count("C")
            baseG = line.count("G")
            aList.extend([baseA, baseT, baseC, baseG])
#            print(aList)
print("effective_genome_size =", sum(aList))

Step7：IGV可视化（Windows本地可视）

将基因组文件以及索引文件，bw文件放在同一文件夹，打开IGV，Genomes-Load Genomes from Files导入基因组文件，File-Load form file导入bw文件即可进行可视化。

Tip：IGV下载地址Downloads | Integrative Genomics Viewer。如果电脑没有安装java（IGV必须依赖JAVA），根据电脑系统选择Java include下载安装即可。

Note：在前面Sam排序转bam文件都可以通过管道一步完成，但是对于服务器性能不强的建议分步完成，避免中间文件出错。以上软件都可以用conda进行安装。