QC for Gamma_H17 (HomSap)

[petrelharp]

PR for new model: #1099
Original paper: https://doi.org/10.1073/pnas.1619508114

Potential issues that might lead to differences between the model implementations:
@chriscrsmith couldn't find the proportions that @xin-huang used in the paper.

QC'er requests: @xin-huang

[xin-huang]

Should I QC this DFE? The proportions is in the supplementary https://www.pnas.org/doi/full/10.1073/pnas.1619508114#supplementary-materials

[petrelharp]

That'd be great!

[xin-huang]

Hi @petrelharp

how to qc a DFE model? Should I follow the review process for demographic models and put the codes in ./stdpopsim/qc/HomSap.py ?

[grahamgower]

This is a dupe of #1054.

[petrelharp]

I'll close this in favor of #1054, but it looks like we need to actually decide on a QC procedure and implement it - in #1055 we said that it should copy the demographic model QC method, but the first step is to write the code to do that. Do you feel up to starting that, @xin-huang?

[xin-huang]

ok, but what are the codes here for ?

[petrelharp]

ok, but what are the codes here for ?

That's doing a comparison of the output of a simulation with a DFE to what's expected in some other way. Instead, we're writing code so that someone can re-implement the DFEs independently and then the code will automatically compare the original implementation to the QC implementation, and error on discrepancies. How about I'll get this started?

[xin-huang]

ok, but what are the codes here for ?

That's doing a comparison of the output of a simulation with a DFE to what's expected in some other way. Instead, we're writing code so that someone can re-implement the DFEs independently and then the code will automatically compare the original implementation to the QC implementation, and error on discrepancies. How about I'll get this started?

sounds good

[petrelharp]

Well, that was annoying, but I finally figured out where the proportions [0.3, 0.7] come from. It's here:

confirmation:

> 1/(1 + 2.31)
[1] 0.3021148

Maybe the code should be changed to reflect this; i.e., instead of

        proportions=[
            0.3,
            0.7,
        ],  # [0.3 and 0.7 were used in Xin's analysis,
        #  but I couldn't find these values in the Huber paper]

do

# p. 2 in supplement says that the total sequence length of synonymous sites LS
# related to the total sequence length LNS by LNS = 2.31 * LS
proportion_synonymous = 1 / (1 + 2.31)
# <snip>
        proportions=[
            proportion_synonymous,
            1 - proportion_synonymous,
        ],

Hm - on second thought, 1 / (1 + 2.31) = 0.3021148036253776, so 0.3 is just as good and nicer to write, but we should put the explanation in.

[clararehmann]

It looks like this is addressed by commit 8def1f7?