MRCA not found

[qspreng...@gmail.com]

I am trying to model the divergence of three population. I use the following command line :

./fsc25221 -t SDKFZItest2.tpl -n 100000 -N 100000 -d -e SDKFZItest2.est -M 0.001 -l 10 -L 40 -q --multiSFS -c3 -B3 

I keep getting the following error: 

Estimating model parameters using 3 batches and 3 threads

Simulations did not converge. MRCA not found! (2 lineages remaining)

Num. steps in coalescent process: 1000001

Bad parameters leading to non-convergent coalescent are written in file SDKFZItest2_bad.par

^[[ASegmentation fault: 11

My PI and I have been looking over the .est and .tpl files and can not find any obvious errors. I have included both the .est and .tpl files along with a graphical representation of the model (parameters labeled).  Does anyone see a glaring error?  

Also, is there a hard limit on coalescence time when estimating large populations? 

Any and all help is much appreciated! 

Best,

Quentin 

[Laurent Excoffier]

Most likely what happens is that some lineages get stuck in some populations after migration has stopped, which is due to bad specification of the model.

As mentioned in the manual, avoid making rules including te same parameter several times, as this is not handled well by fastsimcoal2.

Instead of rules, try to use dummy parameters. For instance if you want that T2>T1, define a parameter DELTA_T and then compute T2=T1+DELTA_T

In this way you'll have a much better control over your model

Also try a simpler model that works and complexify it progressively. Your current model is very complex and it seems very difficult to debug 

Hope it helps

laurent

[Renata Pirani]

Hey everyone....

I have the same mistake as below and I already checked all the suggestions before. 

Estimating model parameters using 12 batches and 1 thread

Simulations did not converge. MRCA not found! (2 lineages remaining)

Num. steps in coalescent process: 1000001

Bad parameters leading to non-convergent coalescent are written in file modelo1_bad.par

Simulations did not converge. MRCA not found! (2 lineages remaining)

Num. steps in coalescent process: 1000001

My command line is: ./fsc26 -t modelo1.tpl -n 100000 -N 250000 -m --removeZeroSFS --multiSFS -e modelo1.est -M -l 10 -E 1 -L 40 -q

And here is my .tpl and .est

//Number of population samples (demes)

3 samples to simulate : 

//Population effective sizes (number of genes) 

84461

NPOP1

NPOP2

//Samples sizes and samples age 

28

22

8 

//Growth rates: negative growth implies population expansion 

0 

0 

0

//Number of migration matrices : 0 implies no migration between demes 

2 

//Migration matrix 0 

0 0 0

0 0 0

0 MIG12 0

//Migration matrix 1 

0 0 0

0 0 0

0 0 0

//historical event: time, source, sink, migrants, new deme size, growth rate, migr mat index 

2 historical event 

TDIV1 0 2 1 RESIZE1 0 1

TDIV2 1 0 1 RESIZE2 0 1 

//Number of independent loci [chromosome] 

1 0 

//Per chromosome: Number of contiguous linkage Block: a block is a set of contiguous loci 

1 

//per Block:data type, number of loci, per gen recomb and mut rates 

FREQ 1 0 3.46e-8 OUTEXP

// Priors and rules file

// *********************

[PARAMETERS]

//#isInt? #name   #dist.#min  #max 

//all Ns are in number of haploid individuals   

1  NPOP1       unif     100  500000   output   

1  NPOP2       unif     100  500000   output  

1  N_ANCAll    unif     100  500000   output

1  N_ANC01     unif     100  500000   output

1  TDIV1     unif     100  500000   output     

1  TDIV2       unif     1000  3000000  output   bounded

0  NMig12  logunif     1e-8   1e2        hide    

[RULES]

[COMPLEX PARAMETERS]

0 RESIZE1 = N_ANC01/NPOP2  hide

0 RESIZE2 = N_ANCAll/N_ANC01  hide

0 MIG12 = NMig12/NPOP2 output 

I really appreciate any suggestions. 

Thank you 

[Laurent Excoffier]

Hi,

if TDIV1 < TDIV2 then some lineages are stuck in populations 1 and 2 without being able to coalesce...

L

[Renata Pirani]

Hey Laurent, thank you!!

Yes it is. So how can I fix it?  

[Renata Pirani]

Thank you!!!

It is working well!!!

[DAVID PALACIO]

how did you fix it?

I hope you could help me

[zhao wei]

Hi, 

How do you fix this problem, can you share it with us?

Best 

[huijie liu]

Hi,

Use fsc27 ( paramInRange ) to constraints, this problem can be easily solved. 

Huijie

[Noelia Guzmán]

Hi everyone, I have been trying all the suggestion but I have the same problem as follow:

Simulations did not converge. MRCA not found! (3 lineages remaining)

Num. steps in coalescent process : 1000001

Num. lineages remaining          : 3
Current time                     : 3.40265e+50
Last historical event processed  : 4

Num. steps in coalescent process: 1000001

Bad parameters leading to non-convergent coalescent are written in file model_final_A_bad.par

My .tpl file and .est file are :

//Number of population samples
5

//Population effective sizes (number of genes)

NPOP0
NPOP1
NPOP2
NPOP3
508908
//Samples sizes
10
10
10
10
10

//Growth rates : negative growth implies population expansion
0
0
0
0
0
//Number of migration matrices : 0 implies no migration between demes
2
//Migration matrix 0
0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0
//Migration matrix 1
0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

//historical event: time, source, sink, migrants, new deme size, new growth rate, migration matrix index
4 historical event
TDIV1 1 4 1 RESIZE1 0 1
TDIV2 0 3 1 RESIZE2 0 1
TDIV3 1 2 1 RESIZE3 0 1
TDIV4 1 3 1 RESIZE4 0 1

//Number of independent loci [chromosome]
1 0
//Per chromosome: Number of contiguous linkage Block: a block is a set of contiguous loci
1

//per Block:data type, number of loci, per generation recombination and mutation rates and optional parameters
FREQ 1 0 2.8e-9 OUTEXP

--------------------------------------------------------------------------------------------------------------------

-------------------------------------------------------------------------------------------------------------------

/ Priors and rules file
// *********************

[PARAMETERS]
//#isInt? #name   #dist.#min  #max
//all Ns are in number of haploid individuals

1  ANCSIZE unif     100  5000000   output
1  NPOP7 unif     100  5000000   output
1  NPOP6 unif     100  5000000   output
1  NPOP5 unif     100  5000000   output
1  NPOP3 unif     100  5000000   output
1  NPOP2 unif     100  5000000   output
1  NPOP1 unif     100  5000000   output
1  NPOP0 unif     100  5000000   output
1  TDIV1 unif     100  10000000   output
1  TDIV2 unif     100  10000000   output
1  TDIV3 unif     TDIV1  10000000  output
1  TDIV4 unif     TDIV3 10000000  output

[COMPLEX PARAMETERS]
0  RESIZE1 = NPOP5/508908 hide
0  RESIZE2 = NPOP6/NPOP3 hide
0  RESIZE3 = NPOP7/NPOP2 hide
0  RESIZE4 = ANCSIZE/NPOP6 hide

My command line is: ./fsc27093 -t model_final_A.tpl -e model_final_A.est --multiSFS -n 100000 -u -L 40 -M -m -0 -c8 -B16 -q

The model is:

I hope you could help me.,

Thanks in advance

[Laurent Excoffier]

Hi,

In the est file, try:
1  TDIV3  unif     TDIV1  10000000  output paramInRange
1  TDIV4  unif     TDIV2 10000000  output  paramInRange

Note thatTDiv2 could still be larger than TDiv4, so that you need also to modify your tpl file to be sure to full coalescence, e.g. as :

TDIV1 4 1 1 RESIZE1 0 1

TDIV2 0 3 1 RESIZE2 0 1

TDIV3 2 1 1 RESIZE3 0 1
TDIV4 3 1 1 RESIZE4 0 1  

Here Pop2 can diverge at any time from Pop1, but it shoudl always leads to a convergence of the coalescent process

Hope it helps

laurent

[Noelia Guzmán]

Thanks Laurent for your help! Its working now!

[Anita Wray]

Hello all, 

I am also coming up with this error and am unsure of what could be causing it. 

My .tpl file and .est file are :

//Parameters for the coalescence simulation program: fastsimcoal.exe
3 samples to simulate :

//Population effective sizes (number of genes)

NMB$
NMQ$
NMC$
//samples sizes
46
44
46

//Growth rates: negative growth implies population expansion
0
0
0

//Number of migration matrices : 0 implies no migration between demes

3

//Migration matrix 0
0 0 0
0 0 0
0 0 0

//Migration matrix 1
0 0 0

MIG1$ 0 MIG1$
0 0 0
//Migration matrix 2

0 0 0
0 0 0
0 0 0
//historical event: time, source, sink, migrants, new deme size, new growth rate, migration matrix index

3 historical event
TMIGSTART$ 0 0 0 1 0 1
TD1$ 2 1 1 1 0 2
TD2$ 1 0 1 1 0 2

//Number of independent loci [chromosome]
1 0
//Per chromosome: Number of contiguous linkage Block: a block is a set of contiguous loci
1
//per Block:data type, number of loci, per generation recombination and mutation rates and optional parameters

FREQ 1 0 1e-9 OUTEXP

// Priors and rules file

// *********************

[PARAMETERS]
//#isInt? #name   #dist.#min  #max
//all Ns are in number of haploid individuals

1 NMC$ unif 200 400000 output
1 NMB$ unif 200 400000 output
1 NMQ$ unif 200 400000 output
1 TD1$ unif 140000 160000 output
1 TD2$ unif 240000 260000 output
1 TMIGSTART$ unif 100 140000 output
0 MIG1$ unif 0.001 0.01 output

[COMPLEX PARAMETERS]

My command line is: ./fsc28-mac -t ${PREFIX}.tpl -e ${PREFIX}.est -m -0 10 -n 10000 -L 75 -s 0 -M

I have a similar model to Renata, except populations 0 and 2 are swapped (pop0 diverges first in evolutionary time, last in fastsimcoal time). I have tried all of the above suggestions and am still stuck (including switching source/sink for TD2$, adding population resizes for the ancestral pop, keeping TD1 and TD2 non-overlapping, etc). Thank you in advance!

Error below:

Simulations did not converge. MRCA not found! (2 lineages remaining)

Num. steps in coalescent process : 1000001

Num. lineages remaining          : 2

Current time                     : 3.40233e+50

Last historical event processed  : 3

Num. steps in coalescent process: 1000001

Bad parameters leading to non-convergent coalescent are written in file unidirectional_admixture_bad.par

[Laurent Excoffier]

Hi,

the problem is that TMIGSTART$  can be larger than  TD1$ and in that case some lineages may be stuck in different demes

Could be safer to write

1 TD1$ unif 140000 160000 output bounded
1 TD2$ unif 240000 260000 output  bounded
1 TMIGSTART$ unif 100 140000 output bounded

best

laurent

[Anita Wray]

Hi Laurent, 

Thank you for your speedy response! To clarify, why would TMIGSTART$ ever be larger than TD1$ if TMIGSTART$ max is 140000 and TD1$ min is 140000? Do the values exceed the min and max provided?

Best, 

Anita

[Laurent Excoffier]

Hi,

if TMIGSTART$ is not explicitly bounded it can exceed 140000  and can thus exceed TD1$

Alternatively you could write:

1 TD2$ unif 240000    260000 output  bounded

1 TMIGSTART$ unif   100 140000 output

1 TD1$    unif  TMIGSTART$   160000 output paramInRange

best

laurent

[Noelia Guzmán]

Hi Dear Laurent, I have been trying new simulations with another parameters and I couldnt found the error. Again the MRCA not found is the error.

This is the .tpl file

//Number of population samples
5
//Population effective sizes (number of genes)
NPOP0
NPOP1
NPOP2
NPOP3
508908
//Samples sizes
10
10
10
10
10
//Growth rates : negative growth implies population expansion
0
0
0
0
0
//Number of migration matrices : 0 implies no migration between demes

3

//Migration matrix 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
//Migration matrix 1

0 0 MIG02 0 0

0 0 0 0 0

MIG20 0 0 0 0

0 0 0 0 0
0 0 0 0 0

//Migration matrix 2

0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
//historical event: time, source, sink, migrants, new deme size, new growth rate, migration matrix index

6 historical event
TINTROG2 0 0 1 1 0 1
TINTROG1 0 0 1 1 0 2
TDIV1 4 1 1 RESIZE1 0 2
TDIV2 0 3 1 RESIZE2 0 2
TDIV3 2 1 1 RESIZE3 0 2
TDIV4 3 1 1 RESIZE4 0 2

//Number of independent loci [chromosome]
1 0
//Per chromosome: Number of contiguous linkage Block: a block is a set of contiguous loci
1
//per Block:data type, number of loci, per generation recombination and mutation rates and optional parameters
FREQ 1 0 2.8e-9 OUTEXP

 and the .est file

// Priors and rules file

// *********************

[PARAMETERS]
//#isInt? #name   #dist.#min  #max
//all Ns are in number of haploid individuals
1  ANCSIZE unif     100  5000000   output
1  NPOP7 unif     100  5000000   output
1  NPOP6 unif     100  5000000   output
1  NPOP5 unif     100  5000000   output
1  NPOP3 unif     100  5000000   output
1  NPOP2 unif     100  5000000   output
1  NPOP1 unif     100  5000000   output
1  NPOP0 unif     100  5000000   output

1  TINTROG2 unif     100  10000000   output bounded
1  TINTROG1 unif     TINTROG2  10000000   output paramInRange
1  TDIV1 unif     TINTROG1  10000000   output paramInRange
1  TDIV2 unif     100  10000000   output bounded

1  TDIV3 unif     TDIV1  10000000  output paramInRange
1  TDIV4 unif     TDIV2 10000000  output paramInRange

0  MIG02        logunif  1e-12 1e-1     output
0  MIG20     logunif  1e-12 1e-1 output

[COMPLEX PARAMETERS]
0  RESIZE1 = NPOP5/508908 hide
0  RESIZE2 = NPOP6/NPOP3 hide
0  RESIZE3 = NPOP7/NPOP2 hide

0  RESIZE4 = ANCSIZE/NPOP6 hide

I see when I plot the par file that there is some problems with the time in which same pops coalesce. How can I fix it?

[Noelia Guzmán]

This is the model