dadi-cli: fit DFE with neutral+gamma distribution

[Chenlu Di]

Hi Ryan, dadi developers and users,

I tried to modify dadi-cli to fit DFE with a neugamma distribution that a p0 proportion being neutral and the rest fitting a gamma distribution. I added a neugamma distribution as following and modified all the integration functions according to the PDF (example as below). I listed all the changes I made in the attached file. 

def neugamma(xx, params):

"""Return Neutral + Gamma distribution"""

"""params = (shape, scale, pneu)

pneu is the proportion of neutral mutations"""

alpha, beta, pneu = params

xx = np.atleast_1d(xx)

out = (1-pneu)*gamma(xx, (alpha, beta))

# Assume gamma < 1e-4 is essentially neutral

out[np.logical_and(0 <= xx, xx < 1e-4)] += pneu/1e-4

# Reduce xx back to scalar if it's possible

return np.squeeze(out)

# For example: in InferDFE.py, adding

if (sele_dist == "neugamma"):

print ("Using integrate_point_pos for neugamma distribution.")

func = cache1d.integrate_point_pos

else:

func = cache1d.integrate

It didn't print the added information line "Using integrate_point_pos for neugamma distribution.". Do you know what could be wrong with these modifications? 

Thank you very much!!

Chenlu

[Chenlu Di]

Hi Ryan,

I made additional changes to the "integrate_point_pos" function to fit an input params without gamma by replacing the first few lines to 

pdf_params, ppos = params[:-1*Npos], params[-1]

gammapos = 0

ppos_l = [ppos]

gammapos_l = [gammapos]

Best,

Chenlu

[Chenlu Di]

Hi Ryan,

I read through the function of "integrate" and found this should be good for neutral+gamma distribution. The point-mass of neutral mutations can be added in adjusting the density function of the selection coefficient. The issue has been solved. 

Thank you!

Best,

Chenlu

[Ryan Gutenkunst]

Hello Chenlu,

I’m sorry I haven’t been more responsive about this. I’m not sure what you mean by "adjusting the density function”. As I mentioned in the first email you sent about this, that’s likely to be problematic. Attached is a example script that uses the must more robust option that adds a point mass at zero explicitly.

Best,
Ryan