搜索
查看: 15335|回复: 6

生信人必须会的软件之bedtools

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2016-9-8 20:11:16 | 显示全部楼层 |阅读模式
bedtools等工具号称是可以代替普通的生物信息学数据处理工程师的!我这里用一个专题来讲解它的用法,其实它能实现的需求,我们写脚本都是可以做的,而且我强烈建议正在学编程的小朋友模仿它的各种功能来增强自己的脚本功力。
好,言归正传,bedtools是非常多的工具的合集,有瑞士军刀的美誉。直接下载二进制版本软件就可以调用全路径来使用,或者把整个文件夹添加到环境变量。
我们先看第一个功能,把alignment的结果文件转为bedgraph格式文件。
参考:http://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html
要实现这个功能,需要用bedtools的genomecov小命令, 有两种方法可以调用:
bedtools genomecov [OPTIONS] [-i|-ibam] -g (iff. -i)
genomeCoverageBed  [OPTIONS] [-i|-ibam] -g (iff. -i)

这个命令本身并不是设计来做格式转换的,bam2bedgraph也只是其中的一个小功能而已,需要加上-bg参数,就可以Report depth in BedGraph format. For details, see: http://genome.ucsc.edu/goldenPath/help/bedgraph.html   
大家观摩我下面给出的测试例子,就明白该功能如何使用了
[Shell] 纯文本查看 复制代码
bedtools genomecov  -bg -i E001-H3K4me1.tagAlign -g mygenome.txt >E001-H3K4me1.bedGraph
bedtools genomecov  -bg -i E001-Input.tagAlign -g mygenome.txt >E001-Input.bedGraph
nohup bedtools genomecov  -bg -ibam BAF180_CT10.unique.sorted.bam >BAF180_CT10.bedGraph &
nohup bedtools genomecov  -bg -ibam BAF180_CT22AM.unique.sorted.bam >BAF180_CT22AM.bedGraph &
nohup bedtools genomecov  -bg -ibam BAF180_CT22.unique.sorted.bam >BAF180_CT22.bedGraph &
nohup bedtools genomecov  -bg -ibam inputCT10sonication.unique.sorted.bam >inputCT10sonication.bedGraph &

首先alignment的文件必须是sort的,然后如果是bed格式的比对文件,用-i 参数来指定输入文件,需要加入参考基因组的染色体大小记录文件(mygenome.txt ),如果是bam格式的比对文件,用-ibam指定输入文件,而且不需要参考基因组的染色体大小记录文件。

然后看看第二个功能,对RNA-seq的比对文件中的比对到各个基因的reads进行计数。
参考:http://www.bio-info-trainee.com/745.html
实现这个功能,用的bedtools的multicov 这个小命令

命令很简单的

bedtools multicov [OPTIONS] -bams BAM1 BAM2 BAM3 … BAMn -bed  <BED/GFF/VCF>

By default, multicov will report the count of alignments in each input BAM file that overlap.

例子:

bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed ivls-of-interest.bed

ivls-of-interest.bed这个文件是必须的,可能需要自己制作,其实用gtf文件也可以的

chr1 0   10000   ivl1

chr1 10000   20000   ivl2

chr1 20000   30000   ivl3

chr1 30000   40000   ivl4

输出结果前三列是坐标,第四列是基因名,跟我们的bed文件一样,只是最后三列是三个样本的计数,是添加上来的!

chr1 0       10000   ivl1    100 2234    0

chr1 10000   20000   ivl2    123 3245    1000

chr1 20000   30000   ivl3    213 2332    2034

chr1 30000   40000   ivl4    335 7654    0

可以看到,它实现的需求,跟htseq这个软件差不多。各种区别,大家可以自己取探索。

接着第三个功能,根据坐标区域来从基因组里面提取fasta序列
参考:#         BED/GFF/VCF +reference --> fasta
#        http://bedtools.readthedocs.io/e ... tools/bamtobed.html
#        http://bedtools.readthedocs.io/e ... tools/getfasta.html
命令也很简单:
[Shell] 纯文本查看 复制代码
bedtools getfasta -fi ~/biosoft/bowtie/hg19_index/hg19.fa  -bed ../macs14_results/highQuality_summits.bed  -fo highQuality.fa
bedtools getfasta -fi ~/biosoft/bowtie/hg19_index/hg19.fa  -bed ../macs14_results/highQuality_peaks.bed  -fo highQuality.fa

我的例子脚本里面用的是bed格式来记录坐标区域,参考基因组用-fi参数指定具体位置,输出的fasta序列文件用-fo参数指定







上一篇:annovar对人类基因组variants注释
下一篇:生信人必须会的软件之samtools
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

1

主题

5

帖子

117

积分

版主

Rank: 7Rank: 7Rank: 7

积分
117
发表于 2016-12-16 17:30:16 | 显示全部楼层


很久之前的之前,一直在学习如何用WES外显子数据call CNV,比较之后,最终选择用conifer与cnvkit,但是其基本原理都来自于以下的bedtools结合DNAcopy,既然这个帖子是说bedtools,我也来分享一下以前从bedtools作者那里学到的一点关于call CNV的分析思路,供大家学习。谢谢

###############################################################
# To make a genome file (for bed tools) using reference genome
###############################################################

awk -v OFS='\t' {'print $1,$2'} hg19.fa.fai > hg19_genomeFile.txt

###############################################################
output hg19.txt
###############################################################
#chr1  195471971
#chr10  130694993
#chr11  122082543
#chr12  120129022
#chr13  120421639


########################################################
# Create a BED file of 5kb windows with 2.5kb overlap
# tiling build 37 (hg19) of the human genome
########################################################
bedtools makewindows -g hg19.txt -w 5000 -s 2500 > hg19.w5k.s2.5k.bedg



########################################################
# Count the number of reads aligned to each overlapping
# 5kb genomic window for each sample.
########################################################
bedtools intersect -a hg19.w5k.s2.5k.bedg \
                   -b /scratch/ucgd/lustre/Projects/WashU/15-05-21_Marth_Quinlan_Breast_Cancer/alignment/B0.bam \
                   -sorted -c -g hg19.txt \
                   > B0.w5k.s2.5k.conc.bedg


########################################################
# restrict depth calls to chrom 1 through X
########################################################
awk 'length($1) <=2' B0.w5k.s2.5k.conc.bedg | grep -v "^MT" | grep -v "^Y" > B0.w5k.s2.5k.conc.chr1-chrX.bedg


########################################################
# Prefix chromosomes with chr
########################################################
awk '{print "chr"$0}' B3.w5k.s2.5k.conc.chr1-chrX.bedg > B3.w5k.s2.5k.conc.chr1-chrX.chromlabels.bedg


########################################################
# Normalize each samples coverage by GC content
# Uses the codachrom package:
#  https://github.com/arq5x/codachrom
########################################################
python codachrom/windowed_coverage.py -w 5000 -s 2500 -c B0.w5k.s2.5k.conc.chr1-chrX.chromlabels.bedg -f hg19.fa


########################################################
# Compute log2 ratios of one sample's (-1) normalized
# coverage versus another's (-2).
########################################################
python codachrom/log_ratio.py -1 B1.w5k.s2.5k.conc.chr1-chrX.chromlabels.bedg.depth.txt -2 B0.w5k.s2.5k.conc.chr1-chrX.chromlabels.bedg.depth.txt > B1vB0.log2.bedg &


# segment the log2 ratios using DNAcopy
# see details in seg.R script

########################################################
# convert DNAcopy output to bedgraph
########################################################
awk '{OFS="\t"; print $3,$4,$5,$7,$6,$1}' segs.b1b0.txt > segs.b1b0.bedg


########################################################
# Identify segments with a mean log2 ratio of > 0.4
# or less than -0.4 and that are at least 15kb in size
########################################################
awk '($4>=0.4 || $4<=-0.4) && ($5 >=5)' segs.b1b0.bedg > segs.b1b0.cnvs.bedg


########################################################
# merge segments that are within 250Kb of one another.
########################################################
sort -k1,1 -k2,2n segs.b3b2.cnvs.bedg | bedtools merge -i - -d 250000 -c 4,5 -o mean,sum > segs.b3b2.cnvs.merged.bedg


########################################################
# get protein coding, known GENCODE genes in
# BED format (downloaded from UCSC)
########################################################
grep protein_coding gencode.bed | grep KNOWN | sort -k1,1 -k2,2n > gencode.prot.known.sorted.bed


########################################################
# collect a gene list for each putative CNV
########################################################
bedtools intersect -a segs.b2b1.cnvs.merged.bedg -b gencode.prot.known.sorted.bed -wa -wb | bedtools groupby -i - -g 1-5 -c 9 -o distinct > segs.b2b1.cnvs.merged.genelist.bedgraph


########################################################
# bgzip and tabix for IGV
########################################################
bgzip segs.b1b0.cnvs.merged.genelist.bedgraph

tabix -p bed segs.b1b0.cnvs.merged.genelist.bedgraph.gz




本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 2 反对 0

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2017-3-22 20:55:24 | 显示全部楼层
1、intersectBed
用来求两个BED或者BAM文件中的overlap,overlap可以进行自定义是整个genome features的overlap还是局部。 加-wa参数可以报告出原始的在A文件中的feature,加-wb参数可以报告出原始的在B文件中的feature, 加-c参数可以报告出两个文件中的overlap的feature的数量,参数-s可以得到忽略strand的overlap,具体案例如下:

案例一:包含着染色体位置的两个文件,分别记为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
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 1 反对 0

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2016-10-13 09:03:42 | 显示全部楼层
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

0

主题

3

帖子

55

积分

注册会员

Rank: 2

积分
55
发表于 2016-10-13 09:28:55 | 显示全部楼层
bedtools ——实用功能之intersect  

功能:取两个间或一个与多个bed文件的交集

使用:bedtools intersect [OPTIONS] -a <FILE> \
                                                   -b  <FILE1,FILE3,FILE3,...,FILEN>

-a 相当于 blast中的query sequence
-b 相当于寻找与query sequence 有交集的subject sequences,2.21.0版本的bedtools 的 subjectsequences 可以为多个

这个功能非常完整,大部分条件都能找到相应的参数来控制:
1.BED/GFF/VCF/和BAM文件都可以做为input file,而且加上 -bed 参数时,输出文件就变成最易懂的bed格式
2. 参数 -s 和 -S 可以控制正负链
所以这个功能的参数也特别多,详细见:http://bedtools.readthedocs.io/e ... ools/intersect.html
当大家遇到这个取序列交集的问题时,像我这样的菜鸟推荐使用 bedtools intersect实现第一次发帖,也不知道有没有描述清楚,见谅哈
回复 支持 反对

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2017-3-22 20:54:39 | 显示全部楼层
BEDTools是可用于genomic features的比较,相关操作及进行注释的工具。而genomic features通常使用Browser Extensible Data (BED) 或者 General Feature Format (GFF)文件表示,用UCSC Genome Browser进行可视化比较。bedtools总共有二三十个工具/命令来处理基因组数据。

比较典型而且常用的功能举例如下:

格式转换,bam转bed(bamToBed),bed转其他格式(bedToBam,bedToIgv);
对基因组坐标的逻辑运算,包括:交集(intersectBed,windowBed),”邻集“(closestBed),补集(complementBed),并集(mergeBed),差集(subtractBed);
计算覆盖度(coverage)(coverageBed,genomeCoverageBed);
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2017-3-22 20:55:51 | 显示全部楼层
3、bedtools multicov.
利用比对的bam文件对每个基因进行read数目计数。
首先注意它需要sort的bam文件及bam的index
bedtools multicov [OPTIONS] -bams BAM1 BAM2 BAM3 … BAMn -bed  <BED/GFF/VCF>
例子:
bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed interest.bed
interest.bed这个文件是必须的,可能需要自己制作,其实用gtf文件也可以的
chr1 0   10000   gene1
chr1 10000   20000   gene2
chr1 20000   30000   gene3
chr1 30000   40000   gene4
输出结果前三列是坐标,第四列是基因名,跟我们的bed文件一样,只是最后三列是三个样本的计数,是添加上来的!
chr1 0       10000   gene1    100 2234    0
chr1 10000   20000  gene2    123 3245    1000
chr1 20000   30000   gene3    213 2332    2034
chr1 30000   40000   gene4    335 7654    0
得到的这个矩阵就可以去用DESeq包来进行差异分析啦!
4.bedtools merge
用于合并位于同一个bed/gff/vcf 文件中的重叠区域。
Bedtools merge  [OPTION] –i <bed/gff/vcf>
-s 必须相同(正负)链的区域才合并(软件默认不考虑正负链特征)
-n 报告合并的区域数量,没有被合并则1
-d 两个独立区域间距小于(等于)该值时将被合并为一个区域
-nms 报告被合并区域的名称
-scores 报告几个被合并特征区域的scores

5、其他
1)pairToPair
比较BEDPE文件搜索overlaps, 类似于pairToBed。
2)bamToBed
将BAM文件转换为BED文件或者BEDPE文件。
bamToBed -i reads.bam
3)windowBed类似于intersectBed, 但是可以指定一个数字,让A中的genome feature增加上下游去和B中的genome features进行overlap。默认情况这个值为1000,可以使用-w加定义,可以用-l指定是上游,用-r指定下游
windowBed -a A.bed -b B.bed -w 5000
windowBed -a A.bed -b B.bed -l 200 -r 20000
4)subtractBed在A中去除掉B中有的genome features
5)coverageBed可以计算深度和覆盖度。如计算基因组任意1Kb的测序read的覆盖度
6)genomeCoverageBed。可以计算给定bam文件在基因组上的覆盖度及每个碱基的深度。
三、软件相关论文:
Quinlan, A.R. & Hall, I.M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841-842 (2010).
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-8-20 03:21 , Processed in 0.045522 second(s), 31 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.