抱歉,您的浏览器无法访问本站
本页面需要浏览器支持(启用)JavaScript
了解详情 >

前言

自己写的好几种算法企图实现bedtools的功能,虽然julia性能足够好,但都难以在效率上达到bedtools的性能,于是最后只能借助轮子了。

代码

using CSV
using DataFrames
using Tables
using GenomicFeatures

function fs2NameTuple(fs::String)
query = DataFrame(CSV.File(fs, delim="\t", header=0))
if size(query)[2] == 5
query = rename!(query, :Column1 => :Chromosome, :Column2 => :Start, :Column3 => :End, :Column4 => :Name, :Column5 => :Score)
# query = sort!(query, [:Chromosome,:Start])
query = Tables.rowtable(query)
result = IntervalCollection([Interval(qy.Chromosome, qy.Start, qy.End, '?', qy.Name) for qy in query], true)
result
elseif size(query)[2] == 4
query = rename!(query, :Column1 => :Chromosome, :Column2 => :Start, :Column3 => :End, :Column4 => :Name)
# query = sort!(query, [:Chromosome,:Start])
query = Tables.rowtable(query)
result = IntervalCollection([Interval(qy.Chromosome, qy.Start, qy.End, '?', qy.Name) for qy in query], true)
result
else
println("please confirm your bed files")
end
end

function getInterval(A::IntervalCollection, B::IntervalCollection)
for i in eachoverlap(A, B)
println(i[1].seqname, "\t", i[1].first, "\t", i[1].last, "\t", i[1].metadata, "\t", i[2].seqname, "\t", i[2].first, "\t", i[2].last, "\t", i[2].metadata)
end
end
getInterval(fs2NameTuple("test1.bed"),fs2NameTuple("test2.bed"))

原理

首先通过CSV将bed文件读入内存,然后将dataFrame转化为nametuple,最后在转化为intervalCollection
通过内置的overlap功能可以几乎瞬间获取所有的交集,性能已经达到了c语言的程度。

评论