Problem simulating selection with the -es and -ej options used in constructing the demography
Closed this issue · 2 comments
Hello,
I'm simulating admixture with selection and am having a problem simulating a selection event starting before the admixture event. I'm going to paste the code I am using below:
Nac=4975360 # current population size of africa
mu=10^-9
locusLength=10000
theta=echo | awk '{print 4*'$Nac'*'$mu'*'$locusLength'}'
rho=echo | awk '{print (4*'$Nac'*'$locusLength'*'$rho_in')}'
totalSampleSize=145
numberOfSimulations=1
sampleSizeAfrica=0
sampleSizeEurope=0
sampleSizeAmerica=145
scaledNeEurope=0.6276
growthRateEurope=-2.674174e-05
scaledNeAmerica=3.2127
growthRateAmerica=-0.006059296
scaledTimeAdmixture=7.263e-05
proportionAdmixture=0.85
scaledTimeSplitAfricaEurope=0.009798
scaledTimeAfricaExpansion=0.1192009
scaledNeAfricaBottleneck=0.0001246
scaledTimeAfricaCrash=0.1192512
scaledNeAfricaAncestral=1.049994
s=995072
s2=1990144
smu=0.01
When age=7.263e-05 or smaller, this works fine:
./msms/bin/msms $totalSampleSize $numberOfSimulations -N $Nac -t $theta -I 3 $sampleSizeAfrica $sampleSizeEurope $sampleSizeAmerica -n 2 $scaledNeEurope -g 2 $growthRateEurope -n 3 $scaledNeAmerica -g 3 $growthRateAmerica -es $scaledTimeAdmixture 3 $proportionAdmixture -ej $scaledTimeAdmixture 3 2 -ej $scaledTimeAdmixture 4 1 -ej $scaledTimeSplitAfricaEurope 2 1 -en $scaledTimeAfricaExpansion 1 $scaledNeAfricaBottleneck -en $scaledTimeAfricaCrash 1 $scaledNeAfricaAncestral -r $recombinationRate $locusLength -SAA $s2 -SaA $s -SI $age 3 0 0 0 -SFC -Smu $smu -Sp 0.5 -Smark -oTrace > $file
When age>7.263e-05 (e.g. age=7.263e-04), this does not work:
./msms/bin/msms $totalSampleSize $numberOfSimulations -N $Nac -t $theta -I 3 $sampleSizeAfrica $sampleSizeEurope $sampleSizeAmerica -n 2 $scaledNeEurope -g 2 $growthRateEurope -n 3 $scaledNeAmerica -g 3 $growthRateAmerica -es $scaledTimeAdmixture 3 $proportionAdmixture -ej $scaledTimeAdmixture 3 2 -ej $scaledTimeAdmixture 4 1 -ej $scaledTimeSplitAfricaEurope 2 1 -en $scaledTimeAfricaExpansion 1 $scaledNeAfricaBottleneck -en $scaledTimeAfricaCrash 1 $scaledNeAfricaAncestral -r $recombinationRate $locusLength -SAA $s2 -SaA $s -SI $age 2 0 0 -SFC -Smu $smu -Sp 0.5 -Smark -oTrace > $file
Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: 2
at at.mabs.model.FrequencyState.getFrequency(FrequencyState.java:56)
at at.mabs.model.selection.SelectionData.setFrequencyToEnd(SelectionData.java:255)
at at.mabs.model.ModelHistroy.simulateSelection(ModelHistroy.java:397)
at at.MSLike.run(MSLike.java:299)
at at.MSLike.runThreads(MSLike.java:424)
at at.MSLike.runMe(MSLike.java:191)
at at.MSLike.main(MSLike.java:541)
at at.MSLike.main(MSLike.java:585)
Note that in the second case, I changes the -SI option to reflect that there should be only 2 populations present. I tried it also with the -SI options with 3 populations present. I get a similar error message:
Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: 3
at at.mabs.model.FrequencyState.getFrequency(FrequencyState.java:56)
at at.mabs.model.selection.SelectionData.setFrequencyToEnd(SelectionData.java:255)
at at.mabs.model.ModelHistroy.simulateSelection(ModelHistroy.java:397)
at at.MSLike.run(MSLike.java:299)
at at.MSLike.runThreads(MSLike.java:424)
at at.MSLike.runMe(MSLike.java:191)
at at.MSLike.main(MSLike.java:541)
at at.MSLike.main(MSLike.java:585)
I've tried this out with simpler demorgaphic models and can't seem to get this to work when I have the -es and ej options used. I can simulate selection in a past epoch if I just use the -ej options, but when I introduce -es and try to create an admixed population, I get the error messages I pasted above.
If you could help me here, I would really appreciate it.
Nandita
So I have finally found the problem, but the final fix is going to be another day or so. However we have a workaround.
First off there is a mistake in the command line. The -SI switch should have the number of demes at that point in time. Hence in this case its 4 when the time is after the -es time. Fixing this however does in fact expose a real bug.
The workaround is to place the 2nd -ej option that is at the same time as the -es option to be just fractionally after the -es option. This won't change the accuracy but will ensure the correct processing order of the switches and should now work fine.
I am closing this due to the fact that we don't intend to directly fix this without the refactor.