ancient pseudo diploid individual and numeric issue

[jqjqyangjq]

Thanks for developing the method.
I am likely facing numeric issues with adding pseudo-diploid data from an ancient individual to the dataset, so in the sfs.gz file, the allele count is always [0,2] or [2,0] or [0,0] (missing).
The warnings I received:
autograd/tracer.py:48: RuntimeWarning: divide by zero encountered in log
return f_raw(*args, **kwargs)
/momi/data/sfs.py:408: RuntimeWarning: invalid value encountered in multiply
ret = -lambd + lambd * np.log(lambd) - scipy.special.gammaln(lambd + 1)

In an out-of-Africa model, this always leads to an unexpectedly large population size of diploid populations that are related to the pseudo-diploid ancient individual.
Alternatively, if I followed the publication, setting the population size of the ancient individual to that of an ancestral lineage, the population size can also be extremely inflated (e.g., 10000000000).

Modifying [0,2] [2,0] to [0,1] [1,0] does not seem to help. specifying model.set_data(sfs, non_ascertained_pops=["ancient_individual"]) does not solve this either.

To me, this seems to be an issue of estimating the population size of the pseudo-diploid individual.

Could you maybe suggest how one should deal with the pseudo-haploid / pseudo-diploid data in this case? Thanks

[jackkamm]

All 3 of these things are needed:

• Use [0,1][1,0] instead of [0,2][2,0]
• Use non_ascertained_pops
• Fix pop size of the ancient individual to be equal to its ancestor, or to any constant arbitrary value (doesn't matter what it is, with just a single allele it should have no affect on the likelihood)

Did you try all 3 of these things together? From your description it sounds like you may have tried them individually/separately. But all 3 are necessary to get a sensible result.

[jqjqyangjq]

Thanks for your suggestions. The version I am using is 2.1.20. I did a quick test similar to what was described in the publication.

1. I am including diploid genomes of high-coverage Western Hunter-gatherer(Loschbour), and haploid genomes, including Malta1 (MA1).
2. I first restricted genotypes to all bi-allelic transversion SNP sites from the TGP project (NYC high-coverage). In the end, I only retained segregating positions between two YRI, two CHB, and high-coverage diploid genomes.
3. Now the allele counts for haploid genomes are either [0,0],[0,1], or [1,0]. and "non_ascertained_pops=['MA1']" is also provided. (I tried to add "Louschbour" in the non-ascertained pops as well)

Without adding haploid genomes, I did not receive any warnings, and all parameters made sense. However, I started having the same warnings again, after adding MA1 (in a model including YRI, CHB, Loschbour, MA1), no matter what values/parameters I provided for the population size of Malta1. The estimates make sense to me, though, including divergence time between CHB and YRI (~85ky), Western Eurasians and Eastern Eurasians(45ky), Western Eurasians and MA1(40ky), and population sizes.

However, after replacing MA1 with another new haploid genome, the estimates of population sizes looked reasonable, but the estimates of split time looked much younger than it should have been (so close to the age of the sample).

I wonder if the warnings would affect the inference at all, and if so, do you have any idea what might be causing the numeric issues?