An Introduction to the gstudio Package

Import from GENEPOP files

The GENEPOP file format, as interpreted by gstudio has the following format.

1. It is a text file and preferrably space delimited.

• The first line has some info for you, it is mostly ignored but should include a bit about what the data are.
• The next \(K\) lines of the file list the names of the loci to be used, one per line.
• The rest of the file contains populations. Each population starts with “Pop” alone on its own line. All individuals are assumed to be in the same population until such time the next “Pop” is reached.
• For each individual, there is an identification value first, and followed by a comma are the genotypes.
• All genotypes are encoded using 3 digits for each allele (e.g., a 3:5 genotype is 003005). Missing data is all zeros ‘000000’ and haploid is just three digits.

The import from read_population() names loci, adds “ID” for the identification column and adds to each population a “Pop-X” designation. Other than that, it is identical to what is below. However, if you have your data in a spreadsheet, there is no need to shove it into a genepop format to import into R.

Text File Import

Options for this function include:

Here are some examples of data files with different kinds of genetic data, each of which exercises the read_population() function in a different way. Hopefully this covers the main types of data being imported, if not, drop me an email Rodney Dyer. Missing genotypes should be missing data or encoded as NA. If you do not have a genotype then leave it blank. There is no reason to use negative numbers or other conventions.

Columns of genotypes are indicated by the required parameter locus.columns so that read_population() knows which columns to treat as locus objects and which to leave as normal data for R. Without this parameter, the data will be read in as character or numeric.

There are some example data files included in the project for you to look at. Depending upon how your computer is set up, they may be placed in different locations. Here is a quick way to find out where the installed folder is for the gstudio package and the location of the ‘data’ folder within it.

system.file("extdata", package = "gstudio")

## [1] "/Library/Frameworks/R.framework/Versions/3.1/Resources/library/gstudio/extdata"

Two Column Data

Here is an example of data where the genotypes are encoded as two columns of data in a csv file.

file <- system.file("extdata", "data_2_column.csv", package = "gstudio")
data <- read_population(file, type = "column", locus.columns = 4:7)
data

##   Population Latitude Longitude  Loc1  Loc2
## 1    Olympia    47.15   -122.89 12:14 15:15
## 2 Bellingham    48.75   -122.49 12:12 15:17
## 3  St. Louis    38.81    -89.98 10:12      
## 4       Ames    42.26    -93.47 10:14      
## 5   Richmond    37.74    -77.16 10:10 17:17

Phased Data

There are times when the gametic phase of the genotypes is important. By default, gstudio will keep alleles sorted in alpha/numeric order. If you need to keep this from happening, pass the optional phased option to read_population(). Notice the differences between this and the previous genotypes.

file <- system.file("extdata", "data_2_column.csv", package = "gstudio")
data <- read_population(file, type = "column", locus.columns = 4:7, phased = TRUE)
data

##   Population Latitude Longitude  Loc1  Loc2
## 1    Olympia    47.15   -122.89 14:12 15:15
## 2 Bellingham    48.75   -122.49 12:12 17:15
## 3  St. Louis    38.81    -89.98 12:10      
## 4       Ames    42.26    -93.47 10:14      
## 5   Richmond    37.74    -77.16 10:10 17:17

AFLP-like data

Genotypes that are ‘aflp’-like are encoded as binary characters (e.g., 0/1) indicating the presence or absence of a particular band.

file <- system.file("extdata", "data_aflp.csv", package = "gstudio")
data <- read_population(file, type = "aflp", locus.columns = c(4, 5))
data

##   Population Latitude Longitude Loc1 Loc2
## 1    Olympia    47.15   -122.89    1    1
## 2 Bellingham    48.75   -122.49    0    0
## 3  St. Louis    38.81    -89.98    1    0
## 4       Ames    42.26    -93.47    0    1
## 5   Richmond    37.74    -77.16    0    0

SNP Minor Allele Data

At times, SNP data is encoded in relation to the number of minor alleles. You can import these data using the type=“snp” option and it will encode them as ‘AA’, ‘AB’, or ‘BB’ with the ‘B’ allele as the minor one.

file <- system.file("extdata", "data_snp.csv", package = "gstudio")
data <- read_population(file, type = "snp", locus.columns = 4:7)
data

##   Population Latitude Longitude snp1 snp2 snp3 snp4
## 1    Olympia    47.15   -122.89  A:B  B:B  A:B  A:B
## 2 Bellingham    48.75   -122.49  A:B  A:A  A:B  B:B
## 3  St. Louis    38.81    -89.98  A:B  A:A          
## 4       Ames    42.26    -93.47  A:A  A:B          
## 5   Richmond    37.74    -77.16  A:A  B:B  A:B  B:B

Zyme-Like Data

Some data is encoded as allozyme genotypes (e.g., 33, 35, 55 for diploid individuals with alleles ‘3’ and ‘5’).

file <- system.file("extdata", "data_zymelike.csv", package = "gstudio")
data <- read_population(file, type = "zyme", locus.columns = 4:7)
data

##   Population Latitude Longitude Loc1 Loc2 Loc3 Loc4
## 1    Olympia    47.15   -122.89  1:2  1:4  1:5  1:5
## 2 Bellingham    48.75   -122.49  1:2  1:2  1:5  1:7
## 3  St. Louis    38.81    -89.98  1:2  1:3          
## 4       Ames    42.26    -93.47  1:1  1:4          
## 5   Richmond    37.74    -77.16  1:2  1:1  2:2  2:2

Pre-Separated Data For Higher Ploidy

file <- system.file("extdata", "data_separated.csv", package = "gstudio")
data <- read_population(file, type = "separated", locus.columns = c(4, 5))
data

##   Population Latitude Longitude  Loc1  Loc2
## 1    Olympia    47.15   -122.89 12:14 15:15
## 2 Bellingham    48.75   -122.49 12:12 15:17
## 3  St. Louis    38.81    -89.98 12:10      
## 4       Ames    42.26    -93.47 10:14      
## 5   Richmond    37.74    -77.16 10:10 17:17

Raw R objects

Saving data once it is in R is trivial and you do it as you would for any other R object. The R object system knows how to serialize its own data using the load() and save() functions.

save(df, file = "MyData.rda")

To load the objects back into the work space, you just do:

load("MyData.rda")

And you can verify that you have data in your work space by listing it.

ls()

##  [1] "data"      "df"        "file"      "loc0"      "loc1"     
##  [6] "loc2"      "loc3"      "loc4"      "loc5"      "loc6"     
## [11] "loci"      "metadata"  "tidy.opts" "x"

Saving as Text

As a default, the function write_population() will write your data file as a comma separated text file with the loci encoded as column separated (see type="separated" above).

write_population(df, file = "~/Desktop/MyData.csv")

Saving as GENEPOP

Raw data can be saved in GENEPOP formats by passing an optional argument mode="genepop" to the write_population() function.

write_population(df, file = "~/Desktop/MyData.txt", mode = "genepop")

Saving for STRUCTURE

Raw data can also be exported for analysis using the program STRUCTURE. Here the optional argument is type="structure"

write_population(df, file = "~/Desktop/MyData.str", mode = "structure")

Data Classes

Often it is important to know which columns of a data set are actually of a particular data type. Here is a simple function that tells you either the name or index of columns in a data.frame that of a specific data type.

column_class(arapat, class = "locus")

## [1] "LTRS" "WNT"  "EN"   "EF"   "ZMP"  "AML"  "ATPS" "MP20"

column_class(arapat, class = "locus", mode = "index")

## [1]  7  8  9 10 11 12 13 14

–

Rarefaction

Rarefaction is a technique used to measure diversity in different populations. It is particularly important for situations where you have different sample sizes. Is there more diversity in the larger population because you sampled more or is it a truly more diverse population? I’ll use the data from the beetle to show how diversity changes with sample sizes and highlight how you can use the rarefaction() function.

In the mainland populations, there are only 36 samples and the allelic diversity is relatively low at the WNT locus.

loci.son <- arapat$WNT[arapat$Species == "Mainland"]
length(loci.son)

## [1] 36

genetic_diversity(loci.son, mode = "Ae")

##      Ae
## 1 1.034

The larger clade on the peninsula has many more individuals and is more diverse.

loci.peninsula <- arapat$WNT[arapat$Species == "Peninsula"]
length(loci.peninsula)

## [1] 252

genetic_diversity(loci.peninsula, mode = "Ae")

##      Ae
## 1 2.025

So is this difference a consequence of the sample sizes or is peninsular Baja California more genetically diverse? To answer this, rarefaction randomly sub-samples the loci.baja data and estimates the value of \(\hat{A}_e\) for samples of size 36 (for our case) to see if the observed diversity differences are due to sampling alone. To visualize the distribution, I throw it into a data.frame and use the ggplot2 functions to make the pretty colored histogram.

Ae.peninsula <- rarefaction(loci.peninsula, mode = "Ae", size = 36)
df <- data.frame(Ae.peninsula)
ggplot(df, aes(x = Ae.peninsula)) + geom_histogram(aes(fill = ..count..), binwidth = 0.05) + 
    scale_fill_gradient("count", low = "#cccccc", high = "#a60000") + theme_bw()

AMOVA Distance

The AMOVA distance metric was first introduced in Excoffier et al. (1992) using restriction fragment encodings. However, a more elegant description of it was described in Smouse & Peakall (1999) using a geometric interpretation. Essentially, the coding of alleles at a locus can be depicted by a vector \(\vec{p}\) whose length is equal to the number of alleles at the locus. The presence of an allele increments the appropriate element of that vector. For two individuals, the squared vector distance between them is:

\[ \delta_{ij}^2 = 2(p_i-p_j)^2 \]

Across loci, these are additive (though see Smouse & Peakall for additional weighing schemes by locus or allele and the relative power of adopting such approaches). Across a set of individuals the squared genetic distance between all individuals can be represented by a square symmetric matrix with a zero diagonal, D. The AMOVA analysis itself is conducting by taking the “Sums of Squared Genetic Distances” for all individuals (SSGD(Total)), within group SSW, and among groups SSA and variance is decomposed following the standard approach for a random effects model. There is essentially no mystery in this approach, though it has been shrouded in obscurity. Dyer et al. (2004) showed how this is just a multivariate linear model amenable to a much broader range of experimental designs than just 1-way and nested designs.

D <- dist_amova(loci)
rownames(D) <- colnames(D) <- as.character(loci)
D

##     A:A A:B A:C A:D B:B B:C B:D C:C C:D D:D
## A:A   0   1   1   1   4   3   3   4   3   4
## A:B   1   0   1   1   1   1   1   3   2   3
## A:C   1   1   0   1   3   1   2   1   1   3
## A:D   1   1   1   0   3   2   1   3   1   1
## B:B   4   1   3   3   0   1   1   4   3   4
## B:C   3   1   1   2   1   0   1   1   1   3
## B:D   3   1   2   1   1   1   0   3   1   1
## C:C   4   3   1   3   4   1   3   0   1   4
## C:D   3   2   1   1   3   1   1   1   0   1
## D:D   4   3   3   1   4   3   1   4   1   0

Bray Curtis Individual Distance

Bray-Curtis Distance (Bray & Curtis 1957) has been primarily used to quantify differences in species composition. It is defined as the total number of species that are unique to either of the two sites standardized by the number of species in both sites.

\[ BC_\delta = \frac{S_i + S_j - 2S_{ij}}{S_i+S_j} \]

where \(S_x\) is the species count and \(S_{ij}\) is the sum of minimum abundances. Lately, this has seen considerable use within individual-based landscape genetic studies. Missing genotypes are set to average allele frequencies, that is to say that every missing genotype is considered to have all the alleles present in the entire population, but with probability equal to their global frequencies. Essentially, this removes the problem like in the situation and does so by taking the non-missing genotype’s genetic distance from the global genetic centroid (it’s cosmic man!). Here is the estimation using two loci.

D <- dist_bray(loci)
rownames(D) <- colnames(D) <- as.character(loci)
D

##       A:A   A:B   A:C   A:D   B:B   B:C   B:D   C:C   C:D   D:D
## A:A 0.000 0.000 0.125 0.125 0.125 0.000 0.000 0.000 0.000 0.000
## A:B 0.000 0.000 0.000 0.000 0.125 0.000 0.000 0.125 0.000 0.125
## A:C 0.125 0.000 0.000 0.125 0.125 0.125 0.125 0.125 0.000 0.000
## A:D 0.125 0.000 0.125 0.000 0.125 0.000 0.125 0.000 0.125 0.125
## B:B 0.125 0.125 0.125 0.125 0.000 0.000 0.000 0.125 0.000 0.125
## B:C 0.000 0.000 0.125 0.000 0.000 0.000 0.125 0.125 0.000 0.000
## B:D 0.000 0.000 0.125 0.125 0.000 0.125 0.000 0.125 0.125 0.125
## C:C 0.000 0.125 0.125 0.000 0.125 0.125 0.125 0.000 0.000 0.125
## C:D 0.000 0.000 0.000 0.125 0.000 0.000 0.125 0.000 0.000 0.125
## D:D 0.000 0.125 0.000 0.125 0.125 0.000 0.125 0.125 0.125 0.000

Euclidean Distance

Euclidean distance is the most straight-forward distance metric available as it is essentially straight-line distance based upon the allele frequencies in each population. It is given by:

\[ d_{eucl} = \sqrt{ \sum_{j=1}^L(p_{ij} - p_{kj})^2 } \]

where \(p_{ij}\) and \(p_{kj}\) are the frequencies of the \(j^{th}\) allele in both the \(i^{th}\) and \(j^{th}\) population. In this and the following distance examples, I am going to take the resulting distance matrix among all pairs of populations and put them into a Neighbor joining tree (via the nj() function from the ape package) as it may be easier to see differences in topologies rather than matrices.

It is perhaps easiest to think of Euclidean distance in x,y coordinate space. This distance can be estimated by stratum.distance() using the optional parameter method=‘eucl’ and it will return a dist matrix.

dist_euclidean(data)

##        101    102     32
## 101 0.0000 0.3335 0.3438
## 102 0.3335 0.0000 0.0874
## 32  0.3438 0.0874 0.0000

Cavalli-Sforza Distance

Another distance approach that is commonly used for microsatellite loci is Cavalli-Sforza distance, \(D_C\) (Cavalli-Sforza and Edwards, 1967). Here population allele frequencies are plot on the surface of a sphere (radius=1) using the square root of the allele frequencies.

\[ D_C = \frac{2}{\pi}\sqrt{(2-2cos\theta)} \]

The genetic distance, \(D_C\) is measured as the chord distance as indicated in Figure. The resulting Neighbor joining tree from this distance is shown in Figure

\[ D_{CS,m,n} = \frac{2}{\pi}\sqrt{ 2 * 1-\sum_{i=1}^{\ell} p_{m,i}*p_{n,i}} \]

dist_cavalli(data)

##        101    102     32
## 101 0.0000 0.2725 0.3368
## 102 0.2725 0.0000 0.2395
## 32  0.3368 0.2395 0.0000

The \(D_{PS}\) Distance

The Bray Curtis distance can also be estimated on groups of individuals. However, when it is done it is often represented as \(D_{PS} = (1 - P_S)\) (where \(P_S\) is the proportion of shared alleles), which I will follow here for simplicity such that the individual distances and the strata distances are not confused. The \(D_{PS}\) distance metric is directly related to Jaccard’s distance as:

\[ D_{PS} = \frac{-J_{\delta(A,B)}}{J_{\delta(A,B)}-2} \]

dist_bray(data)

##        101    102     32
## 101 0.0000 0.1823 0.1798
## 102 0.1823 0.0000 0.2303
## 32  0.1798 0.2303 0.0000

It should be noted here that the function dist_bray() is used for both individual distance AND strata distance depending upon what you pass to it. An individual distance is found by passing it a vector of loci whereas a stratum distance is returned by passing it a data.frame object (that has a \(stratum\) column). In the genetic_distance() function these are also differentiated by mode=“Bray” for the individual one and mode=“Dps” for the stratum.

Jaccard Distance

Jaccard distance is a set-theoretic distance quantifying dissimilarity. Assuming that loci are sets of alleles, the Jaccard dissimilarity between genotypes \(A\) and \(B\) is given by:

\[ J_{\delta(A,B)} = \frac{|A \bigcup B| - |A \bigcap B|}{|A \bigcup B|} \]

dist_jaccard(data)

##        101    102     32
## 101 0.0000 0.3084 0.3048
## 102 0.3084 0.0000 0.3743
## 32  0.3048 0.3743 0.0000

The reverse relationship between \(D_{PS}\) and \(J_{\delta(A,B)}\) is given as:

\[ J_{\delta(A,B)} = \frac{2+D_{PS}}{1 + D_{PS}} \]

Conditional Genetic Distance (cGD; Graph Distance)

Conditional genetic distance, cGD, is a measure based upon conditional genetic covariance and is distinctly different than these other measures as it is not a pair-wise distance metric. Rather it is the distance through a Population Graph topology whose construction is determined by the totality of the data. To estimate cGD from your data, you can use the genetic.distance() function as before and it will do right thing and return you a matrix. However, you should probably look into what a Population Graph really is prior to using it. You will find more information below as well as in the documents for the popgraph package itself. For consistency, the function is shown below BUT this data set with 3 populations is too small to be of interest for a network analysis (how many ways can you connect 3 things…)

dist_cgd(data)

##       101    102     32
## 101 0.000 1.0737 1.2657
## 102 1.074 0.0000 0.4985
## 32  1.266 0.4985 0.0000

Nei’s Genetic Distance

A very common metric of genetic distance, if you think the data you have are due to drift/mutation balance, is that of Nei. The implementation in gstudio of Nei’s distance is based upon the sample size correction from 1978, and is calculated as:

\[ I = \frac{(2N-1)\sum_{l=1}^L\sum_{m=1}^M p_{Alm}p_{Blm}}{\sqrt{ \sum_{l=1}^L(2N\sum_{m=1}^{M_A}p_{Alm}^2-1)(2N\sum_{m=1}^{M_B}p_{Blm}^2 -1)}} \]

for the “Genetic Identity” (where \(p_A\) is the allele frequencies at one population and \(p_B\) are the corresponding frequencies for the other, \(L\) is across loci and \(M\) is across alleles at the \(l^{th}\) locus).

Nei’s (1978) genetic distance, \(D_N\), is:

\[ D_N = -ln(I) \]

dist_nei(data)

##         101      102       32
## 101 0.00000 0.037411 0.049440
## 102 0.03741 0.000000 0.002092
## 32  0.04944 0.002092 0.000000

Comparing strata distance metrics

There are several more genetic distance metrics available and each may have its own set of assumptions. However, it is also true that across these metrics, there is some similarity between them. Just for illustrative purposes, we’ll look at the various strata distances and plot them against each other to look at the correlative structure between alternative measures using the \(ATPS\) locus and the entire beetle data set.

x <- arapat[, c(3, 13, 14)]
summary(x)

##    Population       ATPS          MP20    
##  32     : 19   05:05  :155   05:07  : 64  
##  75     : 11   03:03  : 69   07:07  : 53  
##  Const  : 11   09:09  : 66   18:18  : 52  
##  12     : 10   02:02  : 30   05:05  : 48  
##  153    : 10   07:09  : 14   05:06  : 22  
##  157    : 10   08:08  :  9   (Other):119  
##  (Other):292   (Other): 20   NA's   :  5

Now, we’ll grab the distance metrics

dist.euc <- genetic_distance(x, mode = "Euclidean")
dist.cgd <- genetic_distance(x, mode = "cGD")
dist.nei <- genetic_distance(x, mode = "Nei")
dist.dps <- genetic_distance(x, mode = "Dps")
dist.jac <- genetic_distance(x, mode = "Jaccard")

and then take the upper triangle of each and put them into a \(data.frame\).

df <- data.frame(Euclidean = dist.euc[upper.tri(dist.euc)])
df$cGD <- dist.cgd[upper.tri(dist.cgd)]
df$Nei <- dist.nei[upper.tri(dist.nei)]
df$Dps <- dist.dps[upper.tri(dist.dps)]
df$Jaccard <- dist.jac[upper.tri(dist.jac)]

Before we plot them, there is a bit of cleaning up to do in these data. For populations that have no alleles in common, Nei’s genetic distance will be \(Inf\) (e.g., \(-log(0)\)). Also, with cGD, sets of populations that are independent will also have a infinite distance (e.g., they are not connected so it is impossible to go through the graph from one population to the other). So with these data, we should first remove them.

and then we can plot them against each other and look at their correlations.

df <- df[is.finite(df$Nei), ]
df <- df[is.finite(df$cGD), ]
require(GGally)
ggpairs(df)

As you can see, there is a great deal of correlation between these parameters.

Nei’s \(G_{ST}\) Parameter

Gst(arapat$LTRS, arapat$Population, nperm = 99)

##      Gst     Hs     Ht P
## 1 0.5356 0.2319 0.4992 0

If the loci that you pass is a bunch of loci in a data.frame, it will return the single locus estimate as well as the multilocus estimate (based upon summing the heterosexuality and then estimating it as in Nei, not in just averaging the \(G_{ST}\) values as in Berg & Hamrick (XXXX), see XXXX for more on the differences).

Gst(arapat[, c(3, 7:14)], nperm = 99)

##        Locus    Gst     Hs     Ht  P
## 1       LTRS 0.5356 0.2319 0.4992  0
## 2        WNT 0.4364 0.3682 0.6534  0
## 3         EN 0.4283 0.2569 0.4493  0
## 4         EF 0.5550 0.1942 0.4364  0
## 5        ZMP 0.4724 0.1792 0.3397  0
## 6        AML 0.2954 0.5850 0.8302  0
## 7       ATPS 0.7103 0.2085 0.7197  0
## 8       MP20 0.3088 0.5663 0.8194  0
## 9 Multilocus 0.4544 2.5903 4.7474 NA

Hendrick’s \(G_{ST}^\prime\) Parameter

A correction to Nei’s \(G_{ST}\) was suggested for loci with a lot of alleles. This is because the maximum value for expected heterozygosity is determined by the number of alleles and as such \(G_{ST}\) for high allelic loci is not bound on the interval \([0,1]\) but is maxed out below 1.0. As a consequence, Hedrick

sort(unique(matrix(alleles(arapat$MP20), ncol = 1)))

##  [1] "01" "02" "03" "04" "05" "06" "07" "08" "09" "10" "11" "12" "13" "14"
## [15] "15" "16" "17" "18" "19"

Gst_prime(arapat$MP20, arapat$Population, nperm = 99)

##      Gst     Hs     Ht P
## 1 0.7228 0.5663 0.8194 0

In a similar fashion, the multilocus analog can be found by passing a data.frame to the function (again I am skipping the stratum variable to this function as the data has the strata in a column named ‘Population’).

Gst_prime(arapat[, c(3, 7:14)], nperm = 99)

##        Locus    Gst     Hs     Ht  P
## 1       LTRS 0.7015 0.2319 0.4992  0
## 2        WNT 0.6975 0.3682 0.6534  0
## 3         EN 0.5803 0.2569 0.4493  0
## 4         EF 0.6923 0.1942 0.4364  0
## 5        ZMP 0.5783 0.1792 0.3397  0
## 6        AML 0.7227 0.5850 0.8302  0
## 7       ATPS 0.9023 0.2085 0.7197  0
## 8       MP20 0.7228 0.5663 0.8194  0
## 9 Multilocus 0.6878 0.2674 0.5402 NA

Joost’s \(D_{est}\) Parameter

Following a discussion back-and-forth between Hedrick & Joost, Joost proposed an alternative measure \(D_{EST}\). The estimation of this parameter is found in a similar way as the other structure parameters.

Dest(arapat$MP20, arapat$Population, nperm = 99)

##     Dest     Hs     Ht P
## 1 0.5952 0.5663 0.8194 0

Dest(arapat[, c(3, 7:14)], nperm = 99)

##        Locus   Dest     Hs     Ht  P
## 1       LTRS 0.3496 0.2319 0.4992  0
## 2        WNT 0.4563 0.3682 0.6534  0
## 3         EN 0.2659 0.2569 0.4493  0
## 4         EF 0.3019 0.1942 0.4364  0
## 5        ZMP 0.1999 0.1792 0.3397  0
## 6        AML 0.6037 0.5850 0.8302  0
## 7       ATPS 0.6340 0.2085 0.7197  0
## 8       MP20 0.5952 0.5663 0.8194  0
## 9 Multilocus 0.3629 0.3238 0.5934 NA

Similarities in Parameters

As in genetic distance metrics, there is some similarity in output from these structure parameters. Here is a paired plot of the three parameters as above.

gst <- Gst(arapat)$Gst
gstp <- Gst_prime(arapat)$Gst
dest <- Dest(arapat)$Dest
df <- data.frame(Gst = gst, Gst_prime = gstp, Dest = dest)
require(GGally)
ggpairs(df)

You can tell from these plots which locus has a lot of alleles and which ones do not…

Making Populations

But first, lets create the populations.

freqs <- data.frame(Stratum = rep(c("Pop-1", "Pop-2"), each = 2), Locus = rep("TPI", 
    4), Allele = rep(c("A", "B"), time = 2), Frequency = c(0.75, 0.25, 0.25, 
    0.75))
freqs

##   Stratum Locus Allele Frequency
## 1   Pop-1   TPI      A      0.75
## 2   Pop-1   TPI      B      0.25
## 3   Pop-2   TPI      A      0.25
## 4   Pop-2   TPI      B      0.75

pops <- make_population(freqs, N = 100)
summary(pops)

##  Population        ID          TPI    
##  Pop-1:100   Min.   :  1.0   A:A :64  
##  Pop-2:100   1st Qu.: 25.8   A:B :69  
##              Median : 50.5   B:B :67  
##              Mean   : 50.5            
##              3rd Qu.: 75.2            
##              Max.   :100.0

You can check the creation of the individuals by looking at the resulting allele frequencies. However, it should be noted that you cannot often exactly mimic the allele frequencies due to randomness in selecting alleles (e.g., a fair coin does not give 50/50 with only a few tosses). The following shows they are pretty close.

frequencies(pops, stratum = "Population")

##   Stratum Locus Allele Frequency
## 1   Pop-1   TPI      A     0.770
## 2   Pop-1   TPI      B     0.230
## 3   Pop-2   TPI      A     0.215
## 4   Pop-2   TPI      B     0.785

Isolate Breaking: A Cautionary Tale

As a side note, this is a great place to point out the effects of isolate breaking (e.g., the issue of combining populations when we should not). Here each population is pretty close to no inbreeding (e.g., close to HWE).

Fis(pops[pops$Population == "Pop-1", ])

##   Locus      Fis
## 1   TPI -0.01637

Fis(pops[pops$Population == "Pop-2", ])

##   Locus     Fis
## 1   TPI 0.02237

But if we combine them we see a totally different thing.

Fis(pops)

##   Locus    Fis
## 1   TPI 0.3098

This is a rather poorly appreciated aspect highlighting the case where we do not quite understand the underlying biology yet still forge ahead blindly throwing the data at an analysis.

Mating Entire Populations

The mate() function takes two individuals and produces N offspring. However, we want to take random sets of individuals and make a new population that is the same size. This will take a few lines of code but can easily be accomplished. I’ll wrap it into a function with some comments to illustrate.

mate_population <- function(pop) {
    N <- dim(pop)[1]  # how many to make
    ret <- pop  # Replace loci in place
    for (i in 1:N) {
        # grab 2 random parent indices (no selfing)
        parents <- sample(1:N, size = 2, replace = FALSE)
        # mate them to make 1 offspring
        off <- mate(pop[parents[1], ], pop[parents[2], ], N = 1)
        # add the genotype back to the population.
        ret$TPI[i] <- off$TPI[1]
    }
    return(ret)
}

We can check to see if this works by iterating across generations and looking at allele frequencies. These should bounce around due to drift but should not be all over the place if N is reasonable. Here is a check on that.

T <- 50
df <- data.frame(Time = 1:T, Freq_A = rep(NA, T))
test_freq <- data.frame(Locus = "TPI", Allele = c("A", "B"), Frequency = c(0.5, 
    0.5))
pop <- make_population(test_freq, N = 100)
for (i in 1:T) {
    f <- frequencies(pop)
    df$Freq_A[i] <- f$Frequency[1]
    pop <- mate_population(pop)
}
ggplot(df, aes(x = Time, y = Freq_A)) + geom_line(size = 1.5, color = "darkred") + 
    theme_bw() + ylim(c(0, 1))

Migrating Individuals

This migration model will be rather simple. We will simply replace the genotypes of the first 10 percent of the adults (\(m=0.10\) which is rather large). Since each generation is made randomly, the position in the data.frame is somewhat arbitrary. Moreover, this function is useful only for the 1-locus data set we have created but simple examples can be simplistic and still provide some information.

migrate <- function(pop) {
    t <- pop$TPI[1:10]
    pop$TPI[1:10] <- pop$TPI[101:110]
    pop$TPI[101:110] <- t
    return(pop)
}

Plotting Structure in Drift/Migration Simulations

OK, lets put it together and plot the results.

df <- data.frame(Time = 1:50, Gst = rep(0, 50))
for (t in 1:50) {
    
    # find structure
    x <- Gst(pops, stratum = "Population")
    df$Gst[t] <- x$Gst[1]
    
    # migrate
    pops <- migrate(pops)
    
    # mate
    pops[1:100, ] <- mate_population(pops[1:100, ])
    pops[101:200, ] <- mate_population(pops[101:200, ])
    
}
ggplot(df, aes(x = Time, y = Gst)) + geom_line(size = 1.5, color = "darkred") + 
    theme_bw() + ylim(c(-0.1, 0.4))

So with just a few lines of code, we can create some simulations that show the change in structure as a function of both drift (small \(N\) in these populations) and migration (\(m=0.1\)). The utility of R is that it opens up a wide range of potential simulation work for you and with gstudio you can work with genetic data just like any other kind of data.