我的ChIP-Seq(3):比对和call peaks

比对（1）使用bowtie2建立index

建立index是一劳永逸的事情，还有更简单的方法是从bowtie2的官网直接下载。
这一步很简单，用的是bowtie2-2.3.4，命令为：
bowtie2-bulid ref.fa index
注意：
此处index不仅仅是生成路径，而且要命名index文件的前缀，如本次hg19_index，并不是文件夹的名字。
结果：
产生6个结果文件；耗时：4h；日志文件显示：build a small index
百度了原因，说是小于4G的基因组都会生成small index，但是同样能用而且是后期比对时软件自动识别。说明如下：

small index
值得注意的是，若程序没有完全运行结束，只会产生4个结果文件，会误以为没运行成功，所以要耐心等。

比对（2）：mapping

首先将参考基因组移动到index所在的文件夹，并将名字前缀改为与index一致。以我自己的数据（双端）为例：

bowtie2 -p 2   ###2个线程
-x /refdata_bowtie_index/hg19_index    ###index索引文件的前缀，-x小写
-1 ~/rawdata_trim/clean_R1_paired.fastq.gz   ###R1文件，可多个用逗号隔开，支持fastq.gz
-2 ~/rawdata_trim/clean_R2_paired.fastq.gz   ###R2文件，可多个，与R1对应起来
-S ~/cleandata_bowtie_sam/clean.sam    ###所生成的sam文件的名称，-S大写

结果：得到比对后的sam文件，任务日志会显示比对率等信息。

peak calling

本次选用是MACS2，值得注意的是：根据不同组蛋白修饰类型选择模式，分narrow和broad两种calling，例如H3K27ME3宽峰模式而H3K4ME3窄峰模式（如图）。这是基础知识，一定不要搞错。

narrow and broad
常规的peak calling
macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n fileprefix -B -q 0.01
较宽的peak calling
macs2 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1
较宽的peakcalling自己整理(-B可以不加)：
macs2 callpeak -t ChIP.sam -c Control.sam -f SAM --broad -g hs -n fileprefix -B --broad-cutoff 0.1

参数：
-g 基因组大小，人比较好，有内置的大小，一般是基因组大小的80%，其他物种需要估算。
-f 格式，sam bam都行，但是要告诉它
-n 给你的项目取名字，就是前缀
-B 以bgh等形式储存
-q q值(最小的FDR)的阈值，默认0.05。

以上。

后记：
整理日志文件时发现报错如图

报错截图

还是应该好好读一下软件的文章和说明。挖个坑，解决了问题再来更改。