bedtools的使用技巧(持续更新)

参考链接(注意这个是旧版,但是写的很详细!):https://anjingwd.github.io/AnJingwd.github.io/2017/08/19/bedtools%E4%BD%BF%E7%94%A8%E6%95%99%E7%A8%8B%E8%AF%A6%E8%A7%A3/
参考链接:https://zhuanlan.zhihu.com/p/52322803

技巧一、使用bedtools bamtobed将bam文件转bed文件(基础)

参考:https://uteric.github.io/CHIPSEQ/bamtobed/

bedtools bamtobed -i SRR3022344.sorted.bam >test.bed

第一列:染色体位置
第二列:start
第三列:end
第四列:对应BAM文件的QNAME,包含测序平台,read name等信息
第五列:对应BAM文件的MAPQ,即比对质量
第六列:正负链
注意:
start1和start2起始坐标第一个碱基都为0,所以start=9, end=20表示碱基跨度是从第10位到第20位
染色体用.表示unknown;位置信息用-1表示unknown

image.png

技巧二、使用bedtools intersect计算两个或多个bed中的intersect区域(可接受多个文件类型bed/gff/vcf/bam)

 bedtools intersect [OPTIONS] -a <bed/gff/vcf/bam> -b <bed/gff/vcf/bam>

-wa参数可以报告出原始的在A文件中的feature
-wb参数可以报告出原始的在B文件中的feature
-c参数可以报告出两个文件中的overlap的feature的数量
-wo 返回overlap碱基数
-v 返回非overlap区间
-s 相同链上的feature


#例题请看参考链接!
注意,自己生成测试bed文件,都必须用tab键分割,否则会报错!!

案例一:包含着染色体位置的两个文件,分别记为A文件和B文件。分别来自于不同文件的染色体位置的交集是什么?

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 25

$ bedtools intersect -a A.bed -b B.bed

chr1 15 20

案例二:包含着染色体位置的两个文件,分别记为A文件和B文件。求A文件中哪些染色体位置是与文件B中的染色体位置有overlap.

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 25

$ bedtools intersect -a A.bed -b B.bed -wa

chr1 10 20

案例三:包含着染色体位置的两个文件,分别记为A文件和B文件。求A文件中染色体位置与文件B中染色体位置的交集,以及对应的文件B中的染色体位置.

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 25

$ bedtools intersect -a A.bed -b B.bed -wb

chr1 15 20 chr1 15 25

案例四(经用): 包含着染色体位置的两个文件,分别记为A文件和B文件。求对于A文件的染色体位置是否与文件B中的染色体位置有交集。如果有交集,分别输入A文件的染色体位置和B文件的染色体位置;如果没有交集,输入A文件的染色体位置并以’. -1 -1’补齐文件。

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 25

$ bedtools intersect -a A.bed -b B.bed -loj

chr1 10 20 chr1 15 25

chr1 30 40 . -1 -1

案例五: 包含着染色体位置的两个文件,分别记为A文件和B文件。对于A文件中染色体位置,如果和B文件中染色体位置有overlap,则输出在A文件中染色体位置和在B文件中染色体位置,以及overlap的长度.

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 20

chr1 18 25

$ bedtools intersect -a A.bed -b B.bed -wo

chr1 10 20 chr1 15 20 5

chr1 10 20 chr1 18 25 2

案例六: 包含着染色体位置的两个文件,分别记为A文件和B文件。对于A文件中染色体位置,如果和B文件中染色体位置有overlap,则输出在A文件中染色体位置和在B文件中染色体位置,以及overlap的长度;如果和B文件中染色体位置都没有overlap,则用’. -1-1’补齐文件

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 20

chr1 18 25

$ bedtools intersect -a A.bed -b B.bed -wao

chr1 10 20 chr1 15 20 5

chr1 10 20 chr1 18 25 2

chr1 30 40 . -1 -1

案例七: 包含着染色体位置的两个文件,分别记为A文件和B文件。对于A文件中染色体位置,输出在A文件中染色体位置和有多少B文件染色体位置与之有overlap.

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 20

chr1 18 25

$ bedtools intersect -a A.bed -b B.bed -c

chr1 10 20 2

chr1 30 40 0

案例八(常用): 包含着染色体位置的两个文件,分别记为A文件和B文件。对于A文件中染色体位置,输出在A文件中染色体位置和与B文件染色体位置至少有X%的overlap的记录。

$ cat A.bed

chr1 100 200

$ cat B.bed

chr1 130 201

chr1 180 220

$ bedtools intersect -a A.bed -b B.bed -f 0.50 -wa -wb

chr1 100 200 chr1 130 201

#一个2G的bam文件约需要15-20分钟!!!

注意:当用bedtools intersect 处理大文件时比较耗内存,有效的方法是对A和B文件按照染色体名字(chromosome)和位置(position)排序,重新intersect

bedtools sort [OPTIONS] -i <bed/gff/vcf>
bedtools sort -chrThenSizeA -i test.bed
        -sizeA                  Sort by feature size in ascending order.
        -sizeD                  Sort by feature size in descending order.
        -chrThenSizeA           Sort by chrom (asc), then feature size (asc).
        -chrThenSizeD           Sort by chrom (asc), then feature size (desc).
        -chrThenScoreA          Sort by chrom (asc), then score (asc).
        -chrThenScoreD          Sort by chrom (asc), then score (desc).

技巧三:使用bedtools genomecov染色体和全基因组覆盖度计算(可以用来做深度计算)

单个输入bed文件(-i指定)和genome files;如果输入为bam(-ibam指定)文件,则不需要genome files

bedtools genomecov [OPTIONS] -i <bed/gff/vcf> -g <genome>
 -ibam           The input file is in BAM format
                  Note: BAM _must_ be sorted by position

示例
$ cat ranges-cov-sorted.bed
chr1    4       9
chr1    1       6
chr1    8       19
chr1    25      30
chr2    0       20

$ cat cov.txt  (染色体及每条染色体总碱基数)
chr1    30
chr2    20

bedtools genomecov -i ranges-cov-sorted.bed -g cov.txt
chr1    0       7       30      0.233333 1
chr1    1       20      30      0.666667
chr1    2       3       30      0.1
chr2    1       20      20      1 2
genome  0       7       50      0.14 3
genome  1       40      50      0.8
genome  2       3       50      0.06
#name 覆盖次数 覆盖碱基数 总碱基数   覆盖度
#同时计算单染色体和全基因组覆盖度

如果输入的是-ibam bam 文件,则输出结果为
名称  覆盖次数  该次数下的碱基数 该染色体长度  染色体覆盖度
chr1    0       248951704       248956422       0.999981  
chr1    1       4136    248956422       1.66134e-05
chr1    2       582     248956422       2.33776e-06
chr2    0       242190850       242193529       0.999989
chr2    1       2358    242193529       9.73602e-06
chr2    2       321     242193529       1.32539e-06
chr3    0       198292860       198295559       0.999986
chr3    1       2398    198295559       1.20931e-05
chr3    2       301     198295559       1.51794e-06
chr4    0       190212444       190214555       0.999989
chr4    1       1622    190214555       8.52721e-06
chr4    2       489     190214555       2.57078e-06
chr5    0       181536919       181538259       0.999993
chr5    1       1280    181538259       7.05086e-06
chr5    2       60      181538259       3.30509e-07
chr6    0       170804532       170805979       0.999991
chr6    1       1095    170805979       6.41078e-06
chr6    2       352     170805979       2.06082e-06
chr7    0       159344307       159345973       0.99999


 chr1  248951704+4136+582=248956422 


2G的bam文件大概的时间消耗
real    3m29.649s
user    3m19.071s
sys     0m9.387s
genomecov也会对全基因组的按照不同深度做统计!!(有用啊)

-u Write the original A entry once if any overlaps found in B
注意-u参数,1、使用后相当于进入-wa模式,2、若A与B有相同重复,则会去重

image.png

技巧四、使用bedtools coverage计算染色体给定区间的深度和覆盖度,输入文件可以是bam

bedtools coverage [OPTIONS] -a <bed/gff/vcf> -b <bed/gff/vcf>
其中-a指定interval文件,即你想看的染色体区间  -b指定的是你比对的结果或bed等文件


示例
$ cat A.bed
chr1  0   100
chr1  100 200
chr2  0   100

$ cat B.bed
chr1  10  20
chr1  20  30
chr1  30  40
chr1  100 200

$ bedtools coverage -a A.bed -b B.bed
A的名称 起始 终止 B在A中匹配次数 匹配的总长度 该区域总长度  比例
chr1    0   100    3    30    100    0.3000000
chr1  100 200  1  100 100 1.0000000
chr2  0   100  0  0   100 0.0000000
结果解释:前三列是interval文件的信息,后四列为统计信息
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,711评论 5 468
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 83,932评论 2 376
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,770评论 0 330
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,799评论 1 271
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,697评论 5 359
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,069评论 1 276
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,535评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,200评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,353评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,290评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,331评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,020评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,610评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,694评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,927评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,330评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,904评论 2 341