我的第一次ChIP-seq实践

最近按照Jimmy前辈分享的ChIP-seq实战分析（点击阅读原文），从零开始动手实践了一次，有很大收获，下面分享我的实践过程。

---

1. 软件安装

整个过程基本是从零开始，也就是说服务器没有安装任何所需软件。因为我平时会用到Python，所以第一步安装的是Anaconda，版本是Anaconda3-4.4.0-Linux-x86_64.sh. Ananconda是一个用于科学计算的Python发行版，能够方便解决多版本Python并存（后面会看到）、切换以及各种第三方包安装问题（最大的好处）。

Anaconda和整个ChIP-seq分析没关系，提到它是因为安装Anaconda后可以用BIOCONDA，能够方便安装管理生物信息软件，无需自己解决软件之间依赖关系。

安装本次实践需要的软件，包括

fastqc(0.11.5), sra-tools(2.8.1), bowtie2(2.3.2), samtools(1.5), MACS2(2.1.1.20160309), deeptools(2.5.1)

1.    conda config --add channels defaults

2.    conda config --add channels conda-forge

3.    conda config --add channels bioconda

4.    conda install fastqc

5.    conda install sra-tools

6.    conda install bowtie2

7.    conda install samtools

8.    conda install deeptools

MACS2比较特殊，因为它需要Python2.7的环境才能安装，因为我习惯使用Python3，所以安装的Anaconda版本是Python3.6的。为了安装MACS2，首先创建一个虚拟环境，在python2.7的虚拟环境中使用。

1.    conda create -n py27 python=2

2.    source activate py27

3.    pip install MACS2

用pip没用conda，是因为conda install MACS2需要装好多其他包，我就按照MACS2网站提示pip安装，发现也能用。

---

2. 数据下载

处理过程中需要下载必要的数据,参考生物信息学常见的数据下载（http://www.biotrainee.com/thread-857-1-1.html）

mm10基因组

1.    mkdir -p  ~/reference/genome/mm10

2.    cd ~/reference/genome/mm10

3.    wget http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz

4.    tar zvxf chromFa.tar.gz

5.    cat *.fa > mm10.fa

6.    rm chr*.fa

build index

1.    mkdir -p ~/reference/index/bowtie

2.    cd ~/reference/index/bowtie

3.    nohup time bowtie2-build  ~/reference/genome/mm10/mm10.fa  ~/reference/index/bowtie/mm10 1>mm10.bowtie_index.log 2>&1 &

mm10 refseq注释文件，来自http://genome.ucsc.edu/cgi-bin/hgTables。其中track选择RefSeq Genes，output format选择BED，其他参数默认。 头几行得到结果如下。

---

3. 处理数据

处理的步骤和Jimmy一模一样。 首先需要下载数据，在Jimmy原文中写的清楚，这里也就不赘述了。但是原文截取的原始数据的图好像不太对。我得到的数据如图：

PS：原文截图的确出错了，原文用的是RNA-seq实战的图片，编辑错误。

然后用fastqc看数据质量

1.    ls *fastq|xargs fastqc -t 10

用Bowtie2将fastqc结果比对到mm10基因组上去

1.    bowtie2 -p 6 -3 5 --local -x ~/reference/index/bowtie/mm10 -U ../data/SRR620205.fastq| samtools sort -O bam -o cbx7.bam

2.    bowtie2 -p 6 -3 5 --local -x ~/reference/index/bowtie/mm10 -U ../data/SRR620206.fastq| samtools sort -O bam -o suz12.bam

3.    bowtie2 -p 6 -3 5 --local -x ~/reference/index/bowtie/mm10 -U ../data/SRR620207.fastq| samtools sort -O bam -o RYBP.bam

4.    bowtie2 -p 6 -3 5 --local -x ~/reference/index/bowtie/mm10 -U ../data/SRR620208.fastq| samtools sort -O bam -o IgGold.bam

5.    bowtie2 -p 6 -3 5 --local -x ~/reference/index/bowtie/mm10 -U ../data/SRR620209.fastq| samtools sort -O bam -o IgG.bam

得到比对结果如图。大小和原文已经不一样，可能是index，bowtie2版本不同有关。

然后call peaks，首先要进入py27的虚拟环境中

1.    source activate py27

2.    macs2 callpeak -c IgGold.bam -t suz12.bam -m 10 30 -p 1e-5 -f BAM -g mm -n suz12 2>suz12.masc2.log

3.    macs2 callpeak -c IgGold.bam -t ring1B.bam -m 10 30 -p 1e-5 -f BAM -g mm -n ring1B 2>ring1B.masc2.log

4.    macs2 callpeak -c IgG.bam -t cbx7.bam -m 10 30 -p 1e-5 -f BAM -g mm -n cbx7 2>cbx7.masc2.log

5.    macs2 callpeak -c IgG.bam -t RYBP.bam -m 10 30 -p 1e-5 -f BAM -g mm -n RYBP 2>RYBP.masc2.log

结果如下。和原文结果很接近。但是原文RYBP2不知道是怎么得到的。

PS:我的RYBP2得到的方式，其实原文有说！

然后是bam转换成bw文件，并画图

1.    ls *.bam |while read id

2.    do

3.    samtools index $id $id.bai

4.    done

5.    ls *bam |while read id

6.    do

7.    file=$(basename $id )

8.    sample=${file%%.*}

9.    echo $sample

10.    bamCoverage -b $id -o $sample.bw -p 10 --normalizeUsingRPKM

11.    computeMatrix reference-point --referencePoint TSS -b 10000 -a 10000 -R ~/annotation/CHIPseq/mm10/ucsc.refseq.bed -S $sample.bw --skipZeros -o matrix1_${sample}_TSS.gz --outFileSortedRegions regions1_${sample}_genes.bed

12.    plotHeatmap -m matrix1_${sample}_TSS.gz -out ${sample}.png

13.    done

结果如下。其中跑bamCoverage会报错ring1B.bam does not appear to have an index. You MUST index the file first!

先要跑samtools index，生成.bam.bai文件。这个不知道为什么跑原代码没有生成这个文件。

最后画出profile和heatmap图。

1.    computeMatrix reference-point -p 10 --referencePoint TSS -b 2000 -a 2000 -S *bw -R ~/annotation/CHIPseq/mm10/ucsc.refseq.bed --skipZeros -o tmp4.mat.gz

2.    plotHeatmap -m tmp4.mat.gz -out tmp4.merge.png

3.    plotProfile --dpi 720 -m tmp4.mat.gz -out tmp4.profile.pdf --plotFileFormat pdf --perGroup

4.    plotHeatmap --dpi 720 -m tmp4.mat.gz -out tmp4.merge.pdf --plotFileFormat pdf

得到结果如下

和原文差别好多。后来发现是bamCoverage不需要 --normalizeUsingRPKM这个参数，调整后结果是

现在和原文结果基本一致了。

---

4.总结

本次实践过程，只是单纯把软件安装，数据下载，代码跑了一遍，其中参数的设置及含义也没有仔细查看，这个要在之后的练习中不断学习。 多谢Jimmy分享的教程，能够让读者快速进入实战分析！

编辑：思考问题的熊&&jimmy

你们也赶快实践一下吧~

参考原文：ChIP-seq实战分析

作业1：读 Jimmy 前辈的《Chip-seq 实战分析》后感：新司机驾考不完全指南