做完比對(duì)后，最直接的一個(gè)問(wèn)題就是想知道：你的reads比對(duì)上了哪些區(qū)域沒(méi)比對(duì)上哪些區(qū)域，有多少reads比對(duì)上特定的區(qū)域(該區(qū)域比對(duì)深度)。傳統(tǒng)手藝一般可以使用bedtools cov或者samtools depth來(lái)直接實(shí)現(xiàn)，但是受限其速度，這些工具我個(gè)人用起來(lái)體驗(yàn)不太好，特別是對(duì)文件大小很大的BAM file。今天就給大家推薦一款我覺(jué)得非常好用，速度很快的比對(duì)文件BAM/CAM深度檢測(cè)好用工具，Mosdepth。
工具簡(jiǎn)介
Mosdepth是一種新的命令行工具，用于快速計(jì)算全基因組測(cè)序覆蓋率。它測(cè)量BAM或CRAM文件在基因組中每個(gè)核苷酸位置或基因組區(qū)域的深度?？梢詫⒒蚪M區(qū)域指定為被BED文件覆蓋的指定區(qū)域，或者指定為CNV calling所需的固定大小窗口。Mosdepth使用一種計(jì)算效率高的簡(jiǎn)單算法，使其能夠快速生成覆蓋摘要。Mosdepth比現(xiàn)有工具更快，并且提供了所產(chǎn)生的覆蓋剖面類型的靈活性。
工具特點(diǎn)
總的來(lái)說(shuō)，Mosdepth有以下幾項(xiàng)比較顯著的特點(diǎn):
- 對(duì)比bedtools, samtools 速度快
- 可以根據(jù)給定窗口大小來(lái)計(jì)算平均每個(gè)窗口深度，用于CNV calling
- 能根據(jù)BED文件給出其覆蓋的指定區(qū)域的深度
- 量化輸出，它合并相鄰的堿基，只要它們都同時(shí)落入相同的覆蓋深度的區(qū)域，例如（10X-20X)。
- 能給出每個(gè)染色體和全基因組在特定閾值或以上覆蓋的堿基比例分布。
- 能夠設(shè)定閥值輸出，用于表示在給定閾值處覆蓋每個(gè)區(qū)域中有多少個(gè)堿基。
- 能給出每條染色體和每條染色體特定區(qū)域內(nèi)平均深度的總結(jié)。
工作原理
當(dāng)遇到每條染色體時(shí)，mosdepth會(huì)在染色體的長(zhǎng)度上創(chuàng)建一個(gè)數(shù)組。對(duì)于遇到的每個(gè)起始位置，它會(huì)增加數(shù)組位置的值。對(duì)于每個(gè)停止的位置，它會(huì)減少該位置。由此，特定位置處的深度是其前面的所有陣列位置的累積和（在BEDTools中使用類似的算法，其中單獨(dú)跟蹤開始和停止）。 mosdepth避免重復(fù)計(jì)算重疊的配偶對(duì)，并使用CIGAR操作跟蹤每個(gè)讀取的每個(gè)比對(duì)上的部分。由于這種數(shù)據(jù)結(jié)構(gòu)，可以在不顯著增加運(yùn)行時(shí)間的情況下完成覆蓋分布計(jì)算。下圖顯示了這個(gè)概念：
這個(gè)數(shù)組計(jì)算非?？?。因?yàn)槠錄](méi)有額外的分配或跟蹤對(duì)象，它在概念上也很簡(jiǎn)單。由于這些原因，它比samtools depth更快。
工具安裝
目前最新版是v 0.2.6,安裝方法有很多可以直接下載編譯好的版本:
###下載地址
https://github.com/brentp/mosdepth/releases
或者直接通過(guò)conda來(lái)安裝
conda install mosdepth
使用說(shuō)明
mosdepth 0.2.3
Usage: mosdepth [options] <prefix> <BAM-or-CRAM>
Common Options:
-t --threads <threads> number of BAM decompression threads. (use 4 or fewer) [default: 0]
-c --chrom <chrom> chromosome to restrict depth calculation.
-b --by <bed|window> optional BED file or (integer) window-sizes.
-n --no-per-base dont output per-base depth. skipping this output will speed execution
substantially. prefer quantized or thresholded values if possible.
-f --fasta <fasta> fasta file for use with CRAM files.
Other options:
-F --flag <FLAG> exclude reads with any of the bits in FLAG set [default: 1796]
-i --include-flag <FLAG> only include reads with any of the bits in FLAG set. default is unset. [default: 0]
-x --fast-mode dont look at internal cigar operations or correct mate overlaps (recommended for most use-cases).
-q --quantize <segments> write quantized output see docs for description.
-Q --mapq <mapq> mapping quality threshold [default: 0]
-T --thresholds <thresholds> for each interval in --by, write number of bases covered by at
least threshold bases. Specify multiple integer values separated
by ','.
-R --read-groups <string> only calculate depth for these comma-separated read groups IDs.
-h --help show help
簡(jiǎn)單說(shuō)說(shuō)我平時(shí)常用到的參數(shù)：
-
-t使用線程的數(shù)目，這里推薦少于4個(gè)CPU，因?yàn)槲臋n有說(shuō)到多于4個(gè)線程，速度不會(huì)有提升 -
-c指定具體哪個(gè)染色體去計(jì)算深度 -
-b能根據(jù)BED文件給出其覆蓋的指定區(qū)域的深度,比如你給出的BED文件含有一個(gè)物種對(duì)應(yīng)基因的位置，就能計(jì)算其基因的覆蓋度。 -
-F在計(jì)算覆蓋度時(shí)，排除BAM中含有該FLAG的reads。默認(rèn)是1796，代表著:read unmapped (0x4), not primary alignment (0x100)， read fails platform/vendor quality checks (0x200)， read is PCR or optical duplicate (0x400)。一般設(shè)定為默認(rèn)就好了。 -
-Q比對(duì)最低質(zhì)量的要求，默認(rèn)是0，一般可以按照自己的需要進(jìn)行調(diào)整為5或者10。 -
-x快速模式，推薦給大型的BAM文件使用 -
-T按照閥值（比如1X,2X,10X）輸出達(dá)到指定閥值的核苷酸數(shù)目。
這里可以根據(jù)自己的需要使用對(duì)應(yīng)的參數(shù)。下面給出一些常用的命令行運(yùn)行例子以供大家參考使用:
根據(jù)bed文件來(lái)計(jì)算指定區(qū)間的覆蓋度:
mosdepth -b capture.bed sample-output sample.exome.bam
進(jìn)一步使用-T根據(jù)指定閥值進(jìn)行計(jì)算:
mosdepth --by exons.bed --thresholds 1,10,20,30 $prefix $bam
輸出效果如下:
#chrom start end region 1X 10X 20X 30X
1 11869 12227 ENSE00002234944 358 157 110 0
1 11874 12227 ENSE00002269724 353 127 10 0
1 12010 12057 ENSE00001948541 47 8 0 0
1 12613 12721 ENSE00003582793 108 0 0 0
該工具我個(gè)人體驗(yàn)非常好，又快又好用，輸出結(jié)果簡(jiǎn)單易懂。如果有需要計(jì)算比對(duì)覆蓋需求的小伙伴可以試一試哦。最后給大家附上其github鏈接以便大家進(jìn)一步研究使用該工具:https://github.com/brentp/mosdepth
