ChIP-seq上游分析流程

0. 数据格式及数据存放文件夹说明

获得了测序公司提供的fq.gz文件，每个样本均为双端测序，得到的数据分别命名为“样品名1.fq.gz”与“样品名2.fq.gz”。

下设文件夹包含主要文件说明如下：

/raw 存放公司发回的原始数据。

/clean 存放使用trim_galore软件过滤后的数据。

/qc_raw与/qc_clean 分别对/raw与/clean中的数据进行质控。

/align 使用bowtie2软件将/clean中的序列与hg19参考基因组进行比对，并排序得到的bam文件。

/rmdup 使用samtools markdup -r进行比对结果bam文件的PCR重复序列去除。

/callpeaks 使用macs2软件对去除重复后的bam文件进行callpeaks

1. 使用trim_galore进行过滤

trim_galore.bash脚本位于/raw文件夹中，脚本配置如下。

ls *1.fq.gz | while read id
do
trim_galore -q 25 --phred33 --length 25 -e 0.1 --stringency 4 -j 32 --paired $id $(basename $id "1.fq.gz")2.fq.gz -o ../clean
done

2.使用fastqc与multiqc对/raw与/clean中的数据进行质控

qc_raw.bash脚本位于/qc_raw文件夹中，qc_clean.bash脚本位于/qc_clean文件夹中，配置分别如下。

ls ../raw/*gz | xargs fastqc -t 32 -o ./
multiqc ./

ls ../clean/*gz | xargs fastqc -t 32 -o ./
multiqc ./

3.使用bowtie2软件对/clean中的数据进行比对，并使用samtools对比对得到的bam文件进行排序

分为这么几步：首先对双端测序的两个文件共同做bowtie2比对，得到sam文件。再将sam文件转bam后进行排序，并对排序后的bam文件进行索引建立与质控。脚本位于/align文件夹，配置如下：

ls ../clean/*1.fq.gz | cut -d "_" -f 1 | while read id;
do
        bowtie2 -x /data/ljouyang/reference/index/bowtie2/hg19/hg19 -p 56 -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -S ${id}.sam
done
mv ../clean/*sam ./
ls *sam | while read id
do
        samtools view -bhS -q 20 $id > $(basename $id ".sam").bam
        samtools sort $(basename $id ".sam").bam -@ 32 -o $(basename $id ".sam").sorted.bam
        samtools index -@ 32 $(basename $id ".sam").sorted.bam
        samtools flagstat $(basename $id ".sam").sorted.bam > $(basename $id ".sam").stat
done

cat *stat* | grep % > qc_result.txt
rm *sam

4. 使用samtools中的markdup -r功能对排序后的bam文件进行去除PCR重复

若直接对按染色体排序的bam文件进行PCR重复的去除的话，会提示缺少MC tags，因此需要先转为按照reads排序的bam文件，再添加MC tags，再转回按染色体排序的bam文件，最后去除PCR重复。脚本位于rmdup文件夹，配置如下：

cp ../align/*sorted.bam ./
ls  *.sorted.bam  | while read id
do
samtools sort $id -n -@ 32 -O bam -o $(basename $id ".sorted.bam").nsorted.bam;
samtools fixmate -m -@ 32 -O bam $(basename $id ".sorted.bam").nsorted.bam $(basename $id ".sorted.bam").nmated.bam;
samtools sort $(basename $id ".sorted.bam").nmated.bam -@ 32 -O bam -o $(basename $id ".sorted.bam").allsorted.bam;
samtools markdup -r -@ 32 $(basename $id ".sorted.bam").allsorted.bam $(basename $id ".sorted.bam").rmdup.bam ;
done

ls  *.rmdup.bam  | while read id
do
samtools flagstat $id > $(basename $id ".bam").stat;
done

cat *stat* | grep % > qc_result.txt
rm *.sorted.bam
rm *.nsorted.bam
rm *.nmated.bam
rm *.allsorted.bam

5. 使用macs2对去除PCR重复后的bam文件进行callpeaks

不同批次的测序数据的实验处理不同，无法进行批处理。仅以第三次测序数据为例，脚本如下：

macs2 callpeak -c Hela2-BKDL190828263.rmdup.bam -t Hela1-BKDL190828262.rmdup.bam -f BAM -g hs -n Hela1_rmdup
macs2 callpeak -c U2OS4-BKDL190828265.rmdup.bam -t U2OS3-BKDL190828264.rmdup.bam -f BAM -g hs -n U2OS3_rmdup
macs2 callpeak -c U2OS4-BKDL190828265.rmdup.bam -t U2OS5-BKDL190828266.rmdup.bam -f BAM -g hs -n U2OS5_rmdup
macs2 callpeak -c U2OS4-BKDL190828265.rmdup.bam -t U2OS6-BKDL190828267.rmdup.bam -f BAM -g hs -n U2OS6_rmdup

最后友情宣传生信技能树

生物信息学“义诊”
生物信息学"拍卖会"
全国巡讲:R基础，Linux基础和RNA-seq实战演练 : 预告：12月28-30长沙站
广州珠江新城GEO数据挖掘滚动开班
DNA及RNA甲基化数据分析与课题设计

最后编辑于：2019.11.22 03:00:01