Determine burnin and number of MCMC replicates - K clusters?

[Chiara Duijser]

Hi everyone,

I have a question about inferring the number of genetic cluster present in my dataset using STRUCTURE in R. I have sampled individuals across 4 habitats and want to see how the populations from those habitats are related based on ancestry. After filtering of my dataset I end up with 1930 SNPs. I am using the dartR packages available since my samples are sequenced by Dartseq. 

I have run my analysis using the following code:

out_struc <- gl.run.structure(

  gl5, 

  k.range=1:5, 

  num.k.rep=10, 

  burnin = 1000, [default]

  numreps = 1000, [default]

  exec="./structure/structure", 

  noadmix=FALSE)

Using gl.evanno like below I get the graphs on meanLnP(K) and deltaK that I want to use to determine the correct number of clusters present (see graphs attached)

out_evanno <- gl.evanno(out_struc, plot.out=TRUE) 

However, I am not entirely sure if these graphs are very informative or if I have to increase the burnin since meanLnP(K) is highest around 3/4 but my deltaK is highest around 2 so that confuses me. 

I am currently running the same line of code wil a burnin of 20000 and numreps of 10000 but analyzing this takes a very long time with my computer and I’m still waiting to see these results. Are there any recommendations on what would be the best way to determine the best number of iterations for MCMC burnin (burnin) and the number of MCMC replicates (numreps) for my dataset? I have chosen K 1:5 since I have 4 different habitats/four potentially different populations but read that it's best to go to a higher K than 4 in that case?

Thank you very much for your help!

Kind regards,

Chiara

[Julianna Santos]

Hi Chiara and everyone

Chiara, I do not have an answer to your question, but I'm interested in following that up.

Did you find different results for burnin of 20000 and 10000 reps?

I ran my dataset with

74 individuals, 5475 loci
k = 1:15

10 replicates

First, burnin 20000 and numreps 1000.

Then, I ran it a second time, burnin 200000 and  numreps 100000

I got very different results (see attached -- first run named as 10^4 and second as 10^5).

So my question is, should I always use the results from the maximum number of runs?

Cheers,

[Kate Moffatt]

Hi Julianna, 

I was just wondering how long it took your data set to run in Structure? I have a similar size data set. 73 individuals, 4000 loci and a population of 6. I haven't run my data yet as I am wanting to know if this takes hours or days? Especially for the burnin 200000 and numreps 100000

Could you please let me know how long your analysis took to run? 

Sorry I cannot help you with your question. 

Thanks, 

Kate 

[Julianna Santos]

Hi

I've just realised I provided the wrong parameters in the previous message.

To clarify:

First run:

10,000 Burn-in period
20,000 Reps

Second run:

100000 Burn-in period
200000 Reps

Cheers,

Julianna.

--
You received this message because you are subscribed to a topic in the Google Groups "structure-software" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/structure-software/J4PYKhTpVo0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to structure-softw...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/structure-software/0c560e43-0e13-4f76-9cef-cf232bebe147n%40googlegroups.com.

[Julianna Santos]

Hi Kate

I'm running it on a high-performance computer.

The first run (10,000 Burn-in period and 20,000 Reps) was fairly fast: it took ~2-3 hours. I also ran this (same number of interactions, etc.) from the front end on a remote desktop (16 cores and 64GB RA). It took five days, but at least it did not crash.

The second run (100,000 Burn-in period and 200,000 Rep) was completed in 11 hours on an HPC (I did not try it on a regular computer).

I have another dataset -- substantially larger in terms of loci (>18.000) -- but 62 individuals. It is now completing 24 hours on an HPC.

Happy to answer other questions :)

Cheers,

Julianna.

Julianna Santos | Conservation Biologist

MSc - Zoology

MBA - Project Management

PhD Candidate at The University of Melbourne, Australia

Webpage: https://biodynamos.science.unimelb.edu.au/

M +61 0410 317 628

--
You received this message because you are subscribed to the Google Groups "structure-software" group.
To unsubscribe from this group and stop receiving emails from it, send an email to structure-softw...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/structure-software/d1ed7e57-615a-462a-9017-fb8243ce68c5n%40googlegroups.com.

[Chiara Duijser]

Hi Julianna and Kate,

Thank you for your replies! I have been running it on a normal Windows desktop and it took around 24h for the burnin of 20000 and numreps of 10000 and k.range = 1:5, num.k.rep = 10. I do find different results between low and high burnin and MCMC replicates. I'm not too sure if you should always use the results from the maximum number of runs but would assume it's most accurate due to replication? 

Can I ask you how you determined your burnin and k number to test? Is that based on similar papers or a guideline?

Thank you!

Cheers,

Chiara

Op vrijdag 23 juni 2023 om 11:29:46 UTC+10 schreef Julianna Santos:

[Chiara Duijser]

Hi Julianna,

In the original paper by Pritchard et al. 2000  Inference of Population Structure Using Multilocus Genotype Data it does say "In general, substantial differences between runs can indicate that either the runs should be longer to obtain more accurate estimates or that independent runs are getting stuck in different modes in the parameter space". It might be worth having a look at this if you haven't read the paper already! I'm currently running a new analysis with k.range = 1:5, num.k.rep = 10, burnin = 30000, and numreps = 1000000 to see if anything changes compared to my previous run.

Cheers,

Chiara

Op dinsdag 27 juni 2023 om 11:10:28 UTC+10 schreef Chiara Duijser:

[Julianna Santos]

Hi Chiara

Thank you for following up on this. Yes, I had read that in the paper.

I've got relatively consistent results between runs of each K (num.k.rep = 10) for both -- burnin 10^4 and 10^5 (numreps 2x10^4 and 2x10^5, respectively). Even though 10^5 is, in fact, slightly more consistent for the summary statistics, the difference is not that big. Even comparing the bar plots for each run, both (10^4 or 10^5) produce similar plots among the runs for each K.

Having said that, my question in the first place was more related to the Evanno method results: the optimum K for 10^5 (K = 5 or K = 11) doesn't make sense biologically (note the evanno method figures are misnamed in the first email -- attached the correct versions). The optimum K for 10^4 (K=2 or K=4) does make sense!

So, I would consider the results from the higher burnin and replicates (10^5), but what do I do about the optimum K then?

Ps: answering your previous question: I determined burnin and K based on previous studies and research group standard practices -- note this is the first time I am working with genetics.

Cheers,

You received this message because you are subscribed to a topic in the Google Groups "structure-software" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/structure-software/J4PYKhTpVo0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to structure-softw...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/structure-software/0a5c32ec-4331-4175-8f6c-4a3fd95d78fdn%40googlegroups.com.

[Chiara Duijser]

Hi Julianna,

I looked at your Evanno plots and see what you mean now! It is also the first time for me working with population genetic data and am afraid I don't have the answer for you, sorry! 

Personally I would think it is also important to check if the results make sense biologically. There are R packages to determine the optimal number of K like Structure Harvester but I haven't used these yet. Maybe that could help?

Cheers,

Chiara

Op woensdag 28 juni 2023 om 11:29:14 UTC+10 schreef Julianna Santos:

[Vikram Chhatre]

You can always try using other methods to corroborate your results. For example PCA based methods, Admixture, FastSTR etc. 

Also, appropriate burnin length can be determined based on whether the alpha parameter has stabilized. The manual talks at length about this.

To view this discussion on the web visit https://groups.google.com/d/msgid/structure-software/cc65c6d0-78a2-4554-aae1-d1fa4902599fn%40googlegroups.com.

[Chiara Duijser]

Thank you Vikram for the extra information!

I do have a very basic question, but have been struggling finding correct structureRun output files and using pophelper. @Julianna, looking at your Evanno plots I think you have used pophelper as well? 

When my structureRun is finished I get the following output files (see attached schroonshot; the list goes on but this is just to show you). However, for the online interactive version of pophelper (http://pophelper.com) or for the R function for that matter I need a .txt file. Additionally, I get the following information in my R console that I copy-pasted here:

MCMC completed

 

Inferred ancestry of individuals:

        Label (%Miss) Pop:  Inferred clusters

  1     19.1    (0)    1 :  0.997 0.001 0.001 0.001 0.001 

  2       20    (0)    1 :  0.002 0.002 0.994 0.001 0.001 

  3       21    (0)    1 :  0.002 0.002 0.994 0.001 0.001 

  4     22.1    (0)    1 :  0.997 0.001 0.001 0.001 0.001 

  5       23    (0)    1 :  0.997 0.001 0.001 0.001 0.001 

  6       24    (0)    1 :  0.997 0.001 0.001 0.001 0.001 

  7     25.1    (0)    3 :  0.997 0.001 0.001 0.001 0.001 

  8       26    (0)    3 :  0.997 0.001 0.001 0.001 0.001 

  9     26.1    (0)    3 :  0.997 0.001 0.001 0.001 0.001 

 10       28    (0)    3 :  0.997 0.001 0.001 0.001 0.001 

 11       29    (0)    3 :  0.997 0.001 0.001 0.001 0.001 

 12       30    (0)    2 :  0.001 0.001 0.001 0.001 0.997 

 13     30.1    (0)    3 :  0.997 0.001 0.001 0.001 0.001 

 14     31.1    (0)    2 :  0.001 0.001 0.001 0.001 0.997 

 15     32.1    (0)    2 :  0.001 0.001 0.001 0.001 0.997 

 16       34    (0)    2 :  0.001 0.001 0.001 0.001 0.997 

 17       36    (0)    2 :  0.013 0.351 0.012 0.005 0.618 

 18     36.1    (0)    4 :  0.062 0.719 0.057 0.006 0.156 

 19     37.1    (0)    4 :  0.144 0.726 0.088 0.004 0.039 

 20     38.1    (0)    4 :  0.021 0.824 0.139 0.005 0.012 

 21       39    (0)    4 :  0.001 0.001 0.001 0.001 0.997 

 22     40.1    (0)    4 :  0.098 0.681 0.033 0.006 0.182 

 23       41    (0)    4 :  0.001 0.001 0.001 0.001 0.996 

 

Command line arguments:   C:\Installers\structure_windows_console\console\structure.exe -m gtypes.created.on.2023-06-15.11.38.57.structureRun/gtypes.created.on.2023-06-15.11.38.57.structureRun.k5.r10_mainparams -e gtypes.created.on.2023-06-15.11.38.57.structureRun/gtypes.created.on.2023-06-15.11.38.57.structureRun.k5.r10_extraparams -i gtypes.created.on.2023-06-15.11.38.57.structureRun/gtypes.created.on.2023-06-15.11.38.57.structureRun.k5.r10_data -o gtypes.created.on.2023-06-15.11.38.57.structureRun/gtypes.created.on.2023-06-15.11.38.57.structureRun.k5.r10_out

 

Input File:    gtypes.created.on.2023-06-15.11.38.57.structureRun/gtypes.created.on.2023-06-15.11.38.57.structureRun.k5.r10_data

 

Output File:   gtypes.created.on.2023-06-15.11.38.57.structureRun/gtypes.created.on.2023-06-15.11.38.57.structureRun.k5.r10_out_f

 

Run parameters:

   23 individuals

   1930 loci

   5 populations assumed

   20000 Burn-in period

   100000 Reps

 

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

Estimated Ln Prob of Data   = -27501.2

Mean value of ln likelihood = -24743.4

Variance of ln likelihood   = 5515.7

Mean value of alpha         = 1.0000

Allele frequencies uncorrelated

 

 

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

Proportion of membership of each pre-defined

 population in each of the 5 clusters

 

Given    Inferred Clusters                     Number of

 Pop       1      2      3      4      5      Individuals

 

  1:     0.665  0.001  0.332  0.001  0.001        6

  2:     0.003  0.071  0.003  0.001  0.921        5

  3:     0.997  0.001  0.001  0.001  0.001        6

  4:     0.055  0.492  0.053  0.004  0.397        6

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

Final results printed to file gtypes.created.on.2023-06-15.11.38.57.structureRun/gtypes.created.on.2023-06-15.11.38.57.structureRun.k5.r10_out_f

 

Completed: gl.run.structure 

However, when I translate this to a .txt file it does not work and none of the files in the folder are .txt files either. How can I access the correct output files for using pophelper?

Thank you in advance for any information!

Kind regards,

Chiara 

Op donderdag 29 juni 2023 om 12:09:17 UTC+10 schreef Vikram Chhatre:

[Josh Banta]

Dear Chiara,

I have a tutorial for this:

https://youtu.be/HJgJ4fVJq2s

The tutorial uses the Ubuntu Linux operating system, but any operating system will work.

Best,

Josh

Chiara 

To view this discussion on the web visit https://groups.google.com/d/msgid/structure-software/260bec20-0977-4a95-ac05-add475946a02n%40googlegroups.com.