{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# ipyrad-analysis toolkit: STRUCTURE\n", "\n", "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. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Required software" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# conda install ipyrad -c bioconda\n", "# conda install -c bioconda -c ipyrad structure clumpp\n", "# conda install toyplot -c eaton-lab" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import ipyrad.analysis as ipa\n", "import toyplot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Required input data files\n", "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](https://ipyrad.readthedocs.io/en/latest/API-analysis/cookbook-vcf2hdf5.html). The database file contains the genotype calls information as well as linkage information that is used for subsampling unlinked SNPs and bootstrap resampling." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# the path to your .snps.hdf5 database file\n", "data = \"/home/deren/Downloads/ref_pop2.snps.hdf5\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "### Note: missing data in STRUCTURE analyses:\n", "\n", "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. I recommend not imputing data and simply filtering fairly stringently. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Approximate run times\n", "This example data set should probably be run for a longer burnin and number of reps if it were to be used in a publication. For reference, this data set takes about 2.5 hours to run 12 jobs on a 4 core laptop for a data set with 27 samples and ~125K SNPs. If your data set has more samples or SNPs then it will take longer. If you have 2X as many cores then it will run 2X faster. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Input data file and population assignments\n", "If you are using the \"sample\" input method then population assignments (imap dictionary) are used for for filtering, color coding plots, and for imputation. If you are using the \"kmeans\" imputing method then population assignments are only used for filtering and color coding plots." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# group individuals into populations\n", "imap = {\n", " \"virg\": [\"TXWV2\", \"LALC2\", \"SCCU3\", \"FLSF33\", \"FLBA140\"],\n", " \"mini\": [\"FLSF47\", \"FLMO62\", \"FLSA185\", \"FLCK216\"],\n", " \"gemi\": [\"FLCK18\", \"FLSF54\", \"FLWO6\", \"FLAB109\"],\n", " \"sagr\": [\"CUVN10\", \"CUCA4\", \"CUSV6\"],\n", " \"oleo\": [\"CRL0030\", \"HNDA09\", \"BZBB1\", \"MXSA3017\"],\n", " \"fusi\": [\"MXED8\", \"MXGT4\", \"TXGR3\", \"TXMD3\"],\n", " \"bran\": [\"BJSL25\", \"BJSB3\", \"BJVL19\"],\n", "}\n", "\n", "# require that 50% of samples have data in each group\n", "minmap = {i: 0.5 for i in imap}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Enter data file and params\n", "The `struct` analysis object takes input data as the *.snps.hdf5* file produced by ipyrad. All other parameters are optional. The **imap** dictionary groups individuals into populations and **minmap** can be used to filter SNPs to only include those that have data for at least some proportion of samples in every group. The **mincov** option works similarly, it filters SNPs that are shared across less than some proportion of all samples (in contrast to minmap this does not use imap groupings). \n", "\n", "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 amount of missing data after filtering. These remaining missing values are the ones that will be filled with imputation. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Samples: 27\n", "Sites before filtering: 349914\n", "Filtered (indels): 0\n", "Filtered (bi-allel): 13001\n", "Filtered (mincov): 222081\n", "Filtered (minmap): 112898\n", "Filtered (combined): 226418\n", "Sites after filtering: 123496\n", "Sites containing missing values: 96001 (77.74%)\n", "Missing values in SNP matrix: 142017 (4.26%)\n" ] } ], "source": [ "# init analysis object with input data and (optional) parameter options\n", "struct = ipa.structure(\n", " name=\"test\",\n", " data=data,\n", " imap=imap,\n", " minmap=minmap,\n", " mincov=0.9,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Run STRUCTURE and plot results. \n", "The `burnin` and `numreps` parameters determine the length of the run. For analyses with many samples and with larger values of K you should use much larger values than these. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "struct.mainparams.burnin = 5000\n", "struct.mainparams.numreps = 10000" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parallel connection | oud: 4 cores\n", "[####################] 100% 2:26:57 | running 12 structure jobs \n" ] } ], "source": [ "struct.run(nreps=3, kpop=[2, 3, 4, 5], auto=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Analyze results: Choosing K" ] }, { "cell_type": "code", "execution_count": 83, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
NrepslnPKlnPPKdeltaKestLnProbMeanestLnProbStdev
230.0000000.0000000.000000-254535.7666672023.420259
3325229.90000035261.16666730.892635-229305.8666671141.410147
43-10031.2666671451.8000001.675614-239337.133333866.428568
53-8579.4666670.0000000.000000-247916.6000008537.460208
\n", "
" ], "text/plain": [ " Nreps lnPK lnPPK deltaK estLnProbMean estLnProbStdev\n", "2 3 0.000000 0.000000 0.000000 -254535.766667 2023.420259\n", "3 3 25229.900000 35261.166667 30.892635 -229305.866667 1141.410147\n", "4 3 -10031.266667 1451.800000 1.675614 -239337.133333 866.428568\n", "5 3 -8579.466667 0.000000 0.000000 -247916.600000 8537.460208" ] }, "execution_count": 83, "metadata": {}, "output_type": "execute_result" } ], "source": [ "etable = struct.get_evanno_table([2, 3, 4, 5])\n", "etable" ] }, { "cell_type": "code", "execution_count": 90, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
2345K (N ancestral populations)230000235000240000245000250000255000estLnProbMean010203040deltaK
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# get canvas object and set size\n", "canvas = toyplot.Canvas(width=400, height=300)\n", "\n", "# plot the mean log probability of the models in red\n", "axes = canvas.cartesian(ylabel=\"estLnProbMean\")\n", "axes.plot(etable.estLnProbMean * -1, color=\"darkred\", marker=\"o\")\n", "axes.y.spine.style = {\"stroke\": \"darkred\"}\n", "\n", "# plot delta K with its own scale bar of left side and in blue\n", "axes = axes.share(\"x\", ylabel=\"deltaK\", ymax=etable.deltaK.max() + etable.deltaK.max() * .25)\n", "axes.plot(etable.deltaK, color=\"steelblue\", marker=\"o\");\n", "axes.y.spine.style = {\"stroke\": \"steelblue\"}\n", "\n", "# set x labels\n", "axes.x.ticks.locator = toyplot.locator.Explicit(range(len(etable.index)), etable.index)\n", "axes.x.label.text = \"K (N ancestral populations)\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Analyze results: Barplots" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[K3] 3/3 results permuted across replicates (max_var=0).\n" ] } ], "source": [ "k = 3\n", "table = struct.get_clumpp_table(k)\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# sort list by columns\n", "table.sort_values(by=list(range(k)), inplace=True)\n", "\n", "# or, sort by a list of names (here taken from imap)\n", "import itertools\n", "onames = list(itertools.chain(*imap.values()))\n", "table = table.loc[onames]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
TXWV2LALC2SCCU3FLSF33FLBA140FLSF47FLMO62FLSA185FLCK216FLCK18FLSF54FLWO6FLAB109CUVN10CUCA4CUSV6CRL0030HNDA09BZBB1MXSA3017MXED8MXGT4TXGR3TXMD3BJSL25BJSB3BJVL190.00.51.0
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# build barplot\n", "canvas = toyplot.Canvas(width=500, height=250)\n", "axes = canvas.cartesian(bounds=(\"10%\", \"90%\", \"10%\", \"45%\"))\n", "axes.bars(table)\n", "\n", "# add labels to x-axis\n", "ticklabels = [i for i in table.index.tolist()]\n", "axes.x.ticks.locator = toyplot.locator.Explicit(labels=ticklabels)\n", "axes.x.ticks.labels.angle = -60\n", "axes.x.ticks.show = True\n", "axes.x.ticks.labels.offset = 10\n", "axes.x.ticks.labels.style = {\"font-size\": \"12px\"}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Cookbook" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Advanced: Load existing results

\n", "\n", "Results files can be loaded by providing the `name` and `workdir` combination that leads to the path where your previous results were stored. " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "12 previous results loaded for run [test]\n" ] } ], "source": [ "rerun = ipa.structure(\n", " data=data, \n", " name=\"test\", \n", " workdir=\"analysis-structure\",\n", " imap=imap,\n", " load_only=True,\n", ")" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[K3] 3/3 results permuted across replicates (max_var=0).\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
012
BJSB30.00001.00000.0000
BJSL250.00001.00000.0000
BJVL190.00001.00000.0000
BZBB10.00000.00001.0000
CRL00300.00000.00001.0000
CUCA40.34500.00000.6550
CUSV60.40980.00000.5902
CUVN100.34080.00000.6592
FLAB1091.00000.00000.0000
FLBA1401.00000.00000.0000
FLCK181.00000.00000.0000
FLCK2161.00000.00000.0000
FLMO620.99870.00100.0003
FLSA1851.00000.00000.0000
FLSF331.00000.00000.0000
FLSF471.00000.00000.0000
FLSF541.00000.00000.0000
FLWO61.00000.00000.0000
HNDA090.00000.00001.0000
LALC21.00000.00000.0000
MXED80.17600.69530.1287
MXGT40.15310.80930.0377
MXSA30170.04770.00130.9510
SCCU31.00000.00000.0000
TXGR30.36490.62670.0083
TXMD30.39870.60100.0003
TXWV21.00000.00000.0000
\n", "
" ], "text/plain": [ " 0 1 2\n", "BJSB3 0.0000 1.0000 0.0000\n", "BJSL25 0.0000 1.0000 0.0000\n", "BJVL19 0.0000 1.0000 0.0000\n", "BZBB1 0.0000 0.0000 1.0000\n", "CRL0030 0.0000 0.0000 1.0000\n", "CUCA4 0.3450 0.0000 0.6550\n", "CUSV6 0.4098 0.0000 0.5902\n", "CUVN10 0.3408 0.0000 0.6592\n", "FLAB109 1.0000 0.0000 0.0000\n", "FLBA140 1.0000 0.0000 0.0000\n", "FLCK18 1.0000 0.0000 0.0000\n", "FLCK216 1.0000 0.0000 0.0000\n", "FLMO62 0.9987 0.0010 0.0003\n", "FLSA185 1.0000 0.0000 0.0000\n", "FLSF33 1.0000 0.0000 0.0000\n", "FLSF47 1.0000 0.0000 0.0000\n", "FLSF54 1.0000 0.0000 0.0000\n", "FLWO6 1.0000 0.0000 0.0000\n", "HNDA09 0.0000 0.0000 1.0000\n", "LALC2 1.0000 0.0000 0.0000\n", "MXED8 0.1760 0.6953 0.1287\n", "MXGT4 0.1531 0.8093 0.0377\n", "MXSA3017 0.0477 0.0013 0.9510\n", "SCCU3 1.0000 0.0000 0.0000\n", "TXGR3 0.3649 0.6267 0.0083\n", "TXMD3 0.3987 0.6010 0.0003\n", "TXWV2 1.0000 0.0000 0.0000" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rerun.get_clumpp_table(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Advanced: Add replicates or additional K values

\n", "\n", "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. \n", "\n", "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. " ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "12 previous results loaded for run [test]\n", "Samples: 27\n", "Sites before filtering: 349914\n", "Filtered (indels): 0\n", "Filtered (bi-allel): 13001\n", "Filtered (mincov): 222081\n", "Filtered (minmap): 112898\n", "Filtered (combined): 226418\n", "Sites after filtering: 123496\n", "Sites containing missing values: 96001 (77.74%)\n", "Missing values in SNP matrix: 142017 (4.26%)\n", "Parallel connection | oud: 4 cores\n", "[####################] 100% 3:39:43 | running 13 structure jobs \n" ] } ], "source": [ "# init analysis object with same params as previously\n", "struct = ipa.structure(\n", " name=\"test\",\n", " data=data,\n", " imap=imap,\n", " minmap=minmap,\n", " mincov=0.9,\n", ")\n", "\n", "# use the same params as before \n", "struct.mainparams.burnin = 5000\n", "struct.mainparams.numreps = 10000\n", "\n", "# call run for all K values you want to have 5 finished replicates\n", "struct.run(nreps=5, kpop=[2, 3, 4, 5, 6], auto=True)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
NrepslnPKlnPPKdeltaKestLnProbMeanestLnProbStdev
250.000.000.000000-254950.121675.280397
3525807.6238146.7639.701878-229142.50960.830118
45-12339.147180.722.018760-241481.643556.995749
55-5158.428885.001.413531-246640.066285.676647
653726.580.000.000000-242913.482164.870641
\n", "
" ], "text/plain": [ " Nreps lnPK lnPPK deltaK estLnProbMean estLnProbStdev\n", "2 5 0.00 0.00 0.000000 -254950.12 1675.280397\n", "3 5 25807.62 38146.76 39.701878 -229142.50 960.830118\n", "4 5 -12339.14 7180.72 2.018760 -241481.64 3556.995749\n", "5 5 -5158.42 8885.00 1.413531 -246640.06 6285.676647\n", "6 5 3726.58 0.00 0.000000 -242913.48 2164.870641" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "struct.get_evanno_table([2, 3, 4, 5, 6])" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }