API for intersections

[brentp]

With the current abstractions and API, it's trivial to get a command-line program going that can detect and handle any file type (BAM/CRAM/BED/VCF/etc) and perform intersections:

bedder-rs/src/main.rs

Lines 42 to 45 in fe37df3

	ii.for_each(|intersection| {
	let intersection = intersection.expect("error getting intersection");
	println!("{:?}", intersection);
	});

Users of bedder-rs as a library can then quite simply handle the intersection structs:

bedder-rs/src/intersection.rs

Lines 38 to 51 in fe37df3

	#[derive(Debug)]
	pub struct Intersection {
	/// the Positioned that was intersected
	pub interval: Rc<Position>,
	/// a unique identifier indicating the source of this interval.
	pub id: u32,
	}
	/// An Intersections wraps the base interval and a vector of overlapping intervals.
	#[derive(Debug)]
	pub struct Intersections {
	pub base_interval: Rc<Position>,
	pub overlapping: Vec<Intersection>,
	}

(where Position is an enum of BAM/CRAM/BED/etc).

But how can we make this generally available from the command-line without implementing all of bedtools?

It might be possible to have a very flexible way to do this by allowing user expressions specified from the command-line. Due to the libraries and speeds available, I expect this will be in lua (see https://github.com/brentp/pbr for an example of allowing user expressions), but we can also make it in python with a performance hit.

What will this look like from the **command-line? Here are some on-the-fly musings that I'll revisit and perhaps completely change. Insight/ideas are welcome at any time.

Count intersections

bedder -a my.bed -b other.bed --lua-column "intersections.overlapping.len()"

Note that lua syntax for length of an array is actually #intersections.overlapping (!?) but we'll pretend for now it can use .len().

Filter to intersecting

bedder -a my.bcf -b regions.bed --lua-filter "intersections.overlapping.len() > 0"

Add column to VCF

bedder -a my.vcf -b LCR.bed \
    --info-field '{"Name": "LCR", Description: "how many total unique bases overlapped an LCR", Type: "Integer"}' \
    --lua-column "LCR:intersections.overlapping.union().coverage(intersections.interval)" 

already, this gets quite complicated. But adding a field to a VCF should be one of the more onerous tasks.

With this, perhaps it's good to start by writing union() and coverage() methods in rust since these would be needed and will require more details.

intersect and chop -a intervals to bases covered by -b

This is a common bedtools operation:

bedder -a a.bed -b b.bed \
    --lua-base-filter " ??? "

Now we already have:

1. --lua-filter
2. --lua-column (with multiple uses)
3. --lua-base-filter

[arq5x]

What do you think the performance hit for Python would be?

[brentp]

hard to say. Guessing from e.g. cyvcf2 where python hits every variant, it will depend on how much is done with each interval (variant) after the intersection. Cyvcf2 can iterate over a VCF and print variants with maybe 5% overhead relative to bcftools view (since variant to string is done in C). But once we start messing with INFO fields or format fields the overhead increases.

In all cases there will be a hit for every interval. This is in contrast to the method from @38 (if I understand correctly) that avoids a lot of this overhead.

Once I understand how this intersection stuff will work then I'll prototype a full python API (not just the small expression stuff as in lua here), for now, I hit enough rust issues and design issues.

And, if I implement the lua stuff here, I'll make it such that we can also have python instead.

[38]

Actually I believe one of the most interesting thing in Grass is it actually write python, then transcompile to Rust. There are a binary cache that keeps previously compiled Rust binaries.
The most interesting part is that as user use it more and more, the binary cache will finally becomes a set of binaries that covers all the functionality of bedtools (or even more). So it basically like "let user enumerate their use cases".

Actually I think integrate either python or lua isn't a cheap option at this point, especially most of the filtering function is the hot spot.

[brentp]

I agree there will be large performance hit relative to what Hao has done in grass, but I think it can be manageable, for example we know from cyvcf2 and pysam that we can be within a few percent of C performance despite touching every variant or alignment.

[katiequinlan]

I think performance isn't the top priority, but rather flexibility as that is what bedtools lacks. This comes from the Rust and Python API for new scripts, and I would argue to use Python for filtering on the CLI, simply because it is intuitive for the target community.

Another thing I have been thinking about is how stupidly complex the code in bedtools is for handling the numerous output structures depending on whether it is run without options (just print the intersecting sub-interval), using -wa (print the original A ivl), -wb (add the original B ivl), options like -wo, -c, -loj.

It seems that it might be easier if the output is always a the result of a left outer join (A ivls with no intersections receive NULL defaults for B ivls), and the options add relevant output columns to the output, depending on the options. An exception that ruins this concept is -c, but maybe this strategy could work to simplify the code for most scenarios.

Here is an example

$ cat a.bed
chr1  10  20
chr1  30  40

$ cat b.bed
chr1  12   20

$ bedder -a a.bed -b b.bed
chr1  10  20  chr1 12  20
chr1  30  40  . -1  -1

# filter on -f
$ bedder -a a.bed -b b.bed --filter "intersections.overlapping.a_fraction() > 0.5"
chr1  10  20  chr1 12  20 0.8

# report -f
$ bedder -a a.bed -b b.bed --add "intersections.overlapping.a_fraction() > 0.5"
chr1  10  20  chr1 12  20 0.8
chr1  30  40  . -1  -1 0.0

# report -f and -F
$ bedder -a a.bed -b b.bed \
   --add "intersections.overlapping.a_fraction()" \
   --add "intersections.overlapping.b_fraction()"
chr1  10  20  chr1 12  20 0.8 1.0
chr1  30  40  . -1  -1 0.0 0.0

# report -f, -F, -wo
$ bedder -a a.bed -b b.bed \
   --add "intersections.overlapping.a_fraction()" \
   --add "intersections.overlapping.b_fraction()" \
   --add "intersections.overlapping.overlap()" \   
chr1  10  20  chr1 12  20 0.8 1.0 8
chr1  30  40  . -1  -1 0.0 0.0 0

[arq5x]

Well, that was from @arq5x, but yeah.

[brentp]

@arq5x , I quite like your setup of the problem.
Indeed, we can get quite far with a --filter and a --format (your --add, could also be called --column).

Consider using string formatting (as in luau or python f-strings) to enable:

bedder -a my.bed -b other.bed --format "{a.chrom}\t{a.start}\t{a.stop}\t{intersections.overlapping.len()}"

and for The final example from Aaron above, which was:

 # report -f, -F, -wo
$ bedder -a a.bed -b b.bed \
   --add "intersections.overlapping.a_fraction()" \
   --add "intersections.overlapping.b_fraction()" \
   --add "intersections.overlapping.overlap()" \   
chr1  10  20  chr1 12  20 0.8 1.0 8
chr1  30  40  . -1  -1 0.0 0.0 0

We can have:

bedder -a my.bed -b other.bed \
    --format "{a.chrom}\t{a.start}\t{a.stop}\t{intersections.a_fraction()}\t{intersections.b_fraction()}\t{intersections.overlap()}"

We can have a default so that always "a.chrom\ta.start\ta.end" is already written, in which case, the and --format (equivalent to --add above):

bedder  -a a.bed -b .bed \
    --format "anything that can evaluated and parsed as string..."
   --format "intersections.overlapping.b_fraction()"

Note that the user could also have a single --format with \t.

Then if we want to customize the first 3 columns:

bedder --no-chrom-start-end -a a.bed -b b.bed \
    --format "{a.chrom}\t{clip(a.start, intersections)}\t{clip(a.stop, intersections)}"

the default which outputs chrom, start, end, could be modified with flags -wa, -wb, -wao, -wo`.

I will proceed with this. We will hit limitations or complications when we need to access, for example the INFO field of a VCF. It will need to look something like:

def my_custom_function() {
    columns = []
    for iv in intersections.overlapping {
        # now iv is an enum in rust
        # we need the concrete type in python
        if iv.type() == Type.VCF {
            variant = iv.variant()
             dp = variant.info.get("DP").integer() # again, we convert from enum to python
             columns.append(str(dp))     
        } else if iv.type() == Type.BED {
            # do something
        } else {
            # do something else
        }
    }
    return "\t".join(columns)
}

Presumably, that will be in a script some.py and passed to the command-line and then used as:

bedder --py some.py -a a.bed -b some.vcf.gz \
    --format "my_custom_function()"

[brentp]

We can also use rust templating library (or this one) for fixed fields like chrom, start, stop to reduce overhead

[arq5x]

Awesome, I think this could work well, especially if well documented with lots of examples. What would bs:bases_overlapping() return? The original b intervals or the sub intervals of b that overlap a? For sure we will want a way to return the original intervals and and a way to return the sub intervals that overlap a

[brentp]

Right now bases_overlapping is just a dummy function. It does something, but it's not very useful. I've started to sketch out what the rust API for intersections looks like here:

bedder-rs/src/intersections.rs

Lines 24 to 49 in 6d63a76

	/// IntersectionMode determines what part of overlap is returned.
	/// The f32 is the proportion of overlap required.
	pub enum IntersectionMode {
	// https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html#usage-and-option-summary
	/// Return Chunks of A(B) that overlap B(A)
	Chunks(f32),
	/// Return All of A(B) that overlaps B(A)
	All(f32),
	/// Return A(B) if it does not overlap B(A)
	None(f32),
	}
	impl Intersections {
	pub fn intersections<'a>(
	&'a self,
	a_mode: Option<IntersectionMode>,
	b_mode: Option<IntersectionMode>,
	count: bool,
	) -> Vec<SimpleInterval<'a>> {
	match (a_mode, b_mode) {
	(None, None) => self.a_chunks(None),
	(Some(a_fraction), None) => unimplemented!(),
	(None, Some(b_fraction)) => unimplemented!(),
	(Some(a_fraction), Some(b_fraction)) => unimplemented!(),
	}
	}

it's incomplete, but I'll keep working on it.

[brentp]

I may switch to using bitflags so they are easier to combine.

[brentp]

So, it's more arguments, but I think it makes sense to use enums where arguments are mutally exclusive and bitflags where they can be combined. With that, I have now:

bitflags! {
    /// IntersectionMode indicates requirements for the intersection.
    /// And extra fields that might be reported.
    pub struct IntersectionMode: u8 {
        // https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html#usage-and-option-summary

        /// Default without extra requirements.
        const Empty = 0b00000000;

        /// Return A(B) if it does *not* overlap B(A). Bedtools -v
        const Not = 0b00000001;

        /// Report bases of overlap (-c)
        const Bases = 0b00000010;

        /// Report count of overlaps (-C)
        const Count = 0b00000100;
    }
}

pub enum IntersectionOutput {
    /// Don't report the intersection.
    /// This is commonly used for -b to not report b intervals.
    None,
    /// Report each portion of A that overlaps B
    Part,
    /// Report the whole interval of A that overlaps B
    Whole,
}

pub enum OverlapAmount {
    Bases(u64),
    Fraction(f32),
}

impl Intersections {
    pub fn intersections<'a>(
        &'a self,
        a_mode: IntersectionMode,
        b_mode: IntersectionMode,
        a_output: IntersectionOutput,
        b_output: IntersectionOutput,
        a_requirements: OverlapAmount,
        b_requirements: OverlapAmount,
    ) -> Vec<SimpleInterval<'a>> {
        unimplemented!("Intersections::intersections")
    }
}

So with that, we can, I think replicate most options in bedtools intersect with the exception of strand-based requirements and of -e which is:

Require that the minimum fraction be satisfied for A _OR_ B. In other words, if -e is used with -f 0.90 and -F 0.10 this requires that either 90% of A is covered OR 10% of B is covered. Without -e, both fractions would have to be satisfied.

We could also support this by adding a final argument to intersections that indicates to or the a_requirements and b_requirements as opposed to the default of and.

For strand-based stuff. I propose that can be done by the user (including us when we write the CLI) on the subset of intersections that are returned without regard to strand.

@arq5x , would be good to go over this with you and think about any wholes.

This still ignores, so far, how to actually report intervals, it simply decides how/what is returned from the intersection. But it's a necessary starting point.