九月学徒ChIP-seq学习成果展（6万字总结）（上篇）
学徒第4周是ChIP-seq数据分析实战训练，讲义大纲文末的阅读原文，配套视频在B站：
九月学徒已经结业，表现还不错，学了几个NGS组学数据处理加上部分单细胞，随机安排的文献数据处理图表复现也完成的还不赖，昨天在生信技能树的WGCNA代码就是他写的；重复一篇WGCNA分析的文章（代码版）
因为公众号排版真的是力气活，同样我也是逼着学徒自己排版的，反正这辈子就苦这么一回！“作词作曲”都是他自己！
因为学徒提交作业实在是太长，超过了微信公众号发完限制，所以分成上下文！
首先是上游Chip-Seq
1.文章数据选取
这里紧跟jimmy大神的脚本，选择同样的文章数据，来自于《CARM1 Methylates Chromatin Remodeling Factor BAF155 to Enhance Tumor Progression and Metastasis》
作者首先实验证明了用small haripin RNA来knockout CARM1 只能达到90%的敲除效果，有趣的是，对CARM1的功能影响非常小，说明只需要极少量的CARM1就可以发挥很好的作用。
因此作者通过zinc finger nuclease这种基因组编辑技术设计了100%敲除CARM1的实验材料，这样就能比较CARM1有无时各种蛋白被催化状态了。
其中SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物的一个亚基 BAF155，非常明显的只有在CARM1这个基因完好无损的细胞系里面才能被正常的甲基化。作者证明了BAF155是CARM1这个基因非常好的一个底物， 而且通过巧妙的实验设计，证明了BAF155这个蛋白的第1064位氨基酸(R) 是 CARM1的作用位点。
因为早就有各种文献说明了SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物在癌症的重要作用， 所以作者也很自然想探究BAF155在癌症的功能详情，这里作者选择的是ChIP-seq技术。BAF155是作为SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物的一个组分，必然可以直接或者间接的结合DNA。
而ChIP-seq技术最适合来探究能直接或者间接结合DNA的蛋白的功能，所以作者构造了一种细胞系（MCF7），它的BAF155蛋白的第1064位氨基酸(R) 突变而无法被CARM1这个基因催化而甲基化，然后比较突变的细胞系和野生型细胞系的BAF155的两个ChIP-seq结果，这样就可以研究BAF155是否必须要被CARM1这个基因催化而甲基化后才能行使生物学功能。
作者用me-BAF155特异性抗体+western bloting 证明了正常的野生型MCF7细胞系里面有~74%的BAF155被甲基化。
有一个细胞系SKOV3，可以正常表达除了BAF155之外的其余14种SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物，而不管是把突变的细胞系和野生型细胞系的BAF155混在里面都可以促进染色质重构复合物的组装，所以甲基化与否并不影响这个染色质重构复合物的组装，重点应该研究的是甲基化会影响BAF155在基因组其它地方结合。
2.数据下载
考虑到原始文件过大（主要是练习的服务器配置不够），所以只下载了4个文件，就是WT和MUT的input和IP的数据。
mkdir -p /project/home/lyang/ChIP-seq_test && cd /project/home/lyang/ChIP-seq_test
mkdir rawData && cd rawData
prefetch SRR1042595 -O ./
prefetch SRR1042596 -O ./
prefetch SRR1042597 -O ./
prefetch SRR1042598 -O ./
3.格式转换——单端测序结果
mkdir 1.trans_fastq && cd 1.trans_fastq
ln -s /project/home/lyang/ChIP-seq_test/rawData/*sra ./
fastq-dump --gzip -O ~/ChIP-seq_test/1.trans_fastq SRR1042595.sra
fastq-dump --gzip -O ~/ChIP-seq_test/1.trans_fastq SRR1042596.sra
fastq-dump --gzip -O ~/ChIP-seq_test/1.trans_fastq SRR1042597.sra
fastq-dump --gzip -O ~/ChIP-seq_test/1.trans_fastq SRR1042598.sra
4.质控fastqc
mkdir /project/home/lyang/ChIP-seq_test/2.fasqc && cd /project/home/lyang/ChIP-seq_test/2.fasqc
ln -s ../1.trans_fastq/*.fastq.gz ./
fastqc -t 2 *
得到html报告可以自行查看，或者实验multiqc整合报告
5.用multiqc整合报告
cd /project/home/lyang/ChIP-seq_test/2.fasqc
mkdir multiqc
multiqc *zip -o multiqc/
6.过滤低质量的fq文件
mkdir /project/home/lyang/ChIP-seq_test/3.trim_galore && cd /project/home/lyang/ChIP-seq_test/3.trim_galore
mkdir multiqc
ln -s ../1.trans_fastq/*.fastq.gz ./
nohup trim_galore --phred33 -q 5 --length 36 --stringency 3 --fastqc -o ./ SRR1042595.fastq.gz &
nohup trim_galore --phred33 -q 5 --length 36 --stringency 3 --fastqc -o ./ SRR1042596.fastq.gz &
nohup trim_galore --phred33 -q 5 --length 36 --stringency 3 --fastqc -o ./ SRR1042597.fastq.gz &
nohup trim_galore --phred33 -q 5 --length 36 --stringency 3 --fastqc -o ./ SRR1042598.fastq.gz &
multiqc *zip -o multiqc/
因为我们只有4个样本，所以懒得写循环，也懒得写控制脚本
可以看到与修剪之前相比，%Dups和%GC比例都降低了。
7.bowtie2比对
7.1索引的构建
文中用到的参考基因组是hg19，所以去下载hg19参考基因组，然后用bowtie2构建索引。
mkdir -p /project/home/lyang/refdata/bowtie2/ && cd /project/home/lyang/refdata/bowtie2/
wget -c http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz
gunzip hg19.fa.gz
#构建索引
nohup bowtie2-build --threads 2 -f hg19.fa hg19 &
7.2进行比对得到bam文件
mkdir /project/home/lyang/ChIP-seq_test/4.mapping && cd /project/home/lyang/ChIP-seq_test/4.mapping
ln -s ../3.trim_galore/*trimmed.fq.gz ./
##可以做成一个脚本文件后用nohup bash xxx.bash &挂到后台运行
index=/project/home/lyang/refdata/bowtie2/hg19
ls *_trimmed.fq.gz|while read id
do
bowtie2 -x $index -U $id -S ${id%%_trimmed*}.sam
samtools view -bhS -q 30 ${id%%_trimmed*}.sam > ${id%%_trimmed*}.bam
samtools sort ${id%%_trimmed*}.bam -@ 10 -o ${id%%_trimmed*}.sorted.bam
samtools index -@ 10 ${id%%_trimmed*}.sorted.bam
done
上面的samtools步骤很多，可以管道起来让它简洁一点，不过需要注意的是服务器内存要好！
对bam文件进行QC
ls *sorted.bam |xargs -i samtools index {}
ls *sorted.bam | while read id ;do (samtools flagstat $id > $(basename $id ".sorted.bam").stat);done
查看成功率：
cat *stat*|grep %
结果似乎不错：
11810228 + 0 mapped (100.00% : N/A)
47955599 + 0 mapped (100.00% : N/A)
17790926 + 0 mapped (100.00% : N/A)
42749895 + 0 mapped (100.00% : N/A)
这个百分百比对成功看起来有点不真实。
8.如果有必要，则合并bam文件（这里不需要）
有些数据是同一个样本的测序数据由于被分到了不同的lane中进行测序，所以最后得到的sra文件也分成了数个。
如果存在这种情况则需要对不同的文件进行合并，但是前提是检验了不同bam文件的相关性非常好！
这里的数据单个文件就是一个完整的样本，所以无需合并，但是还是先把代码记录下来，以免日后需要使用：
ls *.bam|sed 's/_[0-9]_trimmed.bam//g' |sort -u |while read id;do samtools merge ../mergeBam/$id.merge.bam $id*.bam ;done
解释：
ls *.bam|sed 's/_[0-9]_trimmed.bam//g' |sort -u 是为了得到所以的样本名，去重以后的样本名
samtools merge ../mergeBam/$id.merge.bam $id*.bam samtools merge后接上merge后的文件名，再加上所有需要merge的$id*.bam文件
9.关于PCR重复的去除
经过bowtie2比对并将bam文件排序，构建索引后得到的文件：
557M Oct 4 10:32 SRR1042595.sorted.bam
4.1M Oct 4 15:55 SRR1042595.sorted.bam.bai
2.6G Oct 4 12:49 SRR1042596.sorted.bam
2.3M Oct 4 15:56 SRR1042596.sorted.bam.bai
889M Oct 4 13:38 SRR1042597.sorted.bam
3.2M Oct 4 15:56 SRR1042597.sorted.bam.bai
2.4G Oct 4 15:36 SRR1042598.sorted.bam
2.3M Oct 4 15:57 SRR1042598.sorted.bam.bai
去除PCR重复并建立索引
mkdir /project/home/lyang/ChIP-seq_test/5.RemovePCRduplicated/ && cd /project/home/lyang/ChIP-seq_test/5.RemovePCRduplicated/
cp ../4.mapping/*sort* ./
ls *sorted.bam | while read id ;do (nohup samtools markdup -r $id $(basename $id ".sorted.bam").rmdup.bam & );done
ls *.rmdup.bam |xargs -i samtools index {}
#对rmdup.bam文件进行QC
ls *.rmdup.bam | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done
QC质控结果：
cat *stat*|grep %
4159545 + 0 mapped (100.00% : N/A)
46290036 + 0 mapped (100.00% : N/A)
9519997 + 0 mapped (100.00% : N/A)
40969756 + 0 mapped (100.00% : N/A)
去除重复前后比较（主要看文件大小）：
ls -lh *bam|cut -d" " -f5-10
230M Oct 4 16:29 SRR1042595.rmdup.bam
557M Oct 4 10:32 SRR1042595.sorted.bam
2.6G Oct 4 16:33 SRR1042596.rmdup.bam
2.6G Oct 4 12:49 SRR1042596.sorted.bam
527M Oct 4 16:30 SRR1042597.rmdup.bam
889M Oct 4 13:38 SRR1042597.sorted.bam
2.3G Oct 4 16:33 SRR1042598.rmdup.bam
2.4G Oct 4 15:36 SRR1042598.sorted.bam
可以看到至少去除了一半的序列。
10.利用macs2来callpeak
10.1 首先了解关于macs2这个软件：
macs2包含一系列的子命令，其中最主要的就是callpeak， 官方提供了使用实例
macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01
一般而言，我们照葫芦画瓢，按照这个实例替换对应部分就行了，介绍一下各个参数的意义
-t: 实验组的输出结果
-c: 对照组的输出结果
-f: -t和-c提供文件的格式，可以是”ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一个。如果不提供这项，就是自动检测选择。
-g: 基因组大小， 默认提供了hs, mm, ce, dm选项， 不在其中的话，比如说拟南芥，就需要自己提供了。
-n: 输出文件的前缀名
-B: 会保存更多的信息在bedGraph文件中，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores
-q: q值，也就是最小的PDR阈值， 默认是0.05。q值是根据p值利用BH计算，也就是多重试验矫正后的结果。
-p：这个是p值，指定p值后MACS2就不会用q值了。
-m: 和MFOLD有关，而MFOLD和MACS预构建模型有关，默认是5：50，MACS会先寻找100多个peak区构建模型，一般不用改，因为我们很大概率上不会懂。
10.2 应用到自己的实例中：
首先确定哪些样本是当作control来消除噪音，哪些是treatment来真正寻找peaks的。
可以看到SRR104259{4，6，8}，SRR1042600都是作为input的control组。
而SRR104259{3，5，7，9}都是作为treatment组。
但是不知道这里control和treatment组是不是存在一一对应的关系，就当作是一一对应的关系往下做吧
mkdir /project/home/lyang/ChIP-seq_test/6.macs2-callpeak && cd /project/home/lyang/ChIP-seq_test/6.macs2-callpeak
ln -s ../5.RemovePCRduplicated/*.bam ./
nohup macs2 callpeak -c SRR1042596.rmdup.bam -t SRR1042595.rmdup.bam -f BAM -g hs -n Xu_MUT_rep2_rmdup &
nohup macs2 callpeak -c SRR1042598.rmdup.bam -t SRR1042597.rmdup.bam -f BAM -g hs -n Xu_WT_rep1_rmdup &
nohup macs2 callpeak -c SRR1042596.sorted.bam -t SRR1042595.sorted.bam -f BAM -g hs -n Xu_MUT_rep2_sort &
nohup macs2 callpeak -c SRR1042598.sorted.bam -t SRR1042597.sorted.bam -f BAM -g hs -n Xu_WT_rep1_sort &
-c 后面加上control组样本名
-t 后面加上treatment组样本名
-f {AUTO,BAM,SAM,BED,ELAND,ELANDMULTI,ELANDEXPORT,BOWTIE,BAMPE,BEDPE}
后面加上format 输出文件的格式
-q QVALUE 后面加上人为设定的q值。默认为0.05
-B, --bdg 加上这个选项则表示需要输出extended fragment pileup和local lambda tracks (two files)
-n NAME 指定输出文件名
-g GSIZE Effective genome size，可以是具体数字，也可以用简写：
'hs' for human (2.7e9), 'mm' for mouse (1.87e9), 'ce' for C. elegans (9e7)
'dm' for fruitfly (1.2e8), Default:hs
10.3 macs2结果输出
macs2的输出文件解读
得到的文件有：
10.3.1 NAME_peaks.xls
虽然后缀名是xls，但实际上就是一个普通的文本文件。包含peak信息的tab分割的文件，前几行会显示callpeak时的命令。输出信息包含：
染色体号
peak起始位点
peak结束位点
peak区域长度
peak的峰值位点（summit position）
peak 峰值的属性（包括pileup峰高和可信度）（pileup height at peak summit, -log10(pvalue) for the peak summit）都是值越高越好
peak的富集倍数（相对于random Poisson distribution with local lambda）
10.3.2 NAME_peaks.narrowPeak
narrowPeak文件是BED6+4格式，可以上传到UCSC浏览。输出文件每列信息分别包含：
1；染色体号
2：peak起始位点
3：结束位点
4：peak name
5：int(-10*log10qvalue)可信度，值越大越好
6 ：正负链
7：fold change
8：-log10pvalue可信度，值越大越好
9：-log10qvalue可信度，值越大越好
10：relative summit position to peak start相对于起始位点来说peaks峰值的位置
bed格式的文件和其他Peak文件唯一不同的是——Peak文件中，由于以1为基，而bed文件是以0为基，所以peak的起始位置需要减1才是bed格式的文件。
10.3.3 NAME_summits.bed
BED格式的文件，包含peak的summits位置，第5列是-log10pvalue。如果想找motif，推荐使用此文件。
能够直接载入UCSC browser，用其他软件分析时需要去掉第一行。？？？这里有点奇怪
10.3.4 NAME_model.r
NAME_model.r能通过 $ Rscript NAME_model.r作图，得到是基于你提供数据的peak模型。
可以看到非常贴心的帮我们写好了脚本，直接运行即可出pdf图。
10.4 macs2背后的统计学原理
参考Chip-seq处理流程
我们对ChIP-seq测序数据寻找peaks的本质就是得到所有测序数据比对在全基因组之后在基因组上测序深度里面寻找比较突出的部分。比如对WES数据来说，各个外显子，或者外显子的5端到3端，理论上测序深度应该是一致的，都是50X~200X，画一个测序深度曲线，应该是近似于一条直线。但对我们的ChIP-seq测序数据来说，在所捕获的区域上面，理论上测序深度是绝对不一样的，应该是近似于一个山峰。
而那些覆盖度高的地方，山顶，就是我们的IP所结合的热点，也就是我们想要找的peaks。我们通常说的ChIP-seq测序的IP，可以是各个组蛋白上各修饰位点对应的抗体，或者是各种转录因子的抗体等等。
热点是这样一些位置，这些位置多次被测得的read所覆盖（我们测的是一个细胞群体，read出现次数多，说明该位置被TF结合的几率大）。那么，read数达到多少才叫多？这就要用到统计检验来分析。假设TF在基因组上的分布是没有任何规律的，那么，测序得到的read在基因组上的分布也必然是随机的，某个碱基上覆盖的read的数目应该服从二项分布。
TF在基因组上的结合其实是一个随机过程，基因组的每个位置其实都有机会结合某个TF，只是概率不一样，说白了，peak出现的位置，是TF结合的热点，而peak-calling就是为了找到这些热点。当n很大，p很小时，二项分布可以近似用泊松分布替代，在这里：
n是测序得到的read总数目，l是单个read的长度，s是基因组的大小。
我们可以算出在某个置信概率（如0.00001）下，随机情况下，某个碱基上可以覆盖的read的数目的最小值，当实际观察到的read数目超过这个值（单侧检验）时，我们认为该碱基是TF的一个结合热点。反过来，针对每一个read数目，我们也可以算出对应的置信概率P。
但是，这只是一个简化的模型，实际情况要复杂好多。比如，由于测序、mapping过程内在的偏好性，以及不同染色质间的差异性，相比全基因组，某些碱基可能内在地会被更多的read所覆盖，这种情况得到的很多peak可能都是假的。
MACS考虑到了这一点，当对某个碱基进行假设检验时，MACS只考虑该碱基附近的染色质区段（如10k），此时，上述公式中n表示附近10k区间内的read数目，s被置为10k。当有对照组实验（Control，相比实验组，没有用抗体捕获TF，或用了一个通用抗体）存在时，利用Control组的数据构建泊松分布，当没有Control时，利用实验组，稍大一点的局部区间（比如50k）的数据构建泊松分布。
read只是跟随着TF一起沉淀下来的DNA fragment的末端，read的位置并不是真实的TF结合的位置。所以在peak-calling之前，延伸read是必须的。不同TF大小不一样，对read延伸的长度也理应不同。我们知道，测得的read最终其实会近似地平均分配到正负链上，这样，对于一个TF结合热点而言，read在附近正负链上会近似地形成“双峰”。
MACS会以某个window size扫描基因组，统计每个window里面read的富集程度，然后抽取（比如1000个）合适的（read富集程度适中，过少，无法建立模型，过大，可能反映的只是某种偏好性）window作样本，建立“双峰模型”。最后，两个峰之间的距离就被认为是TF的长度D，每个read将延伸D/2的长度。
当有对照组实验存在时，MACS会进行两次peak calling。第一次以实验组（Treatment）为实验组，对照组为对照组，第二次颠倒，以实验组为对照组，对照组为实验组。之后，MACS对每一个P计算了相应的FDR（False Discovery Rate）值：
Control-p表示第二次peak calling（颠倒的）得到的置信概率小于P的peak的个数。
Treatment-p表示第一次peak calling得到的置信概率小于P的peak的个数。
FDR综合利用了实验组和对照组的信息，显然，FDR越小越好。
这一段摘抄自plob的教程。
10.5 自己得到的输出结果如下：
#sorted.bam文件得到的结果
3.5K Xu_MUT_rep2_sort_peaks.narrowPeak
5.2K Xu_MUT_rep2_sort_peaks.xls
2.4K Xu_MUT_rep2_sort_summits.bed
4.7K Xu_WT_rep1_sort_peaks.narrowPeak
6.5K Xu_WT_rep1_sort_peaks.xls
3.2K Xu_WT_rep1_sort_summits.bed
#rmdup.bam文件得到的结果
3.6K Xu_MUT_rep2_rmdup_peaks.narrowPeak
5.2K Xu_MUT_rep2_rmdup_peaks.xls
2.5K Xu_MUT_rep2_rmdup_summits.bed
4.8K Xu_WT_rep1_rmdup_peaks.narrowPeak
6.6K Xu_WT_rep1_rmdup_peaks.xls
3.3K Xu_WT_rep1_rmdup_summits.bed
##统计去除PCR重复和没有去重的bam文件找到的peaks数目差别：结果居然也是没有差别
wc -l *bed
43 Xu_MUT_rep2_rmdup_summits.bed
43 Xu_MUT_rep2_sort_summits.bed
58 Xu_WT_rep1_rmdup_summits.bed
58 Xu_WT_rep1_sort_summits.bed
BED包含有3个必须的字段和9个可选字段：
必须字段：
1 chrom - 染色体名字
2 chromStart - 染色体起始位点
3 chromEnd - 染色体终止位点
可选字段：
4 name - 名字
5 score - 分值(0-1000), 用于genome browser展示时上色。
6 strand - 正负链，对于ChIPseq数据来说，一般没有正负链信息。
7 thickStart - 画矩形的起点
8 thickEnd - 画矩形的终点
9 itemRgb - RGB值
10 blockCount - 子元件（比如外显子）的数目
11 blockSizes - 子元件的大小
12 blockStarts - 子元件的起始位点
而peak calling的MACS输出bed文件默认是前5列，其中第五列指的是片段堆积的峰高：
chr1 121485178 121485179 Xu_MUT_rep2_rmdup_peak_1 169.64561
chr10 42385559 42385560 Xu_MUT_rep2_rmdup_peak_2 28.26733
chr10 42387315 42387316 Xu_MUT_rep2_rmdup_peak_3 21.20116
chr10 42392730 42392731 Xu_MUT_rep2_rmdup_peak_4 24.90387
chr10 42599794 42599795 Xu_MUT_rep2_rmdup_peak_5 13.26007
chr11 18127886 18127887 Xu_MUT_rep2_rmdup_peak_6 5.11958
chr12 107104197 107104198 Xu_MUT_rep2_rmdup_peak_7 4.89626
chr16 46392027 46392028 Xu_MUT_rep2_rmdup_peak_8 18.40010
~~~内容过多，无法全部展示~~~
所以在ChIPseeker是画peak coverage的函数covplot要有个weightCol的参数。
11.使用deeptools是进行可视化
deeptools工具介绍
BAM神器--Deeptools使用指南
对于deeptools里的任意子命令，都支持
--help看帮助文档
--numberOfProcessors/-p设置多核处理
--region/-r CHR:START:END处理部分区域。
还有一些过滤用参数部分子命令可用，如ignoreDuplicates,minMappingQuality,samFlagInclude,samFlagExclue.
官方文档见http://deeptools.readthedocs.io/en/latest/index.html, 下面按照用法引入不同的工具。
BAM转换为bigWig或bedGraph
BAM文件是SAM的二进制转换版，bigWig是wig或bedGraph的二进制版，存放区间的坐标轴信息和相关计分(score)，主要用于在基因组浏览器上查看数据的连续密度图，可用wigToBigWig从wiggle进行转换。
为什么要用bigWig? 主要是因为BAM文件比较大，直接用于展示时对服务器要求较大。因此在GEO上仅会提供bw,即bigWig下载，便于下载和查看。如果真的感兴趣，则可以下载原始数据进行后续分析。
deeptools提供bamCoverage和bamCompare进行格式转换。
为了能够比较不同的样本，需要对先将基因组分成等宽分箱(bin)，统计每个分箱的read数，最后得到描述性统计值。
bamCoverage对于单个样本，可以用SES方法进行标准化。
bamCompare 对于两个样本，描述性统计值可以是两个样本的比率，或是比率的log2值，或者是差值。
bamCoverage的基本用法
bamCoverage -e 170 -bs 10 -b ap2_chip_rep1_2_sorted.bam -o ap2_chip_rep1_2.bw
得到的bw文件就可以送去IGV/Jbrowse进行可视化。这里的参数仅使用了-e/--extendReads和-bs/--binSize即拓展了原来的read长度，且设置分箱的大小。其他参数还有
--outFileFormat|-of {bigwig,bedgraph}：指定输出文件的格式(default: bigwig)-b/--bam BAM file：输入文件格式为bam-o/--outFileName FILENAME：输出文件名-e/--extendReads int：拓展了原来的read长度-bs/--binSize int：设置分箱的大小-p/--numberOfProcessors int：设置线程数-r/--region CHR:START:END: 选取某个区域统计--smoothLength: 通过使用分箱附近的read对分箱进行平滑化
如果为了其他结果进行比较，还需要进行标准化，deeptools提供了如下参数：
--scaleFactor: 缩放系数--normalizeUsing {RPKM,CPM,BPM,RPGC,None}: 选择标准化的方法
如果需要以100为分箱，并且标准化到1x，且仅统计某一条染色体区域的正链，输出格式为bedgraph,那么命令行可以这样写：
bamCoverage -e 170 -bs 100 -of bedgraph -r Chr4:12985884:12997458 -b 02-read-alignment/ap2_chip_rep1_1_sorted.bam -o chip.bedgraph
bamCompare和bamCoverage类似，只不过需要提供两个样本，并且采用SES方法进行标准化，于是多了--ratio参数。
把bam文件转为bw文件
cd /project/home/lyang/ChIP-seq_test/5.RemovePCRduplicated
ls *.bam |xargs -i samtools index {} 进行转换前需要先建立index
由于电脑配置问题，无法使用nohup command &这种形式循环运行，所以可以将下面的代码写入一个脚本中，再挂到后台运行。
cat >bam2bw.command
ls *.sorted.bam |while read id;do
bamCoverage -p 10 --normalizeUsing CPM -b $id -o ${id%%.*}.bw
done
ls *.rmdup.bam |while read id;do
bamCoverage -p 10 --normalizeUsing CPM -b $id -o ${id%%.*}.rmdup.bw
done
ctrl+c
nohup bash bam2bw.command &
mkdir bwfile
mv *bw bwfile
载入IGV中进行查看：
查看bam文件
找一个样本SRR1042595.sorted.bam（未去重）和SRR1042595.rmdup.bam（去除PCR重复后）导入IGV中进行查看去除PCR重复前后的变化：
可以看到这个样本没有去除重复前有很多PCR重复。
因为这些片段都是经过PCR重复扩增而来的，所以是一模一样的。因为都是一条母链复制而来的。
查看同一个样本的不同格式
bam，bw，bed文件
可以看到在rmdup后，bam文件上变化挺大的，但是奇怪的是在最后的callpeaks的bed文件居然差别不大。且bw文件也相差不大。
查看bed文件
可以先对bed文件进行处理，这样就可以把结果直接复制到IGV中，得到peaks的坐标地址了。
awk '{print $1":"$2"-"$3}' Xu_WT_rep1_sort_summits.bed
根据上面Xu_WT_rep1_sort_summits.bed文件中的结果去IGV中查看，看到在Xu_WT_rep1样本中有一个peaks。
再找一个Xu_MUT_rep2_rmdup_summits.bed文件中Xu_MUT_rep2样本的peaks。
如果igv载入非常慢的时候，可以考虑点击reload session这个选项：
查看TSS附件信号强度：
背景介绍
在 真核cDNA序列详解 有详细介绍，关于UTR，CDS，内含子，外显子以及TSS之间的关系：
启动子（promoter）：与RNA聚合酶结合并能起始mRNA合成的序列。
转录起始点（TSS）：转录时，mRNA链第一个核苷酸相对应DNA链上的碱基，通常为一个嘌呤。
UTR（Untranslated Regions)：即非翻译区，是信使RNA（mRNA）分子两端的非编码片段。
5'-UTR从mRNA起点的甲基化鸟嘌呤核苷酸帽延伸至AUG起始密码子，3'-UTR从编码区末端的终止密码子延伸至多聚A尾巴（Poly-A）的末端。
下载TSS.bed文件
方法一：[https://sourceforge.net/projects/seqminer/files/Reference%20coordinate/](https://sourceforge.net/projects/seqminer/files/Reference coordinate/)（结果证明这个文件似乎有问题，后面用computeMatrix reference-point调用时总是报错）
方法二：去UCSC官网（推荐在这里下载）
进入table browser工具后进行选择：
group这个选项我不太确定，但是大家都选的这个我就没改了。
track里可以选择NCBI RefSeq或者UCSC gene，这里选了NCBI
group一般不变
table里如果track选择NCBI RefSeq，这里就选择RefSeq；如果track选择UCSC gene，这里就选knownGene
output format根据自己的需求选择
file type returned这里选gzip compressed，这样就可以下载到压缩包格式的输出文件。选text则下载文本格式。
output file一定要写上一个文件名字，如果为空则后面无法下载，而只能在浏览器上查看。
最后点击get output即可。
这里我选择了它默认的whole gene。点击get bed即可开始下载文件。
TSS附件信号强度画图
在tss前后2kb的位置批量画图
mkdir /project/home/lyang/ChIP-seq_test/hg19.tss.bed && cd /project/home/lyang/ChIP-seq_test/hg19.tss.bed
bed=/project/home/lyang/ChIP-seq_test/hg19.tss.bed/hg19.tss.bed
ln -s ../5.RemovePCRduplicated/bwfile ./
## 在tss前后2kb的位置画图
bed=/trainee1/gz15/tss/hg19.tss.bed
for id in bwfile/*bw;
do
file=$(basename $id )
sample=${file%%.bw}
echo $sample
computeMatrix reference-point --referencePoint TSS -p 15 \
-b 2000 -a 2000 \
-R $bed \
-S $id \
--skipZeros -o matrix1_${sample}_TSS_2K.gz \
--outFileSortedRegions regions1_${sample}_TSS_2K.bed
## both plotHeatmap and plotProfile will use the output from computeMatrix
plotHeatmap -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Heatmap_2K.png
plotHeatmap -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Heatmap_2K.pdf --plotFileFormat pdf --dpi 720
plotProfile -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Profile_2K.png
plotProfile -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Profile_2K.pdf --plotFileFormat pdf --perGroup --dpi 720
done
computeMatrix reference-point函数的参数
--referencePoint {TSS,TES,center} The reference point for the plotting could be either the region start (TSS), the region end (TES) or the center of the region
-b/--beforeRegionStartLength INT bp：区域前多少bp
-a/--afterRegionStartLength INT bp：区域后多少bp
-R/--regionsFileName File：指定bed文件
-S/--scoreFileName File：指定bw文件
--skipZeros：默认为False
-o/-out/--outFileName OUTFILENAME：指定输出文件名
--outFileSortedRegions BED file
plotHeatmap和plotProfile函数参数
-m/--matrixFile MATRIXFILE：指定matrix名字（computeMatrix的输出文件名）
-o/-out/--outFileName OUTFILENAME：指定输出文件名，根据后缀自动决定图像的格式
--plotFileFormat {"png", "eps", "pdf", "plotly"，"svg"}：指定输出文件格式
--dpi DPI：指定dpi
--perGroup：默认是画一个样本的所有区域群的图，使用这个选项后则按区域群画所有样本的情况
还可以给多个bed文件来绘图，还可以画genebody的图，因为原理一样，我就不做过多介绍啦。
关于genebody，等回去后有空了看看这篇文献：Gene body methylation can alter gene expression and is a therapeutic target in cancer.
上面的批量代码其实就是为了统计全基因组范围的peak在基因特征的分布情况，也就是需要用到computeMatrix计算，用plotHeatmap以热图的方式对覆盖进行可视化，用plotProfile以折线图的方式展示覆盖情况。
computeMatrix具有两个模式:scale-region和reference-point。前者用来信号在一个区域内分布，后者查看信号相对于某一个点的分布情况。无论是那个模式，都有有两个参数是必须的，-S是提供bigwig文件，-R是提供基因的注释信息。还有更多个性化的可视化选项。
TSS附件信号强度结果解读
上面程序跑完后的结果导出到电脑上后：
profile图（相当于Heatmap图中最上面的折线图单独展示）
Heatmap图
拿一个热图做具体说明：
可以认为在距离TSS上游x处，所有基因normalization后的信号值的均值为0.060。
对了，昨天的九月学徒转录组学习成果展（3万字总结）（上篇） 也是我写的！
一定要继续关注哦，下期更精彩！
学徒写在最后
首先感谢Jimmy老师的教程和代码，基本上只要跟着一步步学下来，肯定能复现漂亮的图，但是其中的原理需要自己仔细研究和领会。
另外，非常感谢jimmy老师对我耐心的指导和引导，当我遇到问题时，能一下了解我遇到的代码问题在哪里，比我百度谷歌半天教程都有用。
学徒已经做的很优秀了，一个月的时间总是短暂的，但学习的脚步不能停下，希望他回去以后能有更多的学习成果跟大家分享！
当然，如果你感兴趣学徒培养，也可以看看招学徒的通知：
你可以先看看我是如何培养学徒的：