Adapt hl.de_novo function

[KoalaQin]

These are smaller functions modified from Julia's work on combining hl.de_novo and Kaitlin’s code.

[ch-kr]

thanks for adding these functions! I have some questions (many because I'm not that familiar with this work) and suggestions

[ch-kr]

rather than running this in this function and in calculate_de_novo_post_prob, why not switch the order and run get_genomic_context first and pass in the three return values to calculate_de_novo_post_prob as function arguments?

[ch-kr]

I also forgot to ask in my initial review -- could you add tests for the new functions?

[KoalaQin]

I just saw your comment on adding tests, I will work on that now.

[KoalaQin]

Back to you!

[ch-kr]

thanks for adding tests! I have some more questions and suggestions

[ch-kr]

makes sense to me to keep the actual code in gnomad_qc, but can you add this documentation back to the function? I took this from Julia's function:

    .. warning::

        This method assumes that the PL and AD fields are present in the genotype fields
        of the child and parents. If they are missing, this method will not work.
        Many of our larger datasets have the PL and AD fields intentionally removed to
        save storage space. If this is the reason that the PL and AD fields are
        missing, the only way to use this method is to set them to their approximate
        values:

        .. code-block:: python

                PL=hl.or_else(PL, [0, GQ, 2 * GQ])
                AD=hl.or_else(AD, [DP, 0])

this is very clear documentation for our users both as a warning that this function needs AD and PL and how to approximate them if they also need to do so

[ch-kr]

can you also add a version of this documentation back to the docstring?

    The simplified version is the same as Hail's methods when using the
    `ignore_in_sample_allele_frequency` parameter. The main difference is that this
    mode should be used when families larger than a single trio are in the dataset, in
    which an allele might be de novo in a parent and transmitted to a child in the
    dataset. This mode will not consider the allele count (AC) in the dataset, and will
    only consider the Phred-scaled likelihoods (PL) of the child and parents,
    allele balance (AB) of the child and parents, the genotype quality (GQ) of the
    child, the depth (DP) of the child and parents, and the population frequency prior.

the explanation that this function behaves the same as Hail's + ignore_in_sample_allele_frequency and which exact fields the function needs is clear and helpful to have

[KoalaQin]

I rewrote the test simpler, I could parametrize more, but I want you to have a look first.

[KoalaQin]

Now I'm looking at the code and data, isn't hl.sum(AD)=DP? Especially it's split. Can we just use parent DP here?

[ch-kr]

this is a good question. I thought I remembered hearing that sum(AD) isn't exactly equal to DP and found this note:

The AD calculation as performed by HaplotypeCaller may not yield exact results because only reads that statistically favor one allele over the other are counted. Due to this fact, the sum of AD may be different than the individual sample depth, especially when there are many non-informative reads.

https://gatk.broadinstitute.org/hc/en-us/articles/360037592411-DepthPerAlleleBySample

could you test to see if sum(AD) == DP here?

[KoalaQin]

I think the problem is that we don't have AD for hom_ref stored, and we're approximating it back by:

AD=hl.or_else(mt.AD, [mt.DP, 0])

The sum(AD) would be DP for those hom ref parent(s) for sure. Do you know where to find the trace (or documentation) that we intentionally removed AD and PL?

[ch-kr]

oh that makes sense. given this function depends on the presence of AD/PL, let's stick with AD anyway, even if DP would be equivalent for our dataset. other people using this function might still have AD and PL in their datasets.

finding the documentation or trail on this topic might take a while -- my instincts are checking slack (both Broad and ATGU), project meeting notes, and joint calling meeting notes

[ch-kr]

thanks for doing the suggested restructuring! the current structure and formatting looks nice. I have some more minor suggestions and some suggestions for the tests

[ch-kr]

I don't think you need this formula since you've added this comment:

Please refer to these sources for more information on the de novo model.

alternatively, if you feel strongly that these equations should stay in this docstring, particularly given your hard work on the formatting (which looks great!), can you move this

Neither Kaitlin's de novo caller nor Hail's de novo method provide a clear
description on how to calculate for de novo calls for hemizygous genotypes in XY
individuals. These equations are included below:

to after this?

Probability of a de novo mutation given the data for hemizygous calls in XY individuals

[KoalaQin]

I revised this to be more logic.

[ch-kr]

this equation extends past the end of the screen for me, and I don't see an option to scroll horizontally -- do you see this as well when you render the docs?

[ch-kr]

also, in these formulas below, you use x to indicate multiplication, but you don't above, e.g.:

can you add the x above or remove them here? also, if you decide to keep a symbol, is there a way to use a dot instead (as in the Hail docs) to save horizontal space?

[KoalaQin]

I don't see that issue on my big screen but I do see it on my laptop but it will resolve if you zoom in the webpage to 90%. (This got better with \cdot.
I saw for the multiplication, we can use either \times or \cdot (as Hail) or just a small space \, , * is not recommended. Hail also used both \, and \cdot. To keep it consistent, I change it to \cdot everywhere except for the DN prior numbers.

[ch-kr]

the \cdot helped a lot, thanks for making the change! I agree with keeping the \times in the de novo prior. an alternative for prior is to display that formula in the same format as the Hail docs:

(not a necessary edit, just something to consider -- happy with whichever you prefer)

[ch-kr]

rather than having so many examples, could you have just a few with the following conditions:

• two examples where calculate_de_novo_post_prob should return the expected P(dn) values
• one example where calculate_de_novo_post_prob will throw an error because freq_prior_expr is out of the expected range?

[ch-kr]

I chatted with Mike, and he pointed me to Kristen's recently written tests as a good example of what our tests should do: https://github.com/broadinstitute/gnomad_methods/blob/main/tests/assessment/test_validity_checks.py.

with my comment above: we should test for valid uses of the function (making sure valid inputs return expected outputs) and incorrect uses of the function (checking that the function returns errors as expected). you've added these suggestions, but it would also be a good idea to add tests for any edge cases you can think of (if possible).

Mike also mentioned that pytest has a built in caplog functionality to capture logger information, which made me realize these functions don't have any loggers. I'm not sure they're necessary, but I wanted to note it since I thought of it

[ch-kr]

Suggested change

	def test_default_get_de_novo_expr_fail_conditions(self):
	def test_default_get_de_novo_expr(self):

since this is a test function, I don't think it needs to include "fail_conditions" explicitly in the name. also, similar to the comment above, I think this should test a few scenarios:

• a case like the one you've included, where multiple fail conditions apply
• a case where only one fail condition is true
• a case where the variant has False for is_de_novo
• a passing case

[KoalaQin]

I'll leave you to resolve the conversations in case I didn't address them as you suggested.

[KoalaQin]

I changed this back because your suggestion will need more changes below.

[ch-kr]

a few more minor suggestions

[ch-kr]

this explanation is great -- it's really clear!

I have a nitpick here: can you combine these into the same paragraph/text block? these sentences logically belong together, so they don't need the extra newline in between

[ch-kr]

the \cdot helped a lot, thanks for making the change! I agree with keeping the \times in the de novo prior. an alternative for prior is to display that formula in the same format as the Hail docs:

(not a necessary edit, just something to consider -- happy with whichever you prefer)

[ch-kr]

I chatted with Mike, and he pointed me to Kristen's recently written tests as a good example of what our tests should do: https://github.com/broadinstitute/gnomad_methods/blob/main/tests/assessment/test_validity_checks.py.

with my comment above: we should test for valid uses of the function (making sure valid inputs return expected outputs) and incorrect uses of the function (checking that the function returns errors as expected). you've added these suggestions, but it would also be a good idea to add tests for any edge cases you can think of (if possible).

Mike also mentioned that pytest has a built in caplog functionality to capture logger information, which made me realize these functions don't have any loggers. I'm not sure they're necessary, but I wanted to note it since I thought of it

[ch-kr]

ah yes, sorry, that wasn't clear. this suggestion was relevant when you wanted to filter de novos from the input dataset upfront. in that scheme, you would pass an argument to this function (e.g., filter_to_candidate_dnms or something similar), and the default for that argument would be True, which means the function would first filter the input dataset to candidate de novo variants before calculating p-values and confidence levels.

however, in your current function design, where you pass in and return only expressions, I don't think this behavior would make sense; I more wanted to suggest this as an alternative option to your previous structure (filtering the HT to variants with GTs that indicate possibility of being a de novo variant).