Correct way to simulate admixture model in fastsimcoal

[高橋大樹]

Dear all

I have a probably simple problem with fsc25221.

I'm trying admixture model of three species, the hybrid origin of one species from other two species, using microsatellite data.

In this scenario, species 1 were originate at time T_2 and derived from (H_PRO) of species 0 genome and (1 - H_PRO) of species 2 genome, and at T_1, species 0 and species 2 fissioned to a population. Without these, no migration between populations was existed.

Current population size of species 0, 1 and 2 are denoted as N-0, N-1, N-2, and population size of ancestral lineage at time T_1 was N-2'.

I want to estimate these 8 parameters (T_1, T_2, H_PRO, N-0, N-1, N-2, N-2', ssr- mutation rate)  by using fsc25221.

But computations were too slow, and do　not progressing.

.tpl-file and .est-file were pasted bellow.

-----------------------------------------------------------------tpl-file-------------------------------------------------------------------------------------------------------------------

//Parameters  for the coalescence simulation program :

3 samples to simulate

//Deme sizes (haploid number of genes)

N_0 

N_1 

N_2

//Sample sizes

96

272

160

//Growth rates

0

0

0

//Number of migration matrices : If 0 : No migration between demes

0

//Historical event: time, source, sink, migrants, new deme size, new growth rate, new migration matrix

3 historical events

T_2 1 0 H_PRO 1 0 0

T_2 1 2 H_REPRO 1 0 0

T_1 0 2 1 N_REL 0 0

//Number of independent chromosome 

1 0

//Number of contiguous linkage blocks

1

//Per Block: Data type, No. of loci, Recombination rate to the right-side locus, plus optional parameters

MICROSAT 18 0.000 STR_MUTATION  0 0

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

--------------------------------------------------------------est-file----------------------------------------------------------------------------------------------------------------------

[PARAMETERS]

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

//all N are in number of diploid individuals

1 N_0 logunif 1000 10000

1 N_1 logunif 3000 20000

1 N_2 logunif 2000 20000

1 N_2' logunif 5000 35000

0 H_PRO unif 0 1

1 T_1 unif 2000 20000

1 T_2 unif 1000 10000

0 STR_MUTATION unif 0.00001 0.0001

[RULES]

T_1 > T_2

N_2' > N_2

[COMPLEX PARAMETERS]

0 H_REPRO = 1 - H_PRO

0 N_REL = N_2' / N_2

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

My question is: what is the correct way to progress the computations under this model?

Thanks a lot for insights.

Best,

Daiki takahashi

[Laurent Excoffier]

Hi,

there are several issues with your simulations, perhaps not exhaustive...

1) If you want to simulate admixture, then instead of using this event

T_2 1 2 H_REPRO 1 0 0

try this

T_2 1 2 1 1 0 0

as it will transfer all remaining lineages in deme 1 to deme 2. Else some lineages from deme 1 will remain in deme 1

2) N_2 is included in N_2' which is not allowed or can lead to bad input file. Look at the par file generated bx fsc, and the simparam file in the result directory. See manual and several discussion items in this group

3) I strongly advise not to use rules. Use another parameterisation. Also look at several discussions on this issue in the group.

For instance just define N_REL as a normal parameter, but define it to vary between 1 and something larger than 1.

4) Given your large sample sizes, simulations will take time indeed...

Hope it helps

laurent

[高橋大樹]

dear. Laurent

Thank you very much for your advice and the program could run successfully.

Daiki takahashi

[guill...@gmail.com]

Dear Laurent,

I'm trying to set a model very similar to this and have a question about:

Wouldn't these 2 events combined:

T_2 1 0 H_PRO 1 0 0

T_2 1 2 1 1 0 0

... force the simulation to set H_PRO to zero? As it reads, deme 2 seems to be the only sink population, receiving all the migrants (reading backwards).

Thanks,

Guillermo Friis

[shiy...@gmail.com]

Hi, Guillermo Friis,

You cannot "reading backwards", because the fsc is coalescent.

[Eve Taylor-Cox]

Hi all, 

I am trying to produce a similar model as in this question. I have followed the suggestions above but I think there is an error as I receive the error:
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 noMig_modelD_bad.par

The model I am trying to create is like so (without migration in the first instance): 

Here are my .tpl and .est files

//Parameters for the coalescence simulation program : fsimcoal2.exe
3 samples to simulate :
//Population effective sizes (number of genes)
NENGLAND
NNIRELAND
NSCOTLAND
//Samples sizes and samples age
138
12
76
//Growth rates  : negative growth implies population expansion
0
0
0
//Number of migration matrices : 0 implies no migration between demes
0
//historical event: time, source, sink, migrants, new deme size, new growth rate, migration matrix index
3 historical event
TADM 1 0 H_PRO 1 0 0
TADM 1 2 1 1 0 0
TDIV2 2 0 1 RESIZE 0 0
//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   MUTRATE

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

[PARAMETERS]
//#isInt? #name #dist.#min #max

//all N are in number of haploid individuals
1       ANCSIZE unif    1000    1000000 output
1       NENGLAND        unif    1000    1000000 output
1       NNIRELAND       unif    1000    1000000 output
1       NSCOTLAND       unif    1000    1000000 output
1       TADM    unif    10      1000    output
0       H_PRO   unif    0       1       output
1       TPLUSDIV        unif    10      10000   output
0       MUTRATE unif    2.9e-9  5.5e-9  output

[RULES]

[COMPLEX PARAMETERS]
0       RESIZE = ANCSIZE/NENGLAND output
0       TDIV1 = TADM+TPLUSDIV output

I imagine it is a simple error - in which case please accept my apologies, I am very new to fastSimCoal2.

Thanks in advace for your help,

Eve

[Eve Taylor-Cox]

Oh! Im so sorry - just seen that I have called it TDIV1 not TDIV2 in the .est file.

Why is it you always spot the error when you post a question?!

Sorry - thanks