R Notebook

For Lab 3 we used data from Camak et al. (2021)¹ to explore population structure using LEA. This week we’re going to run a couple of additional analyses to explore these data further, and we’re going to use the results from these analyses to re-examine some of the conclusions Camak et al. reached. I’ll be asking you not only to produce numerical/graphical results from these analyses, but also to interpret what they mean biologically.

But first I have to introduce the methods, and we’ll use the Protea repens data you’ve already seen to do that.

PCA

First, we’ll run a principal components analysis of the data. To do this we first have to get the data into R using read.structure() from adegenet.²

library(adegenet)
repens <- read.structure("http://darwin.eeb.uconn.edu/eeb348-resources/repens-outliers.stru")

summary(repens)

// Number of individuals: 662
// Group sizes: 12 52 8 36 46 22 41 72 24 39 47 22 12 28 36 42 58 30 35
// Number of alleles per locus: 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
// Number of alleles per group: 247 334 233 325 321 303 316 326 299 295 326 305 292 308 319 328 332 319 323
// Percentage of missing data: 48.99 %
// Observed heterozygosity: 0.13 0.38 0.16 0.28 0.07 0.06 0.15 0.25 0.3 0.02 0.15 0.12 0.13 0.46 0.12 0.23 0.22 0.14 0.34 0.18 0.24 0.16 0.15 0.21 0.25 0.32 0.09 0.18 0.22 0.24 0.18 0.06 0.03 0.08 0.3 0.08 0.42 0.24 0.21 0.14 0.11 0.18 0.42 0.36 0.17 0.44 0.33 0.31 0.16 0.2 0.17 0.34 0.18 0.33 0.17 0.33 0.12 0.14 0.15 0.06 0.24 0.1 0.19 0.2 0.27 0.2 0.12 0.07 0.16 0.23 0.09 0.37 0.16 0.14 0.25 0.15 0.2 0.34 0.34 0.28 0.16 0.26 0.27 0.28 0.23 0.19 0.07 0.09 0.1 0.07 0.13 0.2 0.15 0.36 0.25 0.2 0.1 0.14 0.1 0.37 0.17 0.09 0.49 0.42 0.09 0.45 0.19 0.13 0.15 0.17 0.06 0.31 0.19 0.13 0.1 0.13 0.21 0.34 0.07 0.26 0.16 0.06 0.16 0.06 0.07 0.07 0.23 0.39 0.18 0.15 0.34 0.15 0.11 0.21 0.04 0.38 0.18 0.12 0.27 0.47 0.1 0.13 0.29 0.23 0.24 0.07 0.36 0.1 0.08 0.05 0.2 0.45 0.21 0.25 0.26 0.17 0.22 0.15 0.27 0.12 0.17 0.26 0.22 0.14 0.19 0.2 0.29 0.28 0.18 0.15 0.32 0.04 0.09
// Expected heterozygosity: 0.5 0.46 0.32 0.36 0.2 0.27 0.37 0.34 0.43 0.29 0.3 0.19 0.24 0.48 0.42 0.32 0.33 0.44 0.41 0.47 0.33 0.24 0.47 0.3 0.5 0.32 0.16 0.3 0.38 0.48 0.45 0.22 0.19 0.43 0.35 0.31 0.49 0.46 0.45 0.25 0.16 0.29 0.5 0.49 0.22 0.5 0.5 0.42 0.27 0.49 0.26 0.5 0.35 0.49 0.32 0.46 0.15 0.35 0.24 0.19 0.4 0.5 0.24 0.36 0.5 0.34 0.23 0.45 0.38 0.3 0.16 0.49 0.42 0.18 0.32 0.27 0.48 0.49 0.49 0.47 0.33 0.48 0.35 0.47 0.31 0.44 0.49 0.32 0.28 0.3 0.16 0.28 0.19 0.41 0.49 0.26 0.14 0.32 0.37 0.49 0.48 0.16 0.45 0.45 0.34 0.49 0.24 0.47 0.48 0.3 0.17 0.46 0.36 0.18 0.18 0.18 0.36 0.5 0.37 0.32 0.41 0.15 0.23 0.11 0.19 0.18 0.36 0.42 0.31 0.19 0.41 0.47 0.37 0.25 0.3 0.45 0.49 0.41 0.45 0.5 0.28 0.38 0.46 0.32 0.5 0.49 0.49 0.13 0.19 0.42 0.25 0.48 0.49 0.28 0.41 0.26 0.36 0.19 0.34 0.22 0.39 0.47 0.28 0.18 0.43 0.3 0.3 0.49 0.41 0.22 0.43 0.12 0.34

If we summarize the data, we’ll see that nearly fifty percent of the data is missing, i.e., fifty percent of the individual/locus combinations are empty because we couldn’t score that individual at that locus.³ This amount of missing data isn’t uncommon with genotype-by-sequencing (GBS) data when you don’t have a reference genome, but it is a problem for analysis by principal components. We fill the empty cells by replacing them with the mean allele frequency at the locus in question. This gives us complete data, but it understates the degree to which groupings in the data are different from one another. We also need a data frame that contains the population labels. Here’s how we do those two things:⁴

repens_mod <- tab(repens, freq = TRUE, NA.method = "mean")
repens_df <- data.frame(Population = pop(repens))

Now we run the principal components analysis (PCA) and plot the results. IMPORTANT NOTE: If you leave out the “u” in “colour” you won’t get the results that you want.

library(ggfortify)

repens_pca <- prcomp(repens_mod)
## must be spelled "colour" with a "u"
##
autoplot(repens_pca, data = repens_df, colour = "Population",
         frame.type = "norm") + 
  theme_bw()

With these data the first two principal components don’t account for much of the variation, only a little more than 15%. None of the populations are cleanly distinguished from others, but there are a few that don’t overlap. If you find the ellipses confusing, leave off the frame.type argument.

autoplot(repens_pca, data = repens_df, colour = "Population") + 
  theme_bw()

DAPC

DAPC is an acronym for “discriminant analysis of principal components”. As the name suggests, it combines principal components analysis (like what we’ve just seen) with another statistical technique, discriminant analysis. Principal components analysis is used to reduce the number of dimensions in a data set without trying to tell groups apart. Discriminant analysis takes pre-defined groups and tries to find mathematical functions that can be used to “discriminate” among them.

The first step in a DAPC analysis is to find genetic clusters in the data. The find.clusters() function in adegenet uses a technique known as k-means clustering.⁵ When you run find.clusters(), you’ll be asked how many PCs to retain. You’ll want to run this interactively on the console rather than through an R notebook so that you can answer the questions and see the plots before you answer

clusters <- find.clusters(repens)

In this case I picked 75 for the number of PCs to retain, because that gives me about 80 percent of the total variance. I picked 12 for the number of clusters, since that is the smallest BIC. Since 10-15 are all pretty close, if I were doing this for real, I’d probably look at the results for all of those possibilities to see if one made more sense than another.⁶

Now you’re ready to run dapc(). Specify the same number of PCs that you did with find.clusters(). Again you’ll want to run this interactively, because dapc() will ask you how many discriminant functions to retain. There isn’t a clear cutoff with these data, but there’s a reasonably large jump from 3 to 4, and I picked 3.

repens_dapc <- dapc(repens, pop = pop(repens), n.pca = 75)

Protea repens DAPC

Once you’ve got that, then it’s just a matter of plotting the results.

scatter(repens_dapc)

The numbers in squares are the population labels (numbers from 1-19 in this case), the points correspond to individuals with the color matching the color of the population to which they belong, and the ellipses show the boundaries of where you expect most individuals belonging to that population to fall.

By default, scatter() plots the first two discriminant axes, but I specified three discriminant functions, so we can also plot axis 1 against axis 3 and axis 2 against axis 3 to get a more complete picture.

scatter(repens_dapc, 1, 3)

scatter(repens_dapc, 2, 3)

As you can see, population 8 seems to separate a bit from the remaining populations in these plots.

Notice that if we compare the number of individuals belonging to each population before dapc(), i.e., the number collected in a particular location, it differs from the number belonging to each population after dapc(). This means that some individuals are more similar to those from other individuals from other collection sites than to individuals in the site where they were collected.

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────── tidyverse 1.3.1 ──
✓ tibble  3.1.4     ✓ dplyr   1.0.7
✓ tidyr   1.1.3     ✓ stringr 1.4.0
✓ readr   2.0.1     ✓ forcats 0.5.1
✓ purrr   0.3.4     
── Conflicts ──────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

results <- data.frame(Initial = summary(repens_dapc)$prior.grp.size,
                      Final = summary(repens_dapc)$post.grp.size) %>%
  rename(Population = Initial.Var1,
         Initial = Initial.Freq,
         Final = Final.Freq) %>%
  select(-Final.Var1)
results

If we’re interested, we can determine which individuals have been reclassified.

results <- data.frame(Initial = repens_dapc$grp,
                      Final = repens_dapc$assign) %>%
  filter(Initial != Final)
results

You may or may not need to see which individuals have been classified, but at least you have the code to do so if you think it will be useful.

Project 1

When you submit your project please remember to use this naming convention:

Project1_first name_last name.extension

For example,

Project1_Kent_Holsinger.nb.html

1. Use the data from http://darwin.eeb.uconn.edu/eeb348-resources/gila-trout.stru to perform a principal components analysis and a DAPC of genetic variation in gila trout from southwestern New Mexico. To read the data from this file you’ll need to know that rhere are 154 individuals, the genotype ID is in column 3, and the population ID is in column 1. There is one extra column : column 2. There is not a row of allele names. Individuals appear on two lines rather than one. Be sure to say how many PCs you retained for the DAPC analysis and why you picked that number. Also be sure to say how many clusters you picked and how many discriminant functions you used.

2. Compare your results from this week’s PCA and DAPC with last week’s LEA analysis. Do you see noticeable differences among the patterns of genetic structure depending on the type of analysis you perform?

3. Camak et al. assert that “All lineages have experienced population bottlenecks associated with mortality from drought and severe wildfires.” The data you have don’t allow you to detect any effects of drought or wildfires, but are the patterns of genetic structure you detected consistent with population bottlenecks.⁷

4. Given that we expect populations that are closer to one another to be more genetically similar, is the geographical distribution of the populations illustrated in Figure 1 of Camak et al. consistent with the results of your analyses?

---

1. Camak, D.T., M.J. Osborne, and T.F. Turner. 2021. Population genomics and conservation of Gila trout (Onchorhyncus gilae). Conservation Genetics https://doi.org/10.1007/s10592-021-01355-0.↩︎

2. Remember that there are 662 genotypes (individuals), 173 markers (loci), column 1 contains the labels for genotypes (individuals), column 2 contains the population factor (the population/locality from which a particular sample was collected), there are no optional columns, row 1 contains the marker (locus) names, and genotypes are coded in two rows.↩︎

3. We see this in the line labeled “Percentage of missing data”.↩︎

4. Note: If you happen to have latitudes and longitudes for your collection sites, you could also add those to the data frame. Then you could plot individual points on a map colored by scores on one of the principal component axes.↩︎

5. Given K, this approach finds the center of K clusters such that each individual belongs to one cluster and the distance between every individual and the center of the cluster to which they belong is minimized. See this Wikipedia article for more details.↩︎

6. Depending on your data, the BIC plot may look more like the cross-entropy plot you saw with the gila trout data and LEA last week. In that case, pick a number of clusters that gives you close to the smallest BIC.↩︎

7. This is stretching you a bit, since we won’t talk about genetic drift for a couple of weeks, but you should be able to answer based on your basic understanding of biology.↩︎