bedtools核心用法详解

为什么不得不用bedtools?

  • 速度,当数据到达百万级以上,R和C的速度差别就非常明显了
  • 专业,但凡涉及到region、peak的处理,bedtools都可以胜任

 

目录

bed文件基本处理:导出,排序

用一个query region去overlap另一个database region,并取得属性

 

bed文件基本处理:导出,排序

把R里面的region导出到bed文件

chr <- unlist(lapply(tmp$frag2, function(x) {
    strsplit(x, split = ":|-")[[1]][1]
}))

start <- unlist(lapply(tmp$frag2, function(x) {
    strsplit(x, split = ":|-")[[1]][2]
}))

end <- unlist(lapply(tmp$frag2, function(x) {
    strsplit(x, split = ":|-")[[1]][3]
}))

tmp.bed <- data.frame(chr=chr, start=start, end=end, rsid=current.all.LD.all$RS_Number)

write.table(tmp.bed, file = "cHi-C/tmp.input.bed", row.names = F, col.names = F, quote = F, sep = "\t")

  

基本安装

conda install -c bioconda bedtools

  

排序,构建一个chr.list文件

chr1
chr2
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chr10
chr11
chr12
chr13
chr14
chr15
chr16
chr17
chr18
chr19
chr20
chr21
chr22
chrX
chrY

  

然后sort,不同版本命令略微有差异

cat capture_HiC.curated.bed | bedtools sort -faidx chr.list > capture_HiC.curated.sorted.bed
cat tmp.input.bed | bedtools sort -faidx chr.list > tmp.input.sorted.bed

  

用一个query region去overlap另一个database region,并取得属性

bedtools intersect -a tmp.input.sorted.bed -b capture_HiC.curated.sorted.bed -wo > tmp.output.bed

  

这部分比较细节,需要仔细参考教程:https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html

-a query bed,比如SNP的位置

-b database bed,比如capture Hi-c的位置

-wo 如果有overlap,则原样输出-a和-b的文件信息

 

如果用R的条件规则去判断,估计要花10倍以上的时间。

 

  

 

posted @ 2021-03-03 16:09  Life·Intelligence  阅读(66)  评论(0编辑  收藏
TOP