[sfscode-users] Ne vs Na

[Paula T. <pa...@cs...>]

Hi,

I am using SFS_CODE for simulating data that undergoes a complex demographic scenario. In the end I want to compare the parameters I used in the simulation with some estimates. In my estimates, I have everything scaled by Ne. It seems to me though that the parameters in SFS_CODE are scaled by Na (the ancestral population size, the one used in the burn-in period and specified with -N).

The documentation seems a bit ambiguous about the scaling, but, for example, it says that theta is scaled by Ne (page 7). However, looking at the code, it seems that theta and selection parameters are all updated when the population size is changed. I found the following code lines (line 2500, version 20150815):

                  ppars[pop].Nt *= foo;
                  ppars[pop].THETA *= foo;
                  ppars[pop].INSRATE *= foo;
                  ppars[pop].DELRATE *= foo;
                  if(ppars[pop].RHO > DBL_MIN)
                    ppars[pop].RHO *= foo;
                  else if(ppars[pop].fGC > DBL_MIN)
                    ppars[pop].fGC *= foo;
                  for(r=0; r<gpars.R; r++){
                    ppars[pop].GAMMA[r] *= foo;
                    ppars[pop].lambdaN[r] /= foo;
                    ppars[pop].lambdaP[r] /= foo;
                    ppars[pop].normMean[r] *= foo;
                    ppars[pop].normVar[r] *= (foo*foo);

which, seems to me, means that all parameters (at the input) are in fact scaled by Na and not Ne.

Is this correct or am I missing something?

I think it would be very useful if this was explained in the documentation more clearly. It would also be nice to be able to obtain an Ne (using convertSFS_CODE, for example) from a run with a given demography.

Thank you.

Best regards,
Paula Tataru
--------------------------------------------
Bioinformatics Research Centre,
Aarhus University,
Denmark