search

delt0r / msms

A coalsecent simulator with selection
www.mabs.at/ewing/msms
18 stars 7 forks source link

Problem simulating selection with the -es and -ej options used in constructing the demography #12

Closed ngarud closed 11 years ago

ngarud commented 11 years ago

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

delt0r commented 11 years ago

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.

delt0r commented 11 years ago

I am closing this due to the fact that we don't intend to directly fix this without the refactor.