如何使用deeptools處理BAM數(shù)據(jù)

總體介紹

deeptools是基于Python開發(fā)的一套工具，用于處理諸如RNA-seq, ChIP-seq, MNase-seq, ATAC-seq等高通量數(shù)據(jù)。工具分為四個模塊

• BAM和bigWig文件處理
• 質(zhì)量控制
• 熱圖和其他描述性作圖
• 其他

當然也可以簡單分為兩個部分：數(shù)據(jù)處理和可視化。

對于deeptools里的任意子命令，都支持--help看幫助文檔，--numberOfProcessors/-p設置多核處理，--region/-r CHR:START:END處理部分區(qū)域。還有一些過濾用參數(shù)部分子命令可用，如ignoreDuplicates,minMappingQuality,samFlagInclude,samFlagExclue.

官方文檔見http://deeptools.readthedocs.io/en/latest/index.html, 下面按照用法引入不同的工具。

整體介紹

后續(xù)演示的數(shù)據(jù)來自于Orchestration of the Floral Transition and Floral Development in Arabidopsis by the Bifunctional Transcription Factor APETALA2，如要重復請自行下載比對。

BAM轉換為bigWig或bedGraph

BAM文件是SAM的二進制轉換版，應該都知道。那么bigWig格式是什么？bigWig是wig或bedGraph的二進制版，存放區(qū)間的坐標軸信息和相關計分(score)，主要用于在基因組瀏覽器上查看數(shù)據(jù)的連續(xù)密度圖，可用wigToBigWig從wiggle進行轉換。

bedGraph和wig格式是什么? USCS的幫助文檔稱這兩個格式數(shù)是已經(jīng)過時的基因組瀏覽器圖形軌展示格式，前者展示稀松型數(shù)據(jù)，后者展示連續(xù)性數(shù)據(jù)。目前推薦使用bigWig/bigBed這兩種格式取代前兩者。

為什么要用bigWig? 主要是因為BAM文件比較大，直接用于展示時對服務器要求較大。因此在GEO上僅會提供bw,即bigWig下載，便于下載和查看。如果真的感興趣，則可以下載原始數(shù)據(jù)進行后續(xù)分析。

bigWig更適合展示

deeptools提供bamCoverage和bamCompare進行格式轉換，為了能夠比較不同的樣本，需要對先將基因組分成等寬分箱(bin)，統(tǒng)計每個分箱的read數(shù)，最后得到描述性統(tǒng)計值。對于兩個樣本，描述性統(tǒng)計值可以是兩個樣本的比率，或是比率的log2值，或者是差值。如果是單個樣本，可以用SES方法進行標準化。

bamCoverage的基本用法

bamCoverage -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/Jbrowse進行可視化。 這里的參數(shù)僅使用了-e/--extendReads和-bs/--binSize即拓展了原來的read長度，且設置分箱的大小。其他參數(shù)還有

• --filterRNAstrand {forward, reverse}: 僅統(tǒng)計指定正鏈或負鏈
• --region/-r CHR:START:END: 選取某個區(qū)域統(tǒng)計
• --smoothLength: 通過使用分箱附近的read對分箱進行平滑化

如果為了其他結果進行比較，還需要進行標準化，deeptools提供了如下參數(shù)：

• --scaleFactor: 縮放系數(shù)
• `--normalizeUsingRPKMReads``: Per Kilobase per Million mapped reads (RPKM)標準化
• --normalizeTo1x: 按照1x測序深度(reads per genome coverage, RPGC)進行標準化
• --ignoreForNormalization： 指定那些染色體不需要經(jīng)過標準化

如果需要以100為分箱，并且標準化到1x，且僅統(tǒng)計某一條染色體區(qū)域的正鏈，輸出格式為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參數(shù)。

多樣本分析

這部分內(nèi)容主要分析處理組不同重復間的相關程度，會用到multiBamSummary、plotCorrelation和plotPCA三個模塊。。主要目的是看下對照組和處理組中的組間差異和組內(nèi)相似性。

如果上一步把BAM轉換成BW, 那么multiBamSummary可以用multiBigWigSummary替代

# 統(tǒng)計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

根據(jù)下圖不難發(fā)現(xiàn)，組內(nèi)的不同技術重復間差異性小，而組內(nèi)中的兩個生物學重復看起來只能說還行，但是差異還是小于組間的差異。

散點圖

熱圖

但是看主成分分析結果，總感覺哪里不對勁。不過這僅僅是看總體的分布情況，而不是使用差異peak進行主成分分析，也不知道這樣說對不對。

PCA

peak分布可視化

為了統(tǒng)計全基因組范圍的peak在基因特征的分布情況，需要用到computeMatrix計算，用plotHeatmap以熱圖的方式對覆蓋進行可視化，用plotProfile以折線圖的方式展示覆蓋情況。

computeMatrix具有兩個模式:scale-region和reference-point。前者用來信號在一個區(qū)域內(nèi)分布，后者查看信號相對于某一個點的分布情況。

兩者的區(qū)別

無論是那個模式，都有有兩個參數(shù)是必須的，-S是提供bigwig文件，-R是提供基因的注釋信息。

scale-regions模式

computeMatrix scale-regions \ # 選擇模式
       -b 3000 -a 5000 \ # 感興趣的區(qū)域，-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 \ # 感興趣的區(qū)域，-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_1_TSS.gz \ # 輸出為文件用于plotHeatmap, plotProfile
       --outFileSortedRegions 03-read-coverage/ons1regions1_ap2_chip_re1_1_genes.bed

結果可視化

可視化的方法有兩種，一種是輪廓圖，一種是熱圖。兩則都提供了足夠多的參數(shù)對結果進行細節(jié)上的修改。

plotProfile -m matrix1_ap2_chip_rep1_1_TSS.gz \
              -out ExampleProfile1.png \
              --numPlotsPerRow 2 \
              --plotTitle "Test data profile"
plotHeatmap -m matrix1_ap2_chip_rep1_1_TSS.gz \
      -out ExampleHeatmap1.png \

差不多這個樣子

PS：如果是找到的peak可用R包chipSeeker進行可視化。

順便放一下自己的知識星球，如果你覺得我對你有幫助的話。

知識星球