Allele frequency data

You must specify a matrix of allele frequency data for your samples. (Make sure the data are of class matrix, and that it’s not a data.frame.) I assume that the data consist of bi-allelic SNPs. At each locus, you pick an allele to count across all samples (it doesn’t matter whether it’s randomly chosen or whether it’s always the major or minor allele). The frequency of the counted allele at a locus in a sample is the number of times the counted allele is observed at a locus divided by the number of chromosomes genotyped in that sample. A sample can consist of a single individual or multiple individuals lumped together. So, a successfully genotyped diploid individual heterozygous at a particular locus would have an allele frequency of 0.5. If the sample is a population of 13 haploids, of which 12 have the counted allele at a given locus, the frequency in that sample at that locus would be 12/13.

The matrix of allele frequencies should have one row per sample and one column per locus. Missing data should be denoted with the value NA. An small example allele frequency data matrix is shown below:

An full example allele frequency data matrix is included in the conStruct.data object included with the package.

# load the example data object
data(conStruct.data)

# look at the allele frequency data 
#   for the first 5 populations and 10 loci
conStruct.data$allele.frequencies[1:5,1:10]
#>         loc1 loc2 loc3 loc4 loc5 loc6 loc7 loc8 loc9 loc10
#> sample1 0.00    0    0 0.00    0    0  0.0 0.00 0.80   0.0
#> sample2 0.00    0    0 0.00    0    0  0.0 0.00 0.20   0.0
#> sample3 0.00    0    0 0.05    0    0  0.0 0.05 0.05   0.0
#> sample4 0.00    0    0 0.05    0    0  0.0 0.00 0.05   0.0
#> sample7 0.65    0    0 0.00    0    0  0.1 0.00 0.20   0.8

Geographic sampling coordinates

You must specify a matrix of geographic sampling coordinates, which will be used for plotting the results of the analysis. This should be a matrix with two columns that give the sample x-coordinates (longitude) and y-coordinates (latitude), respectively. The order of rows of the matrix should be the same as the order of the rows of the allele frequencies matrix. If you specify longitude and latitude, they should be in decimal degrees.

A full example sampling coordinate data matrix is included in the conStruct.data object included with the package.

# load the example data object
data(conStruct.data)

# look at the geographic sampling coordinates 
#   for the first 5 populations
conStruct.data$coords[1:5,]
#>           Lon     Lat
#> [1,] -84.4328 42.7828
#> [2,] -84.4168 42.7828
#> [3,] -84.4008 42.7828
#> [4,] -84.3848 42.7828
#> [5,] -84.4328 42.7988

Geographic distance matrix

If you choose to run the spatial model implemented in conStruct, you must specify a matrix of pairwise geographic distance between all samples. If the coordinates of the samples are real locations on Earth (as opposed to simulated coordinates), I recommend calculating pairwise great-circle distance between sampling coordinates (using, e.g., geosphere::distm or geosphere::distGeo).

The order of the samples in the geographic distance matrix should match both that of the geographic coordinates and that of the allele frequency data matrix, and all three matrices should have the same number of rows.

The geographic distance matrix you specify should be the full matrix (that is, not the upper- or lower-triangles), with values of 0 on the diagonal entries.

A full example geographic distance matrix between all samples is in the conStruct.data object included with the package.

# load the example data object
data(conStruct.data)

# look at pariwise geographic distance 
#   between the first 5 populations
conStruct.data$geoDist[1:5,1:5]
#>           [,1]      [,2]      [,3]      [,4]     [,5]
#> [1,] 0.0000000 0.8122983 1.6245967 2.4368950 1.106773
#> [2,] 0.8122983 0.0000000 0.8122983 1.6245967 1.372809
#> [3,] 1.6245967 0.8122983 0.0000000 0.8122983 1.965599
#> [4,] 2.4368950 1.6245967 0.8122983 0.0000000 2.676167
#> [5,] 1.1067733 1.3728092 1.9655990 2.6761670 0.000000

STRUCTURE to conStruct

The program STRUCTURE is one of the most widely used methods for model-based clustering in population genetics. Many existing programs, including plink (v1.9 and above) and PgdSpider, convert data from diverse formats (including .vcf files) into STRUCTURE format. In this section of the vignette, I walk through an example of converting a STRUCTURE format data file into a conStruct format data file.

Example STRUCTURE format dataset, with one row per individual 
and two columns per locus. The first column gives sample names, the 
second refers to the sample locations, and the last 8 columns give 
genotype data for four loci. The numbers in the genotype data refer to 
the allele present at that locus: A1 = `1`, A2 = `2`, missing = `-9`.

conStruct.data <- structure2conStruct(infile = "~/Desktop/myStructureData.str",
                                      onerowperind = TRUE,
                                      start.loci = 3,
                                      start.samples = 1,
                                      missing.datum = -9,
                                      outfile = "~/Desktop/myConStructData")

Example STRUCTURE format dataset, with two rows per individual 
and one column per locus. The first column gives sample names, the 
second refers to the sample locations, and the last 8 columns give 
genotype data for 8 loci. The numbers in the genotype data refer to 
the allele present at that locus: A1 = `1`, A2 = `2`, missing = `0`.

conStruct.data <- structure2conStruct(infile = "~/Desktop/myStructureData.str",
                                      onerowperind = FALSE,
                                      start.loci = 3,
                                      start.samples = 1,
                                      missing.datum = 0,
                                      outfile = "~/Desktop/myConStructData")

pop.data.matrix <- matrix(NA,nrow=4,ncol=ncol(conStruct.data))
for(i in 1:nrow(pop.data.matrix)){
    pop.data.matrix[i,] <- colMeans(
                                conStruct.data[
                                    which(pop.index==i),,
                                    drop=FALSE
                                ],na.rm=TRUE
                            )
}

Microsatellites

This method is designed to run on large datasets consisting of bi-allelic SNPs. If you have a microsatellite dataset and you wish to run conStruct, the first consideration is whether you have sufficient data. You should have more loci than samples in your data matrix (i.e., your data matrix should have more columns than rows).

If that’s the case, the second consideration is how to format your microsat data so that you can run conStruct. There are two standard ways of “SNP-ifying” a microsat dataset.

The first is to lump all microsatellite alleles present at a locus into two categories: “major” and “other”. The “major” allele is the allele that occurs most frequently at a particular locus; all other alleles are put in the “other” bin. You then can create a dataset in which you only report the frequency of the major allele, effectively reducing the number of alleles per locus to 2. This method has the disadvantage of throwing out data, but acknowledges the simplex relationships between alleles at a locus (the sum of the frequencies of all alleles at a locus must be 1).

The second approach, introduced by Cavalli-Sforza, is to split out each allele at a locus into a separate pseudo-locus consisting of only that allele. That is, if you had 4 alleles present in the genotyped sample at a particular locus, at frequencies {0.4,0.3,0.1,0.2}, you would split those out into 4 separate columns in your data matrix (pseudo-loci), with frequencies in the sampled population of {0.4,0.3,0.1,0.2}. This approach has the advantage of not throwing data away, but does not acknowledge the inter-allele dependence structure in frequencies, and therefore introduces some pseudoreplication into the dataset. This pseudoreplication may make you overconfident in your results, as the credible intervals on parameter estimates may be artificially narrow.

I would recommend trying both approaches, and comparing the estimates of pairwise relatedness you get from each to those derived from the raw microsatellite data to see which best recovers the patterns of relatedness in the data. I also recommend running conStruct on datasets SNP-ified using both approaches, and comparing the results.