Efficient intersection of BED and VCF

[apraga]

Hi,

When trying to reproduce the equivalent of bedtools intersect for taking the intersect of a VCF file with a BED file, my naive implementation was to loop over all BED intervals and select matching VCF entries (code below).
This result in a rather slow code while bedtools interesct is very fast. Their algorithm seems to use bins.

Is there a better way to rewrite this using your package to improve the performance ?
Thanks !

fn intersect(vcf_gz: &str, bed_f: &str) -> Result<(), Box<dyn Error>>{
    let mut reader_bed = File::open(bed_f).map(BufReader::new).map(bed::Reader::new)?;

    let mut reader_vcf = File::open(vcf_gz)
        .map(bgzf::Reader::new)
        .map(vcf::Reader::new)?;
    let vcf_tbi = String::from(vcf_gz) + ".tbi";
    let header_vcf = reader_vcf.read_header()?;
    let index = tabix::read(vcf_tbi)?;

    // For each interval in the capture kit
    for result in reader_bed.records::<3>() {
        let record = result?;
        let chr = record.reference_sequence_name().replace("chr", ""); 
        let region = Region::new(chr,
            record.start_position()..=record.end_position());
        let query_vcf = reader_vcf.query(&header_vcf, &index, &region)?;

        // Get matching records
        for result_vcf in query_vcf {
            let record = result_vcf?;
        }
    }
    Ok::<_, Box<dyn std::error::Error>>(())
}