Add function for exploding intervals to loci - #789
Conversation
ch-kr
left a comment
There was a problem hiding this comment.
thanks for adding! a few suggestions
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
ch-kr
left a comment
There was a problem hiding this comment.
a few more minor changes -- sorry about the back and forth on the MT part!
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
ch-kr
left a comment
There was a problem hiding this comment.
thank you for adding additional flexibility! I've added some suggestions to streamline the code (and hopefully didn't break the formatting too much)
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
|
a little late, but I also forgot that we've been requesting that each person who adds or changes code in this repo add tests (https://github.com/broadinstitute/gnomad_methods/tree/main/tests). could you add tests for this new function? |
ch-kr
left a comment
There was a problem hiding this comment.
thanks for adding tests! a few more thoughts
| @pytest.fixture | ||
| def sample_interval_expr(self): | ||
| """Fixture to create a sample interval expression.""" | ||
| return hl.literal( |
There was a problem hiding this comment.
do you need this literal? hl.interval returns an IntervalExpression
There was a problem hiding this comment.
Ah this is very silly on my end. I was using hl.Interval instead of hl.interval before so was wondering why the interval function wasn't working (which is reason for adding hl.literal). Will change this
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
| "Input is a list of IntervalExpressions, so function will return an" | ||
| " ArrayExpression of loci within all input intervals." |
There was a problem hiding this comment.
should these be consistent types? the current behavior means the user passes in a list but receives an ArrayExpression
There was a problem hiding this comment.
I feel an ArrayExpression as input makes more sense in that case. That will be the most common input type in hail isn't it?
| return interval_to_pos_range(intervals_expr).map( | ||
| lambda pos: hl.locus( | ||
| intervals_expr.start.contig, | ||
| pos, | ||
| reference_genome=intervals_expr.start.dtype.reference_genome, | ||
| ) | ||
| ) |
There was a problem hiding this comment.
this can be replaced with a call to _make_loci_array
| ) | ||
|
|
||
| loci_arrays = [_make_loci_array(i) for i in intervals] | ||
| if not flatten: |
There was a problem hiding this comment.
for consistency, this probably should be reordered to be if flatten:
| intervals = intervals.key_by( | ||
| locus=hl.locus( | ||
| intervals[interval_field].start.contig, | ||
| intervals._pos, | ||
| reference_genome=get_reference_genome(intervals[interval_field]), | ||
| ) | ||
| ) |
There was a problem hiding this comment.
this can also be replaced by a call to _make_loci_array
| " intervals overlap." | ||
| ) | ||
| else: | ||
| intervals = intervals.distinct() |
There was a problem hiding this comment.
should we add a warning here that distinct() arbitrarily deduplicates, which means that if the input table was annotated with something like gene or transcript ID, this will pick one at random for duplicated loci?
| if not flatten: | ||
| return hl.array(loci_arrays) | ||
| result = hl.flatten(hl.array(loci_arrays)) | ||
| if deduplicate: | ||
| result = hl.array(hl.set(result)) | ||
| return result |
There was a problem hiding this comment.
should there be a check to make sure the positions are sorted?
| result = hl.array(hl.set(result)) | ||
| return result | ||
|
|
||
| intervals_expr = ( |
There was a problem hiding this comment.
it looks like this expression only gets used once (it doesn't get used if this is a Table below); is this code necessary?
There was a problem hiding this comment.
It also gets used in a table to annotate _pos, but it gets wrapped into interval_to_pos_range either way so maybe I should make that a single call upstream
There was a problem hiding this comment.
I think what I meant when I wrote this comment was that you could directly call _make_loci_array on intervals[interval_field], since the only case that is handled below is if isinstance(intervals, hl.Table). however, if you reorder the code as suggested above, you shouldn't need this if/else at all
Co-authored-by: Katherine Chao <kchao@broadinstitute.org>
ch-kr
left a comment
There was a problem hiding this comment.
a few more comments after revisiting
| .. warning:: | ||
| - Overlapping intervals will produce duplicate loci. Use ``deduplicate=True`` |
There was a problem hiding this comment.
| .. warning:: | |
| - Overlapping intervals will produce duplicate loci. Use ``deduplicate=True`` | |
| .. warning:: | |
| - Overlapping intervals will produce duplicate loci. Use ``deduplicate=True`` |
nit
| chromosomes), as it will create extremely large arrays, which may cause | ||
| performance issues. | ||
|
|
||
| Note that intervals that cross chromosomes are currently not supported. |
There was a problem hiding this comment.
maybe we should raise an error if a user passes an interval that crosses chromosomes?
| " IntervalExpressions!" | ||
| ) | ||
|
|
||
| if isinstance(intervals, hl.Table) and ( |
There was a problem hiding this comment.
seeing this function again with fresh eyes, you could reorder how you tackle the different inputs (Table, IntervalExpression, ArrayExpression). the current function structure checks whether the input is a Table 3 separate times, and reordering would remove the extra checks.
by reordering, I mean this function could handle the case of an input IntervalExpression + return, then the case of an input ArrayExpression + return, and finally the case of an input Table
| result = hl.array(hl.set(result)) | ||
| return result | ||
|
|
||
| intervals_expr = ( |
There was a problem hiding this comment.
I think what I meant when I wrote this comment was that you could directly call _make_loci_array on intervals[interval_field], since the only case that is handled below is if isinstance(intervals, hl.Table). however, if you reorder the code as suggested above, you shouldn't need this if/else at all
| interval_field: Optional[str] = None, | ||
| keep_intervals: Optional[bool] = False, | ||
| deduplicate: bool = True, | ||
| flatten: bool = True, |
There was a problem hiding this comment.
there doesn't seem to be a test covering flatten=False in the test suite
| result = hl.flatten(loci_arrays) | ||
| if deduplicate: | ||
| result = hl.array(hl.set(result)) | ||
| result = hl.sorted(result) |
There was a problem hiding this comment.
we should add a test checking that results are sorted when flatten is True
|
|
||
| intervals = intervals.drop(*fields_to_drop) | ||
|
|
||
| if deduplicate: |
There was a problem hiding this comment.
should the other logic branch (deduplicate=False, keep_intervals=True) case also be covered in tests?
| interval = hl.interval( | ||
| hl.locus("chr1", 100, "GRCh38"), | ||
| hl.locus("chr1", 105, "GRCh38"), | ||
| includes_start=False, | ||
| includes_end=False, | ||
| ) |
There was a problem hiding this comment.
this looks like the same interval that is built into the function above. could this be extracted into a fixture?
| ht = hl.Table.parallelize( | ||
| [{"interval": interval, "gene": "GENE1"}], | ||
| hl.tstruct(interval=hl.tinterval(hl.tlocus("GRCh38")), gene=hl.tstr), | ||
| ) |
There was a problem hiding this comment.
there are a few hl.Table.parallelizes in the tests as well, could you update the fixture or create other fixtures to reduce redundancy?
| # Verify that other fields are still present. | ||
| assert all(hasattr(row, "gene") for row in result) | ||
|
|
||
| def test_explode_interval_expression( |
There was a problem hiding this comment.
the multiple test_explode_interval_expression* tests that only change testing whether the start/end are included/excluded can be collapsed into a single@pytest.mark.parametrized function that uses a table with a few rows that vary the start/end inclusion. this is also true of the test_explode_table_*excludes* tests
Feature Addition:
Added a function in
intervals.pywhich explodes hail intervals into per base loci, keeping same reference genome as original interval.Optionally it can return the original interval that the locus belongs to or drop it from the exploded table.