BAM神器--Deeptools使⽤指南
⽣信媛
尽管以前推送过⼀篇相关的⼯具说明深海⼯具-DeepTools，但是最近处理ChIP数据时，我还是按照⾃⼰学习过程写了⼀遍。各有长处吧。
账簿分类
如何使⽤deeptools处理BAM数据
总体介绍
deeptools是基于Python开发的⼀套⼯具，⽤于处理诸如RNA-q, ChIP-q, MNa-q, ATAC-q等⾼通量数据。⼯具分为四个模块
BAM和bigWig⽂件处理
质量控制
热图和其他描述性作图
其他
当然也可以简单分为两个部分：数据处理和可视化。
对于deeptools⾥的任意⼦命令，都⽀持--help看帮助⽂档，--numberOfProcessors/-p设置多核处理，--region/-r CHR:START:END处理部分区域。还有⼀些过滤⽤参数部分⼦命令可⽤，
如ignoreDuplicates,minMappingQuality,samFlagInclude,samFlagExclue.
官⽅⽂档见adthedocs.io/en/latest/index.html, 下⾯按照⽤法引⼊不同的⼯具。
后续演⽰的数据来⾃于Orchestration of the Floral Transition and Floral
Development in Arabidopsis by the Bifunctional Transcription Factor APETALA2，
free money
如要重复请⾃⾏下载⽐对。
BAM转换为bigWig或bedGraph
BAM⽂件是SAM的⼆进制转换版，应该都知道。那么bigWig格式是什么？bigWig是wig或bedGraph的⼆进制版，存放区间的坐标轴信息和相关计分(score)，主要⽤于在基因组浏览器上查看数据的连续密度图，可⽤wigToBigWig从wiggle进⾏转换。
time after time歌词bedGraph和wig格式是什么? USCS的帮助⽂档称这两个格式数是已经过时的基因组浏览器图形轨展⽰格式，前者展⽰稀松型数据，后者展⽰连续性数据。⽬前推荐使⽤bigWig/bigBed这两种格式取代前两者。
如何收缩脸部毛孔>很郁闷
为什么要⽤bigWig? 主要是因为BAM⽂件⽐较⼤，直接⽤于展⽰时对服务器要求较⼤。因此在GEO上仅会提供bw，即bigWig下载，便于下载和查看。如果真的感兴趣，则可以下载原始数据进⾏后续分析。
deeptools提供bamCoverage和bamCompare进⾏格式转换，为了能够⽐较不同的样本，需要对先将基因组分成等宽分箱(bin)，统计每个分箱的read数，最后得到描述性统计值。对于两个样本，描述性统计值可以是两个样本的⽐率，或是⽐率的log2值，或者是差值。如果是单个样本，可以⽤SES⽅法进⾏标准化。
gariabamCoverage的基本⽤法
tesorobamCoverage -e 170 -bs 10 -b ap2_chip_rep1_2_sorted.bam -o ap2_chip_rep1_2.bw# ap2_chip_rep1_2_sorted.bam 是前期⽐对得到的BAM⽂件
得到的bw⽂件就可以送去IGV/Jbrow进⾏可视化。这⾥的参数仅使⽤了-e/--extendReads和-bs/--binSize即拓展了原来的read长度，且设置分箱的⼤⼩。其他参数还有
--filterRNAstrand {forward, rever}: 仅统计指定正链或负链
--region/-r CHR:START:END: 选取某个区域统计
--smoothLength: 通过使⽤分箱附近的read对分箱进⾏平滑化
如果为了其他结果进⾏⽐较，还需要进⾏标准化，deeptools提供了如下参数：
--scaleFactor: 缩放系数
`—normalizeUsingRPKMReads``: Per Kiloba per Million mapped reads (RPKM)标准化
--normalizeTo1x: 按照1x测序深度(reads per genome coverage, RPGC)进⾏标准化
--ignoreForNormalization：指定那些染⾊体不需要经过标准化
影子英文
如果需要以100为分箱，并且标准化到1x，且仅统计某⼀条染⾊体区域的正链，输出格式为bedgraph，那么命令⾏可以这样写
bamCoverage -e 170 -bs 100 -of bedgraph -r Chr4:12985884:12997458 --normalizeTo1x 100000000 -b 02-read-alignment/ap2_chip_rep1_1_sorted.bam -o chip.bedgraph
成考和自考的区别
bamCompare和bamCoverage类似，只不过需要提供两个样本，并且采⽤SES⽅法进⾏标准化，于是多了--ratio参数。
多样本分析
这部分内容主要分析处理组不同重复间的相关程度，会⽤
到multiBamSummary、plotCorrelation和plotPCA三个模块。。主要⽬的是看下对照组和处理组中的组间差异和组内相似性。
如果上⼀步把BAM转换成BW, 那么multiBamSummary可以⽤multiBigWigSummary
替代
# 统计reads在全基因组范围的情况multiBamSummary bins -bs 1000 --bamfiles 02-read-
alignment/ap2_chip_rep1_1_sorted.bam 02-read-alignment/ap2_chip_rep1_2_sorted.bam 02-read-
alignment/ap2_chip_rep1_3_sorted.bam 02-read-alignment/ap2_chip_rep2_1_sorted.bam 02-read-
alignment/ap2_ctrl_rep1_1_sorted.bam 02-read-alignment/ap2_ctrl_rep1_2_sorted.bam 02-read-
alignment/ap2_ctrl_rep2_1_sorted.bam --extendReads 130 -out treat_results.npz# 散点图plotCorrelation -in
treat_results.npz -o treat_results.png --corMethod spearman -p scatterplot# 热图plotCorrelation -in treat_results.npz -o treat_results_heatmap.png --corMethod spearman -p heatmap# 主成分分析plotPCA -in treat_results.npz -o pca.png
根据下图不难发现，组内的不同技术重复间差异性⼩，⽽组内中的两个⽣物学重复看起来只能说还⾏，但是差异还是⼩于组间的差异。
但是看主成分分析结果，总感觉哪⾥不对劲。不过这仅仅是看总体的分布情况，⽽不是使⽤差异peak进⾏主成分分析，也不知道这样说对不对。
peak分布可视化
为了统计全基因组范围的peak在基因特征的分布情况，需要⽤到computeMatrix计算，
⽤plotHeatmap以热图的⽅式对覆盖进⾏可视化，⽤plotProfile以折线图的⽅式展⽰覆盖情况。
computeMatrix具有两个模式:scale-region和reference-point。前者⽤来信号在⼀个区域内分布，后者查看信号相对于某⼀个点的分布情况。
⽆论是那个模式，都有有两个参数是必须的，-S是提供bigwig⽂件，-R是提供基因的注释信息。
scale-regions模式
computeMatrix scale-regions \ # 选择模式-b 3000 -a 5000 \ # 感兴趣的区域，-b上游，-a下游-R
~/reference/gtf/TAIR10/TAIR10_GFF3_genes.bed \-S 03-read-coverage/ap2_chip_rep1_1.bw \--skipZeros \--outFileNameMatrix 03-read-coverage/matrix1_ap2_chip_rep1_1_scaled.tab \ # 输出为⽂件⽤于plotHeatmap, plotProfile--outFileSortedRegions 03-read-coverage/regions1_ap2_chip_re1_1_genes.bed
reference-point模式
computeMatrix reference-point \ # 选择模式--referencePoint TSS \ # 选择参考点: TES, center-b 3000 -a 5000 \ # 感兴趣的区域，-b上游，-a下游-R ~/reference/gtf/TAIR10/TAIR10_GFF3_genes.bed \-S 03-read-
coverage/ap2_chip_rep1_1.bw \--skipZeros \-out 03-read-coverage/matrix1_ap2_chip_rep1_ \ # 输出为⽂件⽤于plotHeatmap, plotProfile--outFileSortedRegions 03-read-coverage/ons1regions1_ap2_chip_re1_1_genes.bed 结果可视化
可视化的⽅法有两种，⼀种是轮廓图，⼀种是热图。两则都提供了⾜够多的参数对结果进⾏细节上的修改。
海地是哪个国家plotProfile -m matrix1_ap2_chip_rep1_ \-out ExampleProfile1.png \--numPlotsPerRow 2 \--plotTitle 'Test data profile'plotHeatmap -m matrix1_ap2_chip_rep1_ \-out ExampleHeatmap1.png \
PS：如果是找到的peak可⽤R包chipSeeker进⾏可视化。
