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

 

2023年03月20日

快速取类似venn的unique和share peaks

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文件

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文件

human人类的基因组

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小鼠的基因组

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，没有意义的瞎玩。

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，不同版本命令略微有差异

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