\(F\)-statistics with hierfstat and Hickory

What you see here is an R notebook. We’ll talk about how to use them in class on Tuesday. Let’s start with a simple example. We’ll use the data from Isotoma petraea that we discussed in class to illustrate hierfstat and Hickory.

Converting the Isotoma data to a different format

We’ll start by using hierfstat. It needs the data in a different format from what is available on the website. This piece of code reads in the data directly from the website and converts it to the format that hierfstat uses. Other than the cute trick of reading the CSV file directly from the server rather than downloading it first, you can ignore all of the other coding.

library("tidyverse")

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

## The following line clears out any objects in memory. I make it a habit to do this so that I can make sure
## that things left over from the last time I ran R don't confuse things.
##
rm(list = ls())

dat <- read_csv("http://darwin.eeb.uconn.edu/eeb348-resources/isotoma.csv",
                col_types = cols(pop = col_character(), .default = col_integer()))

## If you know R, you'll recognize that I could have expressed this much more succintly, but for anyone who
## cares to follow the logic, this will be easier to see.
##
for (i in 1:nrow(dat)) {
  for (j in 2:ncol(dat)) {
    if (!is.na(dat[i, j])) {
      if (dat[i, j] == 0) {
        dat[i, j] <- 11
      } else if (dat[i, j] == 1) {
        dat[i, j] <- 12
      } else if (dat[i, j] == 2) {
        dat[i, j] <- 22
      }
    }   
  }
}
## This line converts the population abbreviatons to numbers
##
dat$Location <- as.numeric(as.factor(dat$pop))
## The first line puts Location in the first column. The second drops the pop column.
##
dat <- relocate(dat, Location) %>%
  select(-pop)
dat <- as.data.frame(dat)

Using hierfstat

We use wc() to estimate Weir and Cockerham’s \(F\)-statistics.¹

library(hierfstat)

dat_fst <- wc(dat, diploid = TRUE)
dat_fst

$FST
[1] 0.03869543

$FIS
[1] 0.539859

That’s all there is to it. You don’t need to include the diploid = TRUE, but it’s good form to do so. If you happened to try this with the original file you downloaded, you’d get an error because the data aren’t in the right format. The numbers above match what I mentioned in lecture (with a small rounding error in the 4th decimal place for \(F_{IS}\))

Using `Hickory

Using Hickory is pretty straigthforward.² It is, however, a bit more challenging to install than a regular R package. It uses RStan, and you’ll need to refer to this guide to install RStan. In particular, you’ll need to “configure the C++ toolchain”. If that sounds frightening, don’t worry about it. You don’t need to learn to program in C++, and you don’t even have to know how to compile a C++ program. The links will walk you through installing a C++ compiler on your machine. Once you’ve done that, installing Rstan is straigtforward, and then you’re ready to install Hickory. Just do the following:

install.packages("devtools")
install.packages(c("bayesplot"))
devtools::install_github("kholsinger/Hickory", build_vignettes = TRUE)

That last step will take a while, maybe as much as 10 minutes, and you’ll see a bunch of arcane messages in your R console window. So long as you don’t see something that says “Error”, you’re fine. If you do, let me know, and we’ll figure out how to work around it.

Once you’ve got Hickory installed, running the analysis is pretty straightforward.

library("Hickory")

## This line allows Hickory to run four chains simultaneously (assuming that your computer has at least)
## four cores
##
options(mc.cores = parallel::detectCores())

genos <- read_marker_data("http://darwin.eeb.uconn.edu/eeb348-resources/isotoma.csv")
## The "refresh = 0" prevents a lot of messages from being printed in this document. You probably will want
## to leave it out when you run your analyses so that you can see the progress messages and be reassured
## that there's something going on.
##
fit <- analyze_codominant(genos, refresh = 0)

Inference for Stan model: analyze_codominant.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

          mean se_mean    sd     2.5%      25%      50%      75%
f        0.345   0.002 0.102    0.148    0.272    0.344    0.415
theta    0.077   0.002 0.047    0.018    0.043    0.066    0.099
lp__  -121.868   0.194 4.339 -131.204 -124.632 -121.460 -118.723
         97.5% n_eff  Rhat
f        0.549  3659 1.000
theta    0.194   708 1.003
lp__  -114.641   501 1.004

Samples were drawn using NUTS(diag_e) at Sat Sep  4 14:59:39 2021.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

If you compare those numbers with the numbers I reported in lecture, you’ll see that they are very close. They aren’t exactly the same because we’re simulating samples from the posterior distribution rather than calculating it exactly. The figure shows the median (the dot), the 50 percent credible interval (the thick dark blue line), and the 90 percent credible interval (the thin light blue line) for each of the parameters.

Lab exercise #2

We are going to use a subset of the data that Rachel Prunier and collaborators (including me) used to analyze the genetic structure of Protea repens https://doi.org/10.3732/ajb.1600232. You’ll find the data on the course website at http://darwin.eeb.uconn.edu/eeb348-resources/repens-outliers.csv and at http://darwin.eeb.uconn.edu/eeb348-resources/repens-outliers.stru.³

1. Downlad the data in the .stru file. Then work through the following steps to analyze the data with wc().⁴

   1. Convert the data to a form that hierfstat can use by running read.structure() from adegenet.⁵. Note: 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. If everything worked properly, you’ll see the following message:

    Converting data from a STRUCTURE .stru file to a genind object... 
    

   2. Now use wc() to produce estimates of \(F_{IS}\) and \(F_{ST}\).

   3. Now convert the genind object that you used with wc() to a hierfstat data frame (using genind2hierfstat()), and use that dataframe to construct 95 percent confidence intervals for \(F_{IS}\) and \(F_{ST}\)

    boot.vc(repens_df[, 1], repens_df[, -1], diploid = TRUE)$ci
    

   repens_df is the name of the data frame I used. The first argument selects only the pop column (the one with the populations), and the second excludes the pop column, leaving a data frame with only the genotype data.

2. Use Hickory and the data as provided in the .csv file to estimate \(F_{IS}\) and \(F_{ST}\). Report both the posterior mean and the 95 percent credible intervals. Note: The analysis took about 15 minutes on my MacBook. Be sure to plan ahead and leave yourself enough time to commplete the assignment.

3. Is there evidence of inbreeding within populations?

4. Is there evidence of genetic differentiation among populations?

5. Does your answer to #3 or #4 depend on whether you look at results from hierfstat() or Hickory.

6. BONUS QUESTION: Take a look at the help for analyze_codominant() and see if you can explain why the estimates from hierfstat and Hickory look a bit different even though we have quite a bit of data available.

Feel free simply to embed the R code in a copy of this notebook if you’re so inclined. but if you do, be sure to state explicitly how to interpret the output.[^I know how to do it, of course, but I want to make sure you know how to do it too.]

---

1. You can probably guess why the function is called wc().↩︎

2. Of course I’m probably biased. I wrote Hickory.↩︎

3. I’ll explain what the .stru extension refers to next week.↩︎

4. You can also use the trick of simply using the URL above to read the data directly from the course website.↩︎

5. You’ll need to use install.packages() to install adegenet↩︎