操作记录-2020-11-13：精简代码处理ChIP_seq数据

今天准备尝试编写一组精简代码用于处理ChIP_seq数据，希望能成功吧。

1.建立相应目录

对新数据建立对应实验人员（lizexing）、测序类型（ChIP_seq）和日期（2020_11_13）的目录。

# 建立后如下：
(base) zexing@DNA:~/projects/lizexing/ChIP_seq/2020_11_13$

# 新建对应的目录
mkdir raw_data clean_data bam bam_bw bam_sort sam macs2_bdgdiff macs2_callpeak matrix_reference_point matrix_scale_regions fastqc_report MD5_txt scripts_log

2.检查数据完整性

cat md5.txt > check_md5sum.txt && md5sum -c check_md5sum.txt

操纵记录如下

(base) zexing@DNA:~/projects/lizexing/ChIP_seq/2020_11_13/clean_data$ cat *.txt > check_md5sum.txt && md5sum -c check_md5sum.txt
Scr-S-L-KQ_FKDL202604880-1a_1.clean.fq.gz: OK
Scr-S-L-KQ_FKDL202604880-1a_2.clean.fq.gz: OK
shTgm1-2-S-L-KQ_FKDL202604881-1a_1.clean.fq.gz: OK
shTgm1-2-S-L-KQ_FKDL202604881-1a_2.clean.fq.gz: OK
shTgm2-1-S-L-KQ_FKDL202604882-1a_1.clean.fq.gz: OK
shTgm2-1-S-L-KQ_FKDL202604882-1a_2.clean.fq.gz: OK

3.根据需要精简文件名称

操作记录如下：

(base) zexing@DNA:~/projects/lizexing/ChIP_seq/2020_11_13/clean_data$ ll
total 2.7G
drwxrwxr-x  2 zexing zexing 4.0K 11月 16 11:45 .
drwxrwxr-x 15 zexing zexing 4.0K 11月 13 10:36 ..
-rw-rw-r--  1 zexing zexing  476 11月 16 11:31 check_md5sum.txt
-rw-rw-r--  1 zexing zexing  152 11月 16 10:50 MD5_Scr-S-L-KQ_FKDL202604880-1a.txt
-rw-rw-r--  1 zexing zexing  162 11月 16 11:03 MD5_shTgm1-2-S-L-KQ_FKDL202604881-1a.txt
-rw-rw-r--  1 zexing zexing  162 11月 16 11:21 MD5_shTgm2-1-S-L-KQ_FKDL202604882-1a.txt
-rw-rw-r--  1 zexing zexing 401M 11月 16 10:57 Scr_1.clean.fq.gz
-rw-rw-r--  1 zexing zexing 410M 11月 16 11:01 Scr_2.clean.fq.gz
-rw-rw-r--  1 zexing zexing 495M 11月 16 11:11 shTgm1-2_1.clean.fq.gz
-rw-rw-r--  1 zexing zexing 510M 11月 16 11:20 shTgm1-2_2.clean.fq.gz
-rw-rw-r--  1 zexing zexing 423M 11月 16 11:26 shTgm2-1_1.clean.fq.gz
-rw-rw-r--  1 zexing zexing 437M 11月 16 11:27 shTgm2-1_2.clean.fq.gz

4. 在Linux服务器中对ChIP_seq数据进行处理并常规call peak。

vim新建ChIP_seq_script_1将数据质控、比对、格式转换、排序、生成目录、bamCoverage命令转换文件格式和macs2 callpeak综合在一起。

#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program
#     This program is used for ChIP-seq data analysis.
# History
#     2020/11/13       zexing            First release
# 设置变量${dir}为常用目录
dir=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13 # 用户名称和日期需要更改

# 对数据进行质控
fastqc -t 16 -o ${dir}/fastqc_report/ ${dir}/clean_data/*.fq.gz

# 利用for循环进行后续操作
for i in Scr shTgm1-2 shTgm2-1 # 样品名称需要修改
do
# 对数据进行比对
bowtie2 -t -p 16 -x /f/xudonglab/zexing/reference/UCSC_mm10/bowtie2_index/mm10 -1 ${dir}/clean_data/${i}_1.clean.fq.gz -2 ${dir}/clean_data/${i}_2.clean.fq.gz -S ${dir}/sam/${i}.sam

# 对数据进行格式转换
samtools view -@ 16 -S ${dir}/sam/${i}.sam -1b -o ${dir}/bam/${i}.bam

# 对数据进行排序
samtools sort -@ 16 -l 5 -o ${dir}/bam_sort/${i}.bam.sort ${dir}/bam/${i}.bam

# 对数据生成目录
samtools index -@ 16 ${dir}/bam_sort/${i}.bam.sort 

# bamCoverage命令转换文件格式
bamCoverage -p 16 -v -b ${dir}/bam_sort/${i}.bam.sort -o ${dir}/bam_bw/${i}.bam.sort.bw

# 使用macs2进行常规callpeak
macs2 callpeak -t ${dir}/bam_sort/${i}.bam.sort -f BAM -g mm -B -q 0.05 --outdir ${dir}/macs2_callpeak/ -n ${i}

done

在后台运行ChIP_seq_script_1：

nohup bash ChIP_seq_script_1 > ChIP_seq_script_1_log &

5. 使用deeptools软件绘制热图/密度图

1. computeMatrix scale-regions模式计算信号强度并用plotHeatmap/plotProfile作图

scale-regions模式计算的是区域形式，所以指定作图位置的BED或GTF格式文件为macs2 callpeak生成的后缀为peaks.narrowPeak的BED文件。

vim新建ChIP_seq_script_2，脚本如下：

#! /bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program:
#       This program is used for computeMatrix scale-regions.
#History:
# 2020/11/03         zexing              First release
# In the scale-regions mode, all regions in the BED file are stretched or shrunken to the length (in bases) indicated by the user.
# 参数-R 指定作图位置的BED或GTF格式文件，可用#标记同一组区域，默认无。
# 参数-S 输入bigwig文件。
# 参数-o 指定输出为文件名用于plotHeatmap, plotProfile
# 参数-b上游（默认0bp），-a下游（默认0bp）设定感兴趣的区域，如果该区域是基因，则为基因TSS上游或TES下游。
# 参数--skipZeros设定0分区域的处理
# 参数-p 设置线程数

bed=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13/macs2_callpeak
bw=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13/bam_bw
results=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13/matrix_scale_regions

computeMatrix scale-regions \
-R ${bed}/Scr_peaks.narrowPeak ${bed}/shTgm1-2_peaks.narrowPeak ${bed}/shTgm2-1_peaks.narrowPeak \
-S ${bw}/Scr.bam.sort.bw ${bw}/shTgm1-2.bam.sort.bw ${bw}/shTgm2-1.bam.sort.bw \
-o ${results}/matrix_scale_threegroups.gz \
-b 1000 -a 1000 \
-p 16

# 使用plotHeatmap对结果绘制热图并聚类
plotHeatmap -m ${results}/matrix_scale_threegroups.gz \
-o ${results}/scale_threegroups_heatmap.png \
--dpi 750 \
--whatToShow 'heatmap and colorbar' \
--startLabel "Start" \
--endLabel "End" \
--regionsLabel Scr-peak shTgm1-2-peak shTgm2-1-peak \
--samplesLabel Scr shTgm1-2 shTgm2-1
# 使用plotProfile对结果绘制密度图并聚类
plotProfile -m ${results}/matrix_scale_threegroups.gz \
-o ${results}/scale_threegroups_profile.png \
--dpi 750 \
--legendLocation upper-right \
--startLabel "Start" \
--endLabel "End" \
--regionsLabel Scr-peak shTgm1-2-peak shTgm2-1-peak \
--samplesLabel Scr shTgm1-2 shTgm2-1 \
--perGroup

后台运行ChIP_seq_script_2脚本如下

nohup bash ChIP_seq_script_2 > ChIP_seq_script_2_log &

2. computeMatrix reference-point模式计算信号强度并用plotHeatmap/plotProfile作图

reference-point模式计算的是峰值高点模式，所以指定作图位置的BED或GTF格式文件为macs2 callpeak生成的后缀为summits.bed的BED文件。

vim新建ChIP_seq_script_3，脚本如下：

#! /bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program:
#       This program is used for computeMatrix reference_point.
#History:
# 2020/11/03         zexing              First release
# Reference-point refers to a position within a BED region(e.g., the starting point). In this mode, only those genomicpositions before (upstream) and/or after(downstream) of the reference point will be plotted.
# 参数-R 指定作图位置的BED或GTF格式文件，可用#标记同一组区域，默认无。
# 参数-S 输入bigwig文件。
# 参数-o 指定输出为文件名用于plotHeatmap, plotProfile
# 参数-b上游（默认0bp），-a下游（默认0bp）设定感兴趣的区域，如果该区域是基因，则为基因TSS上游或TES下游。
# 参数--skipZeros设定0分区域的处理
# 参数-p 设置线程数

dir=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13
bed=${dir}/macs2_callpeak
bw=${dir}/bam_bw
results=${dir}/matrix_reference_point
computeMatrix reference-point \
-R ${bed}/Scr_summits.bed ${bed}/shTgm1-2_summits.bed ${bed}/shTgm2-1_summits.bed \
-S ${bw}/Scr.bam.sort.bw ${bw}/shTgm1-2.bam.sort.bw ${bw}/shTgm2-1.bam.sort.bw \
-o ${results}/matrix_reference_threegroups.gz \
-b 1000 -a 1000 \
-p 16

# 使用plotHeatmap对结果绘制热图并聚类
plotHeatmap -m ${results}/matrix_reference_threegroups.gz \
-o ${results}/reference_threegroups_heatmap.png \
--dpi 750 \
--whatToShow 'heatmap and colorbar' \
--refPointLabel "Peak_center" \
--regionsLabel Scr-peak shTgm1-2-peak shTgm2-1-peak \
--samplesLabel Scr shTgm1-2 shTgm2-1