比对文件BAM/CAM深度检测好用工具:Mosdepth

做完比对后,最直接的一个问题就是想知道:你的reads比对上了哪些区域没比对上哪些区域,有多少reads比对上特定的区域(该区域比对深度)。传统手艺一般可以使用bedtools cov或者samtools depth来直接实现,但是受限其速度,这些工具我个人用起来体验不太好,特别是对文件大小很大的BAM file。今天就给大家推荐一款我觉得非常好用,速度很快的比对文件BAM/CAM深度检测好用工具,Mosdepth

工具简介

Mosdepth是一种新的命令行工具,用于快速计算全基因组测序覆盖率。它测量BAM或CRAM文件在基因组中每个核苷酸位置或基因组区域的深度。可以将基因组区域指定为被BED文件覆盖的指定区域,或者指定为CNV calling所需的固定大小窗口。Mosdepth使用一种计算效率高的简单算法,使其能够快速生成覆盖摘要。Mosdepth比现有工具更快,并且提供了所产生的覆盖剖面类型的灵活性。

工具特点

总的来说,Mosdepth有以下几项比较显著的特点:

  • 对比bedtools, samtools 速度快
  • 可以根据给定窗口大小来计算平均每个窗口深度,用于CNV calling
  • 能根据BED文件给出其覆盖的指定区域的深度
  • 量化输出,它合并相邻的碱基,只要它们都同时落入相同的覆盖深度的区域,例如(10X-20X)。
  • 能给出每个染色体和全基因组在特定阈值或以上覆盖的碱基比例分布。
  • 能够设定阀值输出,用于表示在给定阈值处覆盖每个区域中有多少个碱基。
  • 能给出每条染色体和每条染色体特定区域内平均深度的总结。

工作原理

当遇到每条染色体时,mosdepth会在染色体的长度上创建一个数组。对于遇到的每个起始位置,它会增加数组位置的值。对于每个停止的位置,它会减少该位置。由此,特定位置处的深度是其前面的所有阵列位置的累积和(在BEDTools中使用类似的算法,其中单独跟踪开始和停止)。 mosdepth避免重复计算重叠的配偶对,并使用CIGAR操作跟踪每个读取的每个比对上的部分。由于这种数据结构,可以在不显着增加运行时间的情况下完成覆盖分布计算。下图显示了这个概念:

这个数组计算非常快。因为其没有额外的分配或跟踪对象,它在概念上也很简单。由于这些原因,它比samtools depth更快。

工具安装

目前最新版是v 0.2.6,安装方法有很多可以直接下载编译好的版本:

###下载地址
https://github.com/brentp/mosdepth/releases

或者直接通过conda来安装

conda install mosdepth

使用说明

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

简单说说我平时常用到的参数:

  • -t使用线程的数目,这里推荐少于4个CPU,因为文档有说到多于4个线程,速度不会有提升
  • -c指定具体哪个染色体去计算深度
  • -b能根据BED文件给出其覆盖的指定区域的深度,比如你给出的BED文件含有一个物种对应基因的位置,就能计算其基因的覆盖度。
  • -F在计算覆盖度时,排除BAM中含有该FLAG的reads。默认是1796,代表着:read unmapped (0x4), not primary alignment (0x100), read fails platform/vendor quality checks (0x200), read is PCR or optical duplicate (0x400)。一般设定为默认就好了。
  • -Q 比对最低质量的要求,默认是0,一般可以按照自己的需要进行调整为5或者10。
  • -x 快速模式,推荐给大型的BAM文件使用
  • -T 按照阀值(比如1X,2X,10X)输出达到指定阀值的核苷酸数目。

这里可以根据自己的需要使用对应的参数。下面给出一些常用的命令行运行例子以供大家参考使用:

根据bed文件来计算指定区间的覆盖度:

mosdepth -b capture.bed sample-output sample.exome.bam

进一步使用-T根据指定阀值进行计算:

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

该工具我个人体验非常好,又快又好用,输出结果简单易懂。如果有需要计算比对覆盖需求的小伙伴可以试一试哦。最后给大家附上其github链接以便大家进一步研究使用该工具:https://github.com/brentp/mosdepth

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 217,657评论 6 505
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,889评论 3 394
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 164,057评论 0 354
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,509评论 1 293
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,562评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,443评论 1 302
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,251评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,129评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,561评论 1 314
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,779评论 3 335
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,902评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,621评论 5 345
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,220评论 3 328
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,838评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,971评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,025评论 2 370
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,843评论 2 354

推荐阅读更多精彩内容

  • 别让抠门毁了你,这7件事必须大手笔 1在打破僵局上必须大手笔 2在打通关系上必须大手笔 3在学习提升上必须大手笔 ...
    阿甘1972阅读 224评论 0 0
  • 文/伊米crystal 昨天很忙,真的没有时间构思文章,抱歉了,看来日更目标完不成了。 趁娃睡下,刚刚又回听了昨天...
    伊米crystal阅读 672评论 5 19
  • 今天是什么日子 起床:4:48 就寝:20:46 天气:晴 心情:喜悦 纪念日:日更第73天,新工作入职第一天 叫...
    江南烟雨小青年阅读 89评论 0 0
  • 日暮斜阳,照目瑙乡晚,正南国春暖, 群芳斗艳,风筝比高。 傣乡竹翠,舞孔雀之屏,恰春风起时,枯叶纷飞,桂花香来。
    黔来客阅读 198评论 0 0
  • 互联网发展至今,与传统行业的结合越来越紧密,很多传统行业的从业者如医生、教师、司机成为互联网产品的用户,而这类用户...
    香9醋阅读 415评论 0 2