bedtools核心用法详解 | Venn | unique | share

 

2023年03月20日

快速取类似venn的unique和share peaks

?

1 2 3 4

module load bedtools2-2.27.1-gcc-5.4.0-iouh4nk bedtools intersect -a HDAC1.merged.narrowPeak -b HDAC2.merged.narrowPeak -wo > HDAC1_2.shared.bed bedtools intersect -a HDAC1.merged.narrowPeak -b HDAC2.merged.narrowPeak -wo -v > HDAC1.uniq.bed bedtools intersect -b HDAC1.merged.narrowPeak -a HDAC2.merged.narrowPeak -wo -v > HDAC2.uniq.bed

 

---

 

为什么不得不用bedtools？

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

 

目录

bed文件基本处理：导出，排序

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

 

bed文件基本处理：导出，排序

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

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

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")

　　

基本安装

?

1	conda install -c bioconda bedtools

　　

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

human人类的基因组

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

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

　　

mouse小鼠的基因组

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

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

　　

手动导出scATAC的peak，这里的不是fragment，而是peak，没有意义的瞎玩。

?

1 2 3 4 5 6

cat *fragments.tsv > all.fragments.tsv only replace the first one sed -e 's/\:/\t/' -e 's/-/\t/' all.fragments.tsv > fragments.all.tsv cat fragments.all.tsv | bedtools sort -faidx mouse.chr.list > fragments.sort.tsv bgzip fragments.sort.tsv tabix -p bed fragments.sort.tsv.gz

　　

　　

 

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

?

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

?

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

　　

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

?

1	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倍以上的时间。