GADMA for 5 population

[SachKrVinni]

Hi, Thanks for explaining the protocol for using GADMA.
With these tutorials, I was able to design a few demographic models for three populations with 24 runs, 1800 iterations.
All runs have one final model, now I am confused which one I should select for the concluding results. It will be a great help if you can provide me the best way to select the closest one?
Now, I was just wondering and need your assistance to do demographic modelling for a 5 population of non-model organisms. I have gone through a few examples (inference_for_four_and_five_populations.html#Inference-for-four-and-five-populations), but my concepts are not clear. Just consider me a newbie, and can you explain to me the basic things I need to do first before running SMAC or BO models?
A detailed, step-by-step help will be greatly appreciated.
I have a VCF file with 28 individuals, 5 populations, and the number of individuals are in order of 5,11,7,3,2.

[noscode]

Hi @SachKrVinni,

It is great to hear that you find the materials helpful. I will do my best to help with your questions.

1. Model selection. You mentioned that you have designed a few demographic models for three populations. Does this mean you built your own models, or did you use the structure model option in GADMA?
   Model comparison is an interesting topic. If the models have the same number of parameters (or the same structure), you can compare them using the log-likelihood values directly. However, if they differ in the number of parameters (which is usually the case), the approach depends on whether your data consists of independent SNPs or not.
   For example, if you selected one SNP per RAD-seq region or performed LD pruning before building your SFS, then the SNPs are unlinked. Otherwise, they are considered linked. It is relatively easy to calculate, you just need number of parameters and log-likelihood values.

• If your data is unlinked, you should use the Akaike Information Criterion (AIC) to compare models and select the best one.
• If your data includes linked SNPs, you will need to perform a block-bootstrap, as explained here. After generating the bootstrapped datasets, you can use the Composite Likelihood AIC (CLAIC) to compare the models. It is not strait forward, but I can help you with that.

2. Inference with five populations is more or less the same as a usual inference using GADMA. It differs in two main aspects:
   a) You need to design and specify your own models that describe possible demographic scenarios, you cannot use structure model.
   b) You must use a different optimisation method: Bayesian optimisation.
   The first point is usually the most challenging. The second point does not require a lot, you just need to specify the usage of another optimisation in GADMA params file. It only requires you to use Python 3.8, as the current implementation of Bayesian optimisation only works with that version.
   You can follow the provided example (Inference for four and five populations), but if you want more accurate results, I recommend using "Global maxeval: 500" for more iterations (note that this will take more time).

Best,
Ekaterina

[SachKrVinni]

Thanks Ekaterina, for explaining the details.
Yes, I am a bit confused in the first step only and didn't know the simplest method to generate a model for my data.

1. Model selection: Yes, I used structure model and SNP linked, means I didn't pruned the data. Although I did filter it with vcftools. I didn't even perform bootstrap, will be of great help if you can help in this regard.
   I ran GADMA as described by you with three pop (only added mutation rate), which generated the best and final model in a different run.
   Which is this:
   cat params_file
   Output directory: /home/Gadma_Test/All_out_dir
   Input data: /home/Gadma_Test/chr1_Merge.vcf,
   /home/Gadma_Test/All_pop_file
   Population labels: [Central, NorthEast, NPG]
   Projections: [10, 22, 14]
   Outgroup: False
   Sequence length: 243826021
   Linked SNP's: True
   Directory with bootstrap: Null
   Engine: moments
   Pts: [30, 40, 50]
   Relative parameters: True
   Theta0: Null
   Mutation rate: 5.3e-09
   Recombination rate: Null
   Time for generation: 31
   Custom filename: Null
   Lower bound: Null
   Upper bound: Null
   Parameter identifiers: Null
   Initial structure: [1, 1, 1]
   Final structure: [1, 1, 1]
   Only sudden: False
   Dynamics: [Sud, Lin, Exp]
   No migrations: False
   Symmetric migrations: False
   Migration masks: Null
   Selection: False
   Dominance: False
   Split fractions: False
   Inbreeding: False
   Ancestral size as parameter: False
   Upper bound of first split: Null
   Upper bound of second split: Null
   Local optimizer: BFGS_log
   Print models' code every N iteration: 0
   Model plot engine: moments
   Draw models every N iteration: 100
   Units of time in drawing: years
   Vmin: 1
   Silence: False
   Verbose: 1
   Number of repeats: 2
   Number of processes: 2
   Resume from: Null
   Only models: False
2. Inference with five populations: a) "You need to design and specify your own models ": Thats the first think I an confused with how can I design my own model for possible demographic histories? Just for learning purposes I used the same demographic models and ran it with 4 pop.
   b) I tried my best to use BO and the parameter file for BO run with 200 eval is here:
   cat bo_optimise_param_file
   Input file: /home/Gadma_Test/chr1_Merge.vcf,
   /home/Gadma_Test/4_pop_file
   Output directory: BO_combination
   Custom filename: demographic_model.py

Num init const: 2
Global optimizer: SMAC_BO_combination
Global maxeval: 200

Number of repeats: 64
Number of processes: 4

If you want to see the log files for both the run. Please let me know will share via email (70MB file can't upload here).

[noscode]

Hi @SachKrVinni,

Okay, it is more clear now.

1. Do you want to compare models with different structures? If you have just this run and structure [1,1,1], then you can take the best by log-likelohood result and that will be it.

2. Model specification might be challenging. Unfortunately, authors of moments recently deleted information about it from the manual. I can suggest several examples: for example, this one or your output of GADMA for three populations — GADMA generates code for the best model for all engines.