BioJulia / BED.jl

MIT License
6 stars 5 forks source link

Cross compatibility with XAM.jl? #24

Open abhinavsns opened 3 months ago

abhinavsns commented 3 months ago

I am interested to iterate over BAM records .bam that have supplementary index files .bai using XAM.jl. Now I have a bed file of interesting regions, I want to iterate over all bam records for each interval in the .bed file. Is there an example available for this?

jonathanBieler commented 3 months ago

Adapted from XAM docs :

using BED, XAM, GenomicFeatures

# Create an interval collection in memory.
regions = open(BED.Reader, "data.bed") do reader
    IntervalCollection(reader,true)
end

reader = open(BAM.Reader, bamfile)

# Iterate over BAM records.
reader = open(BAM.Reader, bamfile, index = bamfile * ".bai")
for region in regions
    for record in eachoverlap(reader, region)
        # `record` overlaps `region`.
        # ...
    end
end
close(reader)
abhinavsns commented 3 months ago

Thanks! I was wondering if this is the most efficient way? or maybe there is a while !eof(reader) version where I only allocate the Record once and read it in-place? Also it would be great if the develop branch could be released for easier installation and compatibility with latest XAM.jl release.