如何使用deeptools處理BAM數(shù)據(jù)
總體介紹
deeptools是基于Python開發(fā)的一套工具，用于處理諸如RNA-seq, ChIP-seq, MNase-seq, ATAC-seq等高通量數(shù)據(jù)。工具分為四個模塊
- BAM和bigWig文件處理
- 質(zhì)量控制
- 熱圖和其他描述性作圖
- 其他
當(dāng)然也可以簡單分為兩個部分：數(shù)據(jù)處理和可視化吱涉。
對于deeptools里的任意子命令逞泄，都支持--help看幫助文檔算色，--numberOfProcessors/-p設(shè)置多核處理怎囚，--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，如要重復(fù)請自行下載比對驳概。
BAM轉(zhuǎn)換為bigWig或bedGraph
BAM文件是SAM的二進(jìn)制轉(zhuǎn)換版赤嚼，應(yīng)該都知道。那么bigWig格式是什么顺又？bigWig是wig或bedGraph的二進(jìn)制版更卒，存放區(qū)間的坐標(biāo)軸信息和相關(guān)計分(score)，主要用于在基因組瀏覽器上查看數(shù)據(jù)的連續(xù)密度圖待榔，可用wigToBigWig從wiggle進(jìn)行轉(zhuǎn)換逞壁。
bedGraph和wig格式是什么? USCS的幫助文檔稱這兩個格式數(shù)是已經(jīng)過時的基因組瀏覽器圖形軌展示格式，前者展示稀松型數(shù)據(jù)锐锣，后者展示連續(xù)性數(shù)據(jù)腌闯。目前推薦使用bigWig/bigBed這兩種格式取代前兩者。
為什么要用bigWig? 主要是因為BAM文件比較大雕憔，直接用于展示時對服務(wù)器要求較大姿骏。因此在GEO上僅會提供bw,即bigWig下載，便于下載和查看斤彼。如果真的感興趣分瘦，則可以下載原始數(shù)據(jù)進(jìn)行后續(xù)分析蘸泻。
deeptools提供bamCoverage和bamCompare進(jìn)行格式轉(zhuǎn)換，為了能夠比較不同的樣本嘲玫，需要對先將基因組分成等寬分箱(bin)悦施，統(tǒng)計每個分箱的read數(shù)，最后得到描述性統(tǒng)計值去团。對于兩個樣本抡诞，描述性統(tǒng)計值可以是兩個樣本的比率，或是比率的log2值土陪，或者是差值昼汗。如果是單個樣本，可以用SES方法進(jìn)行標(biāo)準(zhǔn)化鬼雀。
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進(jìn)行可視化顷窒。 這里的參數(shù)僅使用了-e/--extendReads和-bs/--binSize即拓展了原來的read長度，且設(shè)置分箱的大小源哩。其他參數(shù)還有
-
--filterRNAstrand {forward, reverse}: 僅統(tǒng)計指定正鏈或負(fù)鏈 -
--region/-r CHR:START:END: 選取某個區(qū)域統(tǒng)計 -
--smoothLength: 通過使用分箱附近的read對分箱進(jìn)行平滑化
如果為了其他結(jié)果進(jìn)行比較鞋吉，還需要進(jìn)行標(biāo)準(zhǔn)化，deeptools提供了如下參數(shù)：
-
--scaleFactor: 縮放系數(shù) - `--normalizeUsingRPKMReads``: Per Kilobase per Million mapped reads (RPKM)標(biāo)準(zhǔn)化
-
--normalizeTo1x: 按照1x測序深度(reads per genome coverage, RPGC)進(jìn)行標(biāo)準(zhǔn)化 -
--ignoreForNormalization： 指定那些染色體不需要經(jīng)過標(biāo)準(zhǔn)化
如果需要以100為分箱璧疗，并且標(biāo)準(zhǔn)化到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方法進(jìn)行標(biāo)準(zhǔn)化，于是多了--ratio參數(shù)却音。
多樣本分析
這部分內(nèi)容主要分析處理組不同重復(fù)間的相關(guān)程度改抡，會用到multiBamSummary、plotCorrelation和plotPCA三個模塊系瓢。阿纤。主要目的是看下對照組和處理組中的組間差異和組內(nèi)相似性。
如果上一步把BAM轉(zhuǎn)換成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)的不同技術(shù)重復(fù)間差異性小欠拾，而組內(nèi)中的兩個生物學(xué)重復(fù)看起來只能說還行，但是差異還是小于組間的差異骗绕。
但是看主成分分析結(jié)果藐窄，總感覺哪里不對勁。不過這僅僅是看總體的分布情況酬土，而不是使用差異peak進(jìn)行主成分分析荆忍，也不知道這樣說對不對。
peak分布可視化
為了統(tǒng)計全基因組范圍的peak在基因特征的分布情況，需要用到computeMatrix計算刹枉，用plotHeatmap以熱圖的方式對覆蓋進(jìn)行可視化叽唱，用plotProfile以折線圖的方式展示覆蓋情況。
computeMatrix具有兩個模式:scale-region和reference-point微宝。前者用來信號在一個區(qū)域內(nèi)分布棺亭，后者查看信號相對于某一個點的分布情況。
無論是那個模式芥吟，都有有兩個參數(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
結(jié)果可視化
可視化的方法有兩種种樱，一種是輪廓圖蒙袍，一種是熱圖。兩則都提供了足夠多的參數(shù)對結(jié)果進(jìn)行細(xì)節(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進(jìn)行可視化害幅。
順便放一下自己的知識星球，如果你覺得我對你有幫助的話岂昭。
