modelling a genetic barrier

[Yorgos]

Dear popABC developers,

I would also like to congratulate you on creating this really useful
software! My question is the following:

I want to model the genetic barrier among, say, 4 populations with the
following topology (in order to visualise the topology properly,
please change to courier font):

|| Neanc7
|| |
t3|| ----------------
|| | |
|| Neanc5 Neanc6
|| | |
|| | |
t2|| --------- ---------
|| | | | |
|| | | | |
|| | | | |
t1|| Neanc1 Neanc2 Neanc3 Neanc4
|| | | | |
|| | | | |
|| | | | |
\/ Ne1 Ne2 Ne3 Ne4

This genetic barrier was supposedly present between t1 and t3 in the
past (which is translated as migration happening only between Neanc1
and Neanc2 on one side and Neanc3 and Neanc4 on the other), but was
lost from t1 onwards (which could be represented by a migration matrix
without any zero non-diagonal elements between any pair of present-day
populations).

My question is: does popABC support this kind of topology (e.g.
connecting Ne1 and Neanc1 without a coalescent event) or is there a
better way to model this?

In case my approach is correct, would I then just have to fill in a
migration matrix like the following one?

Ne1 Ne2 Ne3 Ne4

Ne1 t1 x t1 x t1 x t1 x
t2 x t2 x t2 x t2 x
t3 x t3 x t3 x t3 x

Ne2 t1 x t1 x t1 x t1 x
t2 x t2 x t2 x t2 x
t3 x t3 x t3 x t3 x

Ne3 t1 x t1 x t1 x t1 x
t2 x t2 x t2 x t2 x
t3 x t3 x t3 x t3 x

Ne4 t1 x t1 x t1 x t1 x
t2 x t2 x t2 x t2 x
t3 x t3 x t3 x t3 x

Thank you for you time!

Yorgos

[Joao Sollari Lopes]

Hi Yorgos,

The population model that you need to set is this one (best viewed with courier font)

    ||              Nea5
    ||               |
  t5||        -------------
    ||        |           |
    ||      Nea4          |
    ||        |           |
  t4||     ------         |
    ||     |    |         |
    ||     |    |        Nea3
    ||     |    |         |
  t3||     |    |    ----------
    ||     |    |    |        |
    ||     |    |    |       Nea2
    ||     |    |    |        |
  t2||     |    |    |    ---------
    ||     |    |    |    |       |
    ||     |    |    |    |      Nea1
    ||     |    |    |    |       |
  t1||     |    |    |    |    ------
    ||     |    |    |    |    |    |
    ||     |    |    |    |    |    |
    \/    Ne1  Ne2  Ne3  Ne4  Ne5  Ne6

where you have to set 2 ghost-populations (pop with 0 samples) to set the barrier formation and disappearance events.
This model is coded in the program as

    ||     |             
  t5||     |<---------   
  t4||     |<----    |   
  t3||     |    |    |<----   
  t2||     |    |    |    |<----
  t1||     |    |    |    |    |<----
    ||     |    |    |    |    |    |
    \/     0    1    2    3    4    5

which has the following population-joining events (from present to past)

t1:5->4; t2:4->3; t3:3->2; t4:1->0; t5:2->0

The migration matrix should be something like this:

              pop0 pop1 pop2 pop3 pop4 pop5

pop0  now->t1 0    0.34 0.33 0.33 0    0
      t1->t2  0    1    0    0    0    0
      t2->t3  0    0.34 0.33 0.33 0    0
      t3->t4  0    0.5  0.5  0    0    0
      t4->t5  0    0    1    0    0    0
     
pop1  now->t1 0.34 0    0.33 0.33 0    0
      t1->t2  1    0    0    0    0    0
      t2->t3  0.34 0    0.33 0.33 0    0
      t3->t4  0.5  0    0.5  0    0    0
      t4->t5  0    0    0    0    0    0

pop2  now->t1 0.33 0.33 0    0.34 0    0
      t1->t2  0    0    0    1    0    0
      t2->t3  0.33 0.33 0    0.34 0    0
      t3->t4  0.5  0.5  0    0    0    0
      t4->t5  1    0    0    0    0    0

pop3  now->t1 0.33 0.33 0.34 0    0    0
      t1->t2  0    0    1    0    0    0
      t2->t3  0.33 0.33 0.34 0    0    0
      t3->t4  0    0    0    0    0    0
      t4->t5  0    0    0    0    0    0

pop4  now->t1 0    0    0    0    0    0
      t1->t2  0    0    0    0    0    0
      t2->t3  0    0    0    0    0    0
      t3->t4  0    0    0    0    0    0
      t4->t5  0    0    0    0    0    0
   
pop5  now->t1 0    0    0    0    0    0
      t1->t2  0    0    0    0    0    0
      t2->t3  0    0    0    0    0    0
      t3->t4  0    0    0    0    0    0
      t4->t5  0    0    0    0    0    0


So, in the end you will have the following model
     
    ||            Nea5
    ||             |  
  t5||        ------------
    ||        |          |
    ||       Nea6        |
    ||        |         Nea5
  t4||     ------        |
  t3||     |    |     ------
    ||     |    |     |    |
    ||     |    |     |    |
  t2||     |    |  B  |    |
    ||     |    |  A  |    |
    ||     |    |  R  |    |
    ||     |    |  R  |    |
    ||     |    |  I  |    |
    ||     |    |  E  |    |
  t1||     |    |  R  |    |

    ||     |    |     |    |
    ||     |    |     |    |
    \/    Ne1  Ne2   Ne3  Ne4

where t1 is when the barrier disappears and t2 when the barrier is created. After t1 and before t2 there is free gene flow between all the populations. Unfortunately it is not possible to shut down migration in popabc, so you end up having to consider gene flow previously to the formation of the barrier (the reasons are further discussed in https://groups.google.com/forum/?fromgroups#!topic/popabc/rWe1otYypGw; also, credits for using ghost populations to change migration patterns go all to Shu-ping!)

Note that:
Ne1 is the size of Population 1
Ne2 is the size of Population 2
Ne3 is the size of Population 3
Ne4, Nea1 and Nea2 are the sizes of Population 4 up to t1, t2 and t3, respectively
m1 is the portion of migrants of Population 1
m2 is the portion of migrants of Population 2
m3 is the portion of migrants of Population 3
m4, ma1 and ma2 are the portions of migrants of Population 4 up to t1, t2 and t3, respectively

And also, you should ignore the summary statistics that are related to the ghost populations, since these don't add anymore information.

Finally, here is the prior file you should use:

--
[niter] [gentime] 7 [nloci]

[loci scalars]

[loci types]

2 5 4 4 3 3 2 1 0 2 0

[Ne1 prior]
[Ne2 prior]
[Ne3 prior]
[Ne4 prior]
1 1 1
1 1 1

[Nea1 prior]
[Nea2 prior]
[Nea3 prior]
[Nea4 prior]
[Nea5 prior]

[t1 prior]
[t2 prior]
[t3 prior]
[t4 prior]
[t5 prior]

[m1 prior]
[m2 prior]
[m3 prior]
[m4 prior]
0
0

[ma1 prior]
[ma2 prior]
[ma3 prior]
[ma4 prior]

[STRs mutation prior]
[SNPs mutation prior]

[STRs recombination prior]
[SNPs recombination prior]

1

0    0.34 0.33 0.33 0    0
0    1    0    0    0    0
0    0.34 0.33 0.33 0    0
0    0.5  0.5  0    0    0
0    0    1    0    0    0
     
0.34 0    0.33 0.33 0    0
1    0    0    0    0    0
0.34 0    0.33 0.33 0    0
0.5  0    0.5  0    0    0
0    0    0    0    0    0

0.33 0.33 0    0.34 0    0
0    0    0    1    0    0
0.33 0.33 0    0.34 0    0
0.5  0.5  0    0    0    0
1    0    0    0    0    0

0.33 0.33 0.34 0    0    0
0    0    1    0    0    0
0.33 0.33 0.34 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    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    0    0    0    0    0
0    0    0    0    0    0
0    0    0    0    0    0
0    0    0    0    0    0
--

Ok. Hope this reply doesn't arrive too late.
In any case you should take a look at "Bayesian SSC" (if you haven't already).

Best,
Joao

[Milvago]

Hi Joao, Shu-Ping's approach is very clever and thank you Joao for
describing a model with a temporal barrier to gene-flow. My analysis
requires that after t1 (when barrier fade off) only within sister-
pairs there is exchange of migrants. Can I just change the migration
matrix (using this model) to state this scenario of exchange? or do I
need to modify the tree-model topology?

Second, I want to create a simple competing model to the one mentioned
above. Thus, I am wondering if the model depicted on Fig 9 in popABC
documentation showing tree branching scheme 13 for 4 populations can
be use as an alternative model to compete Joao's model. In fact, I
want that model to represent restricted gene-flow among NON-sister-
pairs and further active exchange within sister-pairs. Can this be
done by only changing the migration matrix? Do I need to consider
using the ghost population approach? If so, could you suggest how to
built this alternative model.

Thank you so much for your effort to support the use of popABC,

Raul

On Apr 5, 8:06 am, Joao Sollari Lopes <j.sollari.lo...@gmail.com>
wrote:

> barrier (the reasons are further discussed inhttps://groups.google.com/forum/?fromgroups#!topic/popabc/rWe1otYypGw;

[Joao Sollari Lopes]

Just some minor amendments to this post.

1) The simplified model is actually:

    ||            Nea5
    ||             | 
  t5||        ------------
    ||        |          |
    ||       Nea4        |
    ||        |         Nea3

  t4||     ------        |
  t3||     |    |     ------
    ||     |    |     |    |
    ||     |    |     |    |
  t2||     |    |  B  |    |
    ||     |    |  A  |    |
    ||     |    |  R  |    |
    ||     |    |  R  |    |
    ||     |    |  I  |    |
    ||     |    |  E  |    |
  t1||     |    |  R  |    |
    ||     |    |     |    |
    ||     |    |     |    |
    \/    Ne1  Ne2   Ne3  Ne4

2) In the detailed model, Ne4 and Nea2 are the sizes of Population 4 up to t2 and t3, respectively. And Nea1 can be discarded since it refers to an "Ancient Ghost population" (!).

3) Finally, in the prior file, [ma1 prior] should be '0'. Again this refers to the "Ancient Ghost population" and we are not considering migration between these populations.

Joao

[Joao Sollari Lopes]

An easier model to consider, that I think best describes what Yorgo wanted (and I miss it in the first place) is the following population model (best viewed with courier font)

    ||              Nea4
    ||               |
  t4||        -------------
    ||        |           |
    ||      Nea3          |
    ||        |           |
  t3||     ------         |

    ||     |    |         |
    ||     |    |       Nea2
    ||     |    |         |
  t2||     |    |    ---------
    ||     |    |    |       |
    ||     |    |    |     Nea1
    ||     |    |    |       |
  t1||     |    |    |    ------
    ||     |    |    |    |    |
    ||     |    |    |    |    |
    ||     |    |    |    |    |
    \/    Ne1  Ne2  Ne3  Ne4  Ne5

where there is only 1 ghost-population (Pop5) to set the barrier disappearance event.

This model is coded in the program as

    ||     |            

  t4||     |<---------  

    ||     |         |
  t3||     |<----    |  
    ||     |    |    |
  t2||     |    |    |<----  
    ||     |    |    |    |
  t1||     |    |    |    |<----
    ||     |    |    |    |    |
    \/     0    1    2    3    4

which has the following population-joining events (from present to past)

t1:4->3; t2:3->2; t3:1->0; t4:2->0

This time we will consider no migration in ancient populations Nea2 and Nea3. So the migration matrix should be something like this:

              pop0 pop1 pop2 pop3 pop4

pop0  now->t1 0    0.34 0.33 0.33 0

      t1->t2  0    1    0    0    0

      t2->t3  0    1    0    0    0

      t3->t4  0    0    0    0    0
    

pop1  now->t1 0.34 0    0.33 0.33 0

      t1->t2  1    0    0    0    0

      t2->t3  1    0    0    0    0

      t3->t4  0    0    0    0    0

pop2  now->t1 0.33 0.33 0    0.34 0

      t1->t2  0    0    0    1    0

      t2->t3  0    0    0    0    0

      t3->t4  0    0    0    0    0

pop3  now->t1 0.33 0.33 0.34 0    0

      t1->t2  0    0    1    0    0

      t2->t3  0    0    0    0    0

      t3->t4  0    0    0    0    0

pop4  now->t1 0    0    0    0    0

      t1->t2  0    0    0    0    0

      t2->t3  0    0    0    0    0

      t3->t4  0    0    0    0    0

      t4->t5  0    0    0    0    0

So, in the end we have the following model:
    
    ||            Nea5
    ||             | 
  t4||        ------------
    ||        |    B     |
    ||       Nea4  A     |
    ||        |    R    Nea3
  t3||     ------  R     |
  t2||     |    |  I  ------

    ||     |    |  E  |    |
  t1||     |    |  R  |    |
    ||     |    |     |    |
    ||     |    |     |    |
    \/    Ne1  Ne2   Ne3  Ne4

where t1 is when the barrier disappears.

Ne1 is the size of Population 1
Ne2 is the size of Population 2
Ne3 is the size of Population 3

Ne4 and Nea1 are the sizes of Population 4 up to t1 and t2, respectively
m1 is migration Pop2->Pop1
m2 is migration Pop1->Pop2
m3 is migration Pop4->Pop3
m4 and ma1 are mirgation Pop3->Pop4 up to t1 and t2, respectively

Finally, here is the prior file you should use:

--

[niter] [gentime] 5 [nloci]

[loci scalars]

[loci types]

2 4 3 3 2 1 0 2 0

[Ne1 prior]
[Ne2 prior]
[Ne3 prior]
[Ne4 prior]
1 1 1

[Nea1 prior]
[Nea2 prior]
[Nea3 prior]
[Nea4 prior]

[t1 prior]
[t2 prior]
[t3 prior]
[t4 prior]

[m1 prior]
[m2 prior]
[m3 prior]
[m4 prior]
0

[ma1 prior]
0
0

[STRs mutation prior]
[SNPs mutation prior]

[STRs recombination prior]
[SNPs recombination prior]

1

0    0.34 0.33 0.33 0

0    1    0    0    0

0    1    0    0    0
0    0    0    0    0

    
0.34 0    0.33 0.33 0

1    0    0    0    0

1    0    0    0    0

0    0    0    0    0

0.33 0.33 0    0.34 0

0    0    0    1    0
0    0    0    0    0

0    0    0    0    0

0.33 0.33 0.34 0    0
0    0    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    0    0
0    0    0    0    0
0    0    0    0    0
--

Joao

[Joao Sollari Lopes]

1.
So, if I got you right you want a population model where the barrier does not fade (best viewed in courirer font):

    ||            Nea4
    ||             | 
  t4||        ------------
    ||        |          |
    ||       Nea3        |
    ||        |         Nea2
  t3||     ------        |
  t2||     |    |     ------
    ||     |    |     |    |
    ||     |    |     |    |
    ||     |    |  #  |    |
  t1||     |    |  #  |    |
    ||     |    |  #  |    |
    \/    Ne1  Ne2 # Ne3  Ne4
   
If that is the case you need to set the following population model:

    ||            Nea4
    ||              |
  t4||        ------------
    ||        |          |
    ||      Nea3         |
    ||        |          |
  t3||     ------        |

    ||     |    |        |
    ||     |    |      Nea2
    ||     |    |        |
  t2||     |    |    ---------
    ||     |    |    |       |
    ||     |    |    |     Nea1
    ||     |    |    |       |
  t1||     |    |    |    ------
    ||     |    |    |    |    |
    ||     |    |    |    |    |
    \/    Ne1  Ne2  Ne3  Ne4  Ne5

where you have to set 1 ghost-populations (Pop5) to set the barrier formation event.

This model is coded in the program as

    ||     |             

  t4||     |<---------

    ||     |         |
  t3||     |<----    |
    ||     |    |    |   
  t2||     |    |    |<----   
    ||     |    |    |    |
  t1||     |    |    |    |<----
    ||     |    |    |    |    |
    \/     0    1    2    3    4

which has the following population-joining events (from present to past)

t1:4->3; t2:3->2; t3:1->0; t4:2->0

 
The migration matrix should be something like this:

              pop0 pop1 pop2 pop3 pop4

pop0  now->t1 0    1    0    0    0
      t1->t2  0    0.34 0.33 0.33 0
      t2->t3  0    0.5  0.5  0    0
      t3->t4  0    0    1    0    0

pop1  now->t1 1    0    0    0    0
      t1->t2  0.34 0    0.33 0.33 0
      t2->t3  0.5  0    0.5  0    0

      t3->t4  0    0    0    0    0

pop2  now->t1 0    0    0    1    0
      t1->t2  0.33 0.33 0    0.34 0
      t2->t3  0.5  0.5  0    0    0
      t3->t4  1    0    0    0    0

pop3  now->t1 0    0    1    0    0
      t1->t2  0.33 0.33 0.34 0    0

      t2->t3  0    0    0    0    0

      t3->t4  0    0    0    0    0

pop4  now->t1 0    0    0    0    0

      t1->t2  0    0    0    0    0

      t2->t3  0    0    0    0    0

      t3->t4  0    0    0    0    0

So,

Ne1 is the size of Population 1
Ne2 is the size of Population 2
Ne3 is the size of Population 3

Ne4 and Nea1 are the sizes of Population 4 up to t1 and t2, respectively

m1 is the portion of migrants of Population 1
m2 is the portion of migrants of Population 2
m3 is the portion of migrants of Population 3

m4 and ma1 are the portions of migrants of Population 4 up to t1 and t2, respectively

Again, you should ignore the summary statistics that are related to the ghost populations.

Here is the prior file you should use:

--
[niter] [gentime] 5 [nloci]

[loci scalars]

[loci types]

2 4 3 3 2 1 0 2 0

[Ne1 prior]
[Ne2 prior]
[Ne3 prior]
[Ne4 prior]
1 1 1

[Nea1 prior]
[Nea2 prior]
[Nea3 prior]
[Nea4 prior]

[t1 prior]
[t2 prior]
[t3 prior]
[t4 prior]

[m1 prior]
[m2 prior]
[m3 prior]
[m4 prior]
0

[ma1 prior]
[ma2 prior]
[ma3 prior]

--

2.
I'm not sure if I got you correctly. If I am not mistaken you want to assume a model with a barrier formation from the origin of the sister-pairs to the present:

    ||           Nea3
    ||            | 
  t3||        -----------
    ||        |    #    |
    ||       Nea2  #    |
    ||        |    #    |
  t2||     ------  #   Nea1
    ||     |    |  #    |
  t1||     |    |  #  ------
    ||     |    |  #  |    |
    ||     |    |  #  |    |
    ||     |    |  #  |    |
    \/    Ne1  Ne2 # Ne3  Ne4

If that is the case, topology 16 is a good starting point (topology 13 would be fine too, the only difference between them is the order of spliting of the sister-pairs. So I'll just use 16 for consistency):

    ||     |             
  t3||     |<---------

    ||     |         |   
  t2||     |<----    |
    ||     |    |    |   
  t1||     |    |    |<----   
    ||     |    |    |    |
    \/     0    1    2    3

Then you need to have a migration matrix such as the one above:
       
              pop0 pop1 pop2 pop3

pop0  now->t1 0    1    0    0

      t1->t2  0    1    0    0  

      t2->t3  0    0    0    0

pop1  now->t1 1    0    0    0
      t1->t2  1    0    0    0  
      t2->t3  0    0    0    0

pop2  now->t1 0    0    0    1

      t1->t2  0    0    0    0  

      t2->t3  0    0    0    0

pop3  now->t1 0    0    1    0  

      t1->t2  0    0    0    0  

      t2->t3  0    0    0    0  

So, the prior file should be something like this:

--
[niter] [gentime] 4 [nloci]

[loci scalars]

[loci types]

1 16

[Ne1 prior]
[Ne2 prior]
[Ne3 prior]
[Ne4 prior]

[Nea1 prior]
[Nea2 prior]
[Nea3 prior]

[t1 prior]
[t2 prior]
[t3 prior]

[m1 prior]
[m2 prior]
[m3 prior]
[m4 prior]

0
0

[STRs mutation prior]
[SNPs mutation prior]

[STRs recombination prior]
[SNPs recombination prior]

1

0    1    0    0
0    1    0    0  
0    0    0    0

1    0    0    0
1    0    0    0  
0    0    0    0

0    0    0    1
0    0    0    0  
0    0    0    0

0    0    1    0  
0    0    0    0  

0    0    0    0  
--

Note that now:
m1 is actually the migration rate Pop2->Pop1
m2 is actually the migration rate Pop1->Pop2
m3 is actually the migration rate Pop4->Pop3
m4 is actually the migration rate Pop3->Pop4

Hope that helps,
Joao