Structure v.2.3.4 is a standard tool for examining population genetic structure based on allele frequencies within and among populations. Although many new implementations of the structure algorithm have been developed in recent years offering improvements to speed, the classic tool offers a number of useful options that keep it relevant to this day.
# conda install ipyrad -c conda-forge -c bioconda
# conda install structure clumpp -c ipyrad
# conda install toyplot -c conda-forge
import ipyrad.analysis as ipa
import toyplot
Your input data should be a .snps.hdf
database file produced by ipyrad. If you do not have this you can generate it from any VCF file following the vcf2hdf5 tool tutorial. The database file contains the genotype calls information as well as linkage information that is used for subsampling unlinked SNPs and bootstrap resampling.
# the path to your .snps.hdf5 database file
SNPS = "/tmp/oaks.snps.hdf5"
# download example hdf5 dataset (158Mb, takes ~2-3 minutes)
URL = "https://www.dropbox.com/s/x6a4i47xqum27fo/virentes_ref.snps.hdf5?raw=1"
ipa.download(url=URL, path=SNPS);
file already exists
Structure infers the values of missing data while it runs the MCMC chain. No imputation is required, but it will perform more accurately if there is less missing data and when base calls are more accurate. For RAD-seq data I generally recommend not imputing data and simply filtering fairly stringently (e.g., <10-20% missing).
This example data set is fairly large in terms of the number of SNPs, but not that big in number of samples. Both factors affect run times. If you have many processors available you can parallelize runs so that many replicates across many K values can be run simultaneously. Still, expect run times to take many hours or even days.
The IMAP population assignments are used to select which samples from the dataset to include in the analysis, and also can be used to apply filtering to require that data is distributed across all relevant groups, using minmap settings. For more information on IMAP see the general ipa
tutorial.
# group individuals into populations
IMAP = {
"virg": ["LALC2", "SCCU3", "FLSF33", "FLBA140"],
"mini": ["FLSF47", "FLMO62", "FLSA185", "FLCK216"],
"gemi": ["FLCK18", "FLSF54", "FLWO6", "FLAB109"],
"sagr": ["CUVN10", "CUCA4", "CUSV6"],
"oleo": ["CRL0030", "HNDA09", "BZBB1", "MXSA3017"],
"fusi": ["MXED8", "MXGT4", "TXGR3", "TXMD3"],
"bran": ["BJSL25", "BJSB3", "BJVL19"],
}
# require that 50% of samples have data in each group
MINMAP = {i: 0.5 for i in IMAP}
Here we init a structure analysis tool instance. The name
and workdir
settings tell it where to save the output files. The imap
and minmap
settings are used to subselect samples and filter for missing data based on groups. The mincov
applies filters for missing data based on all samples. The minmaf
filter removes sites with a minor allele frequency cutoff.
When you init the object it will load the data and apply filtering. The printed output tells you how many SNPs were removed by each filter and the remaining number of missing data after filtering. By default the ipa.structure
tool samples a random set of unlinked SNPs for each replicate analysis, which represents a single SNP from each RAD locus. Thus, the replicate analyses in ipa.structure represent variation over different random samples of SNPs.
# init analysis object with input data and (optional) parameter options
struct = ipa.structure(
name="test",
workdir="analysis-structure",
data=SNPS,
imap=IMAP,
minmap=MINMAP,
mincov=0.95,
minmaf=0.05,
)
25 previous results loaded for run [test] Samples: 26 Sites before filtering: 1182005 Filtered (indels): 0 Filtered (bi-allel): 26588 Filtered (mincov): 958870 Filtered (minmap): 717351 Filtered (subsample invariant): 592354 Filtered (minor allele frequency): 56920 Filtered (combined): 1809076 Sites after filtering: 46962 Sites containing missing values: 25241 (53.75%) Missing values in SNP matrix: 25241 (2.07%) SNPs (total): 46962 SNPs (unlinked): 15887
The burnin
and numreps
parameters determine the length of the run. You can access other optional parameters of structure from the mainparams attribute.
struct.mainparams.burnin = 50000
struct.mainparams.numreps = 500000
See the general ipa tutorial for more details on parallelization. Here we select 5 K values to distribute tests over, and select 5 replicates for each test, for a total of 25 analyses. If you have 25 cores available then all jobs will start, else the extra jobs will be queued and start when the others finish. The progress bar tracks finished jobs, so it can jump from 0 to 50 or 100% in big strides.
struct.ipcluster['cores'] = 25
struct.run(nreps=5, kpop=[2, 3, 4, 5, 6], auto=True)
[####################] 100% 20:31:17 | running 25 structure jobs
The get_evanno_table
function will summarize finished results from your structure tool. (You can reload results later by re-initing a structure tool with the same name and workdir; see bottom of this notebook.) You can then use the plotting call below to visualize the deltaK value that is recommended for choosing the best K by the Evanno method. I'm of the camp that believes you shouldn't read into these model fitting exercises too much. Show results for multiple K values as seems fitting to describe diversity in your system.
There is an optional argument you can use called max_var_multiple
. This will filter out any replicate runs that have a estLnProbStdev that is N times greater than the replicate run with the smallest value. In other words, if some replicates have very large variance it may represent their failure to converge, and so they will be excluded. Try values like 2 or 3 for this option to explore its effect on your results.
etable = struct.get_evanno_table([2, 3, 4, 5])
etable
Nreps | lnPK | lnPPK | deltaK | estLnProbMean | estLnProbStdev | |
---|---|---|---|---|---|---|
2 | 5 | 0.000e+00 | 0.000e+00 | 0.000 | -2.984e+05 | 1.710e+03 |
3 | 5 | 3.355e+04 | 2.758e+04 | 46.695 | -2.648e+05 | 5.907e+02 |
4 | 5 | 5.966e+03 | 1.102e+06 | 4484.424 | -2.589e+05 | 2.458e+02 |
5 | 5 | -1.096e+06 | 0.000e+00 | 0.000 | -1.355e+06 | 2.450e+06 |
# get canvas object and set size
canvas = toyplot.Canvas(width=400, height=300)
# plot the mean log probability of the models in red
axes = canvas.cartesian(ylabel="estLnProbMean")
axes.plot(etable.estLnProbMean * -1, color="darkred", marker="o")
axes.y.spine.style = {"stroke": "darkred"}
# plot delta K with its own scale bar of left side and in blue
axes = axes.share("x", ylabel="deltaK", ymax=etable.deltaK.max() + etable.deltaK.max() * .25)
axes.plot(etable.deltaK, color="steelblue", marker="o");
axes.y.spine.style = {"stroke": "steelblue"}
# set x labels
axes.x.ticks.locator = toyplot.locator.Explicit(range(len(etable.index)), etable.index)
axes.x.label.text = "K (N ancestral populations)"
# choose a K value and run clumpp to permute results
k = 4
table = struct.get_clumpp_table(k)
[K4] 5/5 results permuted across replicates (max_var=0).
# sort list by columns
table.sort_values(by=list(range(k)), inplace=True)
# or, sort by a list of names (here taken from imap)
import itertools
onames = list(itertools.chain(*IMAP.values()))
table = table.loc[onames]
# build barplot
canvas = toyplot.Canvas(width=500, height=250)
axes = canvas.cartesian(bounds=("10%", "90%", "10%", "45%"))
axes.bars(table)
# add labels to x-axis
ticklabels = [i for i in table.index.tolist()]
axes.x.ticks.locator = toyplot.locator.Explicit(labels=ticklabels)
axes.x.ticks.labels.angle = -60
axes.x.ticks.show = True
axes.x.ticks.labels.offset = 10
axes.x.ticks.labels.style = {"font-size": "12px"}
Results files can be loaded by providing the name
and workdir
combination that leads to the path where your previous results were stored.
rerun = ipa.structure(
data=data,
name="test",
workdir="analysis-structure",
imap=imap,
load_only=True,
)
12 previous results loaded for run [test]
rerun.get_clumpp_table(3)
[K3] 3/3 results permuted across replicates (max_var=0).
0 | 1 | 2 | |
---|---|---|---|
BJSB3 | 0.0000 | 1.0000 | 0.0000 |
BJSL25 | 0.0000 | 1.0000 | 0.0000 |
BJVL19 | 0.0000 | 1.0000 | 0.0000 |
BZBB1 | 0.0000 | 0.0000 | 1.0000 |
CRL0030 | 0.0000 | 0.0000 | 1.0000 |
CUCA4 | 0.3450 | 0.0000 | 0.6550 |
CUSV6 | 0.4098 | 0.0000 | 0.5902 |
CUVN10 | 0.3408 | 0.0000 | 0.6592 |
FLAB109 | 1.0000 | 0.0000 | 0.0000 |
FLBA140 | 1.0000 | 0.0000 | 0.0000 |
FLCK18 | 1.0000 | 0.0000 | 0.0000 |
FLCK216 | 1.0000 | 0.0000 | 0.0000 |
FLMO62 | 0.9987 | 0.0010 | 0.0003 |
FLSA185 | 1.0000 | 0.0000 | 0.0000 |
FLSF33 | 1.0000 | 0.0000 | 0.0000 |
FLSF47 | 1.0000 | 0.0000 | 0.0000 |
FLSF54 | 1.0000 | 0.0000 | 0.0000 |
FLWO6 | 1.0000 | 0.0000 | 0.0000 |
HNDA09 | 0.0000 | 0.0000 | 1.0000 |
LALC2 | 1.0000 | 0.0000 | 0.0000 |
MXED8 | 0.1760 | 0.6953 | 0.1287 |
MXGT4 | 0.1531 | 0.8093 | 0.0377 |
MXSA3017 | 0.0477 | 0.0013 | 0.9510 |
SCCU3 | 1.0000 | 0.0000 | 0.0000 |
TXGR3 | 0.3649 | 0.6267 | 0.0083 |
TXMD3 | 0.3987 | 0.6010 | 0.0003 |
TXWV2 | 1.0000 | 0.0000 | 0.0000 |
You can continue an analysis with the same name
and workdir
by setting additional replicates or values of K and calling .run()
again. Here I will increase the number of replicates per K value from 3 to 5, and run one additional K value. Be sure to use all of the same parameter and filtering values that you used in the previous run or you might cause unexpected problems.
Here because we already finished 3 replicates for K=2,3,4,5 it will run 2 more for each of those, and it will run 5 replicates for K=6 since we do not have any finished replicates of those yet. You can see which result files exist for a named analysis object by accessing the .result_files
attribute, or by looking in the working directory. To overwrite existing files instead of adding more replicates you can use force=True
in the run command. You could also simply create a new object with a different name.
# init analysis object with same params as previously
struct = ipa.structure(
name="test",
data=data,
imap=imap,
minmap=minmap,
mincov=0.9,
)
# use the same params as before
struct.mainparams.burnin = 5000
struct.mainparams.numreps = 10000
# call run for all K values you want to have 5 finished replicates
struct.run(nreps=5, kpop=[2, 3, 4, 5, 6], auto=True)
12 previous results loaded for run [test] Samples: 27 Sites before filtering: 349914 Filtered (indels): 0 Filtered (bi-allel): 13001 Filtered (mincov): 222081 Filtered (minmap): 112898 Filtered (combined): 226418 Sites after filtering: 123496 Sites containing missing values: 96001 (77.74%) Missing values in SNP matrix: 142017 (4.26%) Parallel connection | oud: 4 cores [####################] 100% 3:39:43 | running 13 structure jobs
struct.get_evanno_table([2, 3, 4, 5, 6])
Nreps | lnPK | lnPPK | deltaK | estLnProbMean | estLnProbStdev | |
---|---|---|---|---|---|---|
2 | 5 | 0.00 | 0.00 | 0.000000 | -254950.12 | 1675.280397 |
3 | 5 | 25807.62 | 38146.76 | 39.701878 | -229142.50 | 960.830118 |
4 | 5 | -12339.14 | 7180.72 | 2.018760 | -241481.64 | 3556.995749 |
5 | 5 | -5158.42 | 8885.00 | 1.413531 | -246640.06 | 6285.676647 |
6 | 5 | 3726.58 | 0.00 | 0.000000 | -242913.48 | 2164.870641 |