help with posteriors

[Kritika Garg]

Hi All,
         I am using popabc to differentiate between migration and no migration scenario for my data set. I am using 9 microsatellite markers for this. For now I am using simple 2 pop models. I am doing the no migration model first. My prior files are as follows:

1000000 2 2 9

1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

m m m m m m m m m

0

1 10 100000
1 10 100000

1 10 100000

1 1000 1e+06

0
0


2 0.0001 0 0.0005 0
0

0
0

----------------------------------------------------------------------------------------
PopABC - Mark Beaumont & Joao Lopes                                             01/05/09

>no_iterations, generation_time, no_populations, no_loci

>escalar per locus (autosome - 1; X-linked - 0.75; Y-linked or mitDNA - 0.25)

>type of DNA data (s - sequence; m - microssatelites)

>topology:       0 - uniform distribution;
                 3 - uniform distribution (and choose a Model marker).

>ne1 params:     1 - uniform distribtuion;
                 2 - generalized gamma distribution.
>ne2 params

>neanc1 params

>t1 params:      1 - uniform distribtuion;
                 2 - generalized gamma distribution.

>mig1 params:    0 - zero migration;
                 1 - uniform distribtuion;
                 2 - generalized gamma distribution;
                 3 - uniform distribution (on number of migrations);
                 4 - generalized gamma distribution (on number of migrations).
                 [for 3 and 4 real mig rate is calculated as nmig/Ne]
>mig2 params

>mutM params:    0 - zero mutation;
                 1 - lognormal distribution: (mean of mean(log10); stdev of mean(log10);
                 mean of Sdev(log10); stdev of stdev(log10). Stdev truncated at 0.
                 2 - normal distribution: (mean of mean; stdev of mean; mean of Sdev;
                 stdev of stdev. Stdev truncated at 0.
>mutS params

>recM params:    0 - zero mutation;
                 1 - lognormal distribution: (mean of mean(log10); stdev of mean(log10);
                 mean of Sdev(log10); stdev of stdev(log10). Stdev truncated at 0.
                 2 - normal distribution: (mean of mean; stdev of mean; mean of Sdev;
                 stdev of stdev. Stdev truncated at 0.
>recS params
----------------------------------------------------------------------------------------
Tree topology:

      ||       PopA1
      ||         |
      ||         |
    t1||     ---------
      ||     |       |
      ||     |       |
      \/   Pop1    Pop2

I have a few questions regarding the out put.
1) The mutation rate I have given is taken as a constant in my reject file. I am not sure why there is no variation in microsatellite mutation rate.

Output reject file:

top	avMutM	sdMutM	avRecM	sdRecM	t1
1	0.0001	0.0005	0	0	75389.7
1	0.0001	0.0005	0	0	84598.5
1	0.0001	0.0005	0	0	68973.1
1	0.0001	0.0005	0	0	40349
1	0.0001	0.0005	0	0	20459.7
1	0.0001	0.0005	0	0	21897.8
1	0.0001	0.0005	0	0	69626.9
1	0.0001	0.0005	0	0	9783.84
1	0.0001	0.0005	0	0	11719.4
1	0.0001	0.0005	0	0	30760.3
1	0.0001	0.0005	0	0	135635
1	0.0001	0.0005	0	0	2323.35
1	0.0001	0.0005	0	0	84041.4
1	0.0001	0.0005	0	0	15331.2
1	0.0001	0.0005	0	0	33732.9
1	0.0001	0.0005	0	0	13914.2
1	0.0001	0.0005	0	0	20896
1	0.0001	0.0005	0	0	69444.3
1	0.0001	0.0005	0	0	8421.47
1	0.0001	0.0005	0	0	74364.3
1	0.0001	0.0005	0	0	163115
1	0.0001	0.0005	0	0	210900
1	0.0001	0.0005	0	0	44407.9
1	0.0001	0.0005	0	0	2485.12
1	0.0001	0.0005	0	0	44820.8
1	0.0001	0.0005	0	0	12439.9
1	0.0001	0.0005	0	0	12941.8

2) To evaluate if my priors are exhaustive I generated the prior-posterior plots for my data set. My problem is I don't get good posteriors for NeA and tev. Can you please suggest how I should vary my priors to get better posteriors. I have attached the graphs for the above mentioned parameters.

Thank you.

Regards,
Kritika

[Joao Sollari Lopes]

Hi Kritika,

Regarding question 1), what is saved in the .dat file are the values taken by the hyper priors of the mutation rate. This is fine. If you want the results for each locus, you have to run popabc as:

>simulate toy.prs toy.ssz toy.sst outfile 0 1 0

or run popabc using the menu interface and choose to create file with mutation rate when asked.
This will create a .mut file with a column for each locus with the mutation rate used for each simulation

Regarding question 2), are you using the regression-step? if so sometimes it helps to logit transform the data before performing this step. In the R script choose transf="logit" when using function makepd4(). If you're not using the regression-step, use it.
In any case, it seems that if you're not considering migration the estimated of your splitting time are of almost 0, this is characteristic of large migration rates (or that the two populations are actually just a single one). As for NeA, it seems that the posterior is squashed into the upper boundary of the prior. Consider using a bigger upper boundary.

Hope this helps,
Joao