为什么不得不用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 = " ")
基本安装
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倍以上的时间。