---
title: R-code for 'Heterozygosity at neutral and immune loci does not influence neonatal
  mortality due to microbial infection in Antarctic fur seals'
author: Vivienne Litzke, Meinolf Ottensmann, Jaume Forcada, Louise Heitzmann & Joseph
  I. Hoffman
output:
  html_document:
    df_print: paged
header-includes:
- \usepackage{pdflscape}
- \newcommand{\blandscape}{\begin{landscape}}
- \newcommand{\elandscape}{\end{landscape}}
---

```{r setup, include=FALSE, message=FALSE, warning=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

This document provides all the `R code` used in our paper. Both the Rmarkdown file and the data can be downloaded from the accompanying GitHub repository on (https://github.com/vlitzke/HeterozygosityPupSurvival) as a zip archive containing all the files. We recommend to download or clone this [GitHub repository](https://github.com/vlitzke/HeterozygosityPupSurvival) in order to access the documentation together with all the files that are needed to repeat analyses shown in this document. Just click on the link above and then on the green box `Clone or download`. In order to function properly, the same structure of folders must be kept. If you have any questions, don't hesitate to contact vivienne.litzke[at]mail.mcgill.ca

The data originates from samples collected from a colony of Antarctic fur seals (*Arctocephalus gazella*) at Bird Island, South Georgia between the years of 2000 and 2014. We investigated the effects of neutral and immune gene heterozygosity on early mortality due to bacterial infection using the `inbreedR` package.[^1]

[^1]: Stoffel, M. A., Esser, M., Kardos, M., Humble, E., Nichols, H., David, P., & Hoffman, J. I. (2016). inbreedR: an R package for the analysis of inbreeding based on genetic markers. Methods in Ecology and Evolution, 7(11), 1331-1339.

- If you have downloaded the project from github then you will see that:
- The raw data required are located in the folder `data/`
- The *Arctocephalus gazella* transcriptome[^2] may be downloaded [here](http://ramadda.nerc-bas.ac.uk/repository/entry/show/Polar+Data+Centre/NERC-BAS+Datasets/Genomics/Transcriptomes/Arctocephalus_gazella/Arctocephalus_gazella_transcripts.fasta?entryid=synth%3A2d2268fe-907c-45b0-a493-0a6cab8642e6%3AL1RyYW5zY3JpcHRvbWVzL0FyY3RvY2VwaGFsdXNfZ2F6ZWxsYS9BcmN0b2NlcGhhbHVzX2dhemVsbGFfdHJhbnNjcmlwdHMuZmFzdGE%3D) and saved as `arc_gaz_transcriptome.fasta` in `data`.
- This pipeline invokes the `MIcroSAtellite identification tool` for primer identification. Click on the following link for details and how to install it: [MISA](http://pgrc.ipk-gatersleben.de/misa/misa.html)[^3]. 
- Primer development was conducted using [primer3](http://primer3.sourceforge.net/releases.php)[^4].
- Additionally, the `R packages` listed below are required and may be installed on your system.

[^2]: Humble, E., Thorne, M.A., Forcada, J. & Hoffman, J.I., (2016). Transcriptomic SNP discovery for custom genotyping arrays: impacts of sequence data, SNP calling method and genotyping technology on the probability of validation success. BMC research notes, 9(1), p.418.
[^3]: Thiel, T., (2003). MISA—Microsatellite identification tool. Website http://pgrc. ipk-gatersleben.
[^4]: Untergasser, A., Nijveen, H., Rao, X., Bisseling, T., Geurts, R. & Leunissen, J.A., (2007). Primer3Plus, an enhanced web interface to Primer3. Nucleic acids research, 35(suppl_2), pp.W71-W74.

***

### Download packages and libraries

In order to repeat analyses presented in this manuscript a number of packages that extend the functionalities of base `R` are required. These can be installed using the code shown below.

```{r install, comment="install packages", eval=FALSE}
install.packages('inbreedR')
install.packages("readxl")
install.packages("ggplot2")
install.packages("gridExtra")
install.packages("stringi")
install.packages("adegenet")
install.packages("AICcmodavg")
install.packages("reshape2")
install.packages("kableExtra")
source("https://bioconductor.org/biocLite.R")
biocLite("qvalue")
```

```{r, comment="Load packages", message=FALSE}
library(inbreedR)
library(readxl)
library(magrittr)
library(ggplot2)
library(grid)
library(gridExtra)
library(AICcmodavg)
library(Matrix)
library(lme4)
library(qvalue)
library(adegenet)
library(reshape2)
library(kableExtra)
```

In order to use `inbreedR`, the working format is typically an *individual x loci* matrix, where rows represent individuals and every two columns represent a single locus. If an individual is heterozygous at a given locus, it is coded as 1, whereas a homozygote is coded as 0, and missing data are coded as NA. 

The first step is to read the data from an excel file. Our original table includes, plate number, well number, species, id, year, health status (represented by a binomial with 0 for healthy and 1 for infected), birth weight, and the following markers (a and b for alleles). 

```{r read table, message=FALSE, warning=FALSE}
## read data
seals <-  readxl::read_xlsx("data/data.xlsx")[1:234, ]
## express alleles as numerals
seals[8:ncol(seals)] <- lapply(seals[8:ncol(seals)], as.numeric) 
```

Here's an example of what the data frame looks like:
```{r seal header}
head(seals[1:6,4:11])
```

Since demographic data is present in the beginning of our data frame, we will start our new genotype file from the 8th column onwards. The function `convert_raw` converts a common format for genetic markers (two columns per locus) into the `inbreedR` working format. Afterwards, `check_data` allows us to test whether the genotype data frame has the correct format for subsequent analyses that use `inbreedR` functions.

```{r convert, results='hide', error=FALSE, message=FALSE, comment=FALSE}
seals_geno <- convert_raw(seals[8:ncol(seals)])
invisible(lapply(seals_geno, table, useNA = "always"))
check_data(seals_geno, num_ind = 234, num_loci = 61)
```

***

\pagebreak

### Separate markers
Divide the neutral and immune markers from their respective columns in the adjusted inbreedR format, and compute standard multilocus heterozygosity (sMLH).[^5] 

[^5]: Coltman, D. W. and J. Slate. 2003. Microsatellite measures of inbreeding: a meta-analysis. Evolution 57:971–983.
```{r}
## subset markers based on type
imm_msats <- seals_geno[, 1:13]
neut_msats <- seals_geno[, 14:61]

## estimate sMLH
all_het <- sMLH(seals_geno)
neutral_het <- sMLH(neut_msats)
immune_het <- sMLH(imm_msats)
```

### Create and reshape dataframes

Take out id, health, marker types, and birth weight as variables. Tack on the year. 
```{r message=FALSE, warning=FALSE}
sealdf <- data.frame(id = seals[[4]],  health = factor(seals[[6]]), 
                       All = all_het, Neutral = neutral_het, Immune = immune_het)

sealdf_reshaped <- reshape2::melt(sealdf)

birthweight <- as.numeric(as.character(seals[["Birthweight"]]))

sealdf_weight <- data.frame(id = seals[[4]],  health = factor(seals[[6]]), birthweight, 
                             All = all_het, Neutral = neutral_het, Immune = immune_het)

sealdf_year <- cbind(sealdf_weight, year = seals[[5]]) 
```

***

### Calculate *g*~2~
*g*~2~ is a proxy for identity disequilibrium. It is a measure of two-locus disequilibrium, which quantifies the extent to which heterozygosities are correlated across pairs of loci.[^6] This allows us to take a look at our neutral marker heterozygosity to determine if there is variation in inbreeding in the population. 

[^6]:David, P., Pujol, B., Viard, F., Castella, V., & Goudet, J. (2007). Reliable selfing rate estimates from imperfect population genetic data. Molecular ecology, 16(12), 2474-2487.

```{r neutral g2, echo= FALSE, message=FALSE, warning=FALSE, results='hide'}
## calculate g2
g2_neut <- g2_microsats(neut_msats, nperm = 9999, nboot = 9999)

## format output
g2_neut_bs <- data.frame(bs = g2_neut$g2_boot,
                            lcl = g2_neut$CI_boot[[1]],
                            ucl = g2_neut$CI_boot[[2]],
                            g2  = g2_neut$g2,
                            p = g2_neut$p_val)
```

Plot the distribution of g2 estimates:
```{r neutral g2 plot, warning=FALSE}
g2_neut_bs_hist <-
  ggplot2::ggplot() +
  theme_classic() +
  geom_histogram(binwidth = 0.00025, data = g2_neut_bs, aes(x = bs),
                 color = "#0294A5",
                 fill = "#0294A5") +
  geom_errorbarh(data = g2_neut_bs,
                 aes(y = 1600, x = g2, xmin = lcl, xmax = ucl),
                 color = "black", size = 0.7, linetype = "solid") +
  geom_linerange(data = g2_neut_bs,
                 aes(ymin = 0, ymax = 1600, x = g2),
                 linetype = 'dotted') +
  theme(text = element_text(size = 12),
        panel.border = element_blank(),
        strip.background =element_rect(fill = "white", colour = "white"),
        strip.text = element_text(colour = 'white'),
        plot.margin = grid::unit(c(2,2,2,2), 'mm')) +
  facet_wrap(~p) +
  ylab(" ") +
  ylab("Counts") +
  labs(x = expression(italic(g)["2"])) +
  ggtitle("a)") +
  scale_y_continuous(expand = c(0,0), limits = c(0,1800)) +
  scale_x_continuous(limits = c(-0.0015, 0.0035),
                     breaks = seq(-0.0015, 0.0035, 0.002),
                     expand = c(0,0)) 
```

```{r g2_neut_bs_hist, echo=FALSE, warning=FALSE, error=FALSE}
ggsave(g2_neut_bs_hist, filename = "g2_neut_bs_hist.tiff",
     units = "in", dpi = 300, width = 5, height = 5)
```

***

### Plot heterozygosity among marker sets
In order to visualize sMLH for all, neutral, and immune markers, create the following box-plot:
``` {r}
het_plot <- 
  ggplot(data = sealdf_reshaped, aes(x=health, y=value, fill = variable)) + 
  stat_boxplot(aes(x = health, y = value), 
               geom='errorbar', linetype=1, width=0.5) +
  geom_boxplot( aes(x = health, y = value),outlier.shape=1) +    
  geom_jitter(shape=16, position=position_jitter(0.4), size = .1) +
  theme_classic() +
  theme(legend.position = "none",
        panel.border = element_blank(),
        strip.background = element_blank(),
        text = element_text(size = 12),
        plot.margin = grid::unit(c(2,2,2,2), 'mm')) +
  xlab("Infection status") +
  ylab("sMLH") +
  ggtitle("b)") +
  scale_fill_manual(values =  c("#A79C93", "#0294A5", "#C1403D")) +
  facet_wrap(~variable, nrow = 1) +
  scale_y_continuous(limits = c(0.3, 2),
                     expand = c(0,0))
```

```{r het_plot, echo=FALSE}
ggsave(het_plot, filename = "het_plot.tiff",
     units = "in", dpi = 300, width = 5, height = 5)
```

***

\pagebreak

### Estimating heterozygosity for individual loci
As we have previously looked at genome-wide effects, it may be of interest to look for local effects. Therefore, we wanted to examine the heterozygosity for each locus. First, define the confidence interval:

```{r loci het, message=FALSE, results='hide'}
confidence_interval <- function(vector) {
  ## standard deviation
  vec_sd <- sd(vector)  
  ## sample size
  n <- length(vector)     
  ## sample mean
  vec_mean <- mean(vector)    
  ## error according to t distribution
  error <- qt((.95 + 1)/2, df = n - 1) * vec_sd / sqrt(n)       
  ## confidence interval as a vector
  result <- c("lower" = vec_mean - error, "upper" = vec_mean + error)   
  return(result)
}
```

Calculate the heterozygosity for each locus, and use a regression on infection status:
```{r}
## calcaute sMLH
het_per_locus <- apply(seals_geno, 2, sMLH)
## add factors 
df <-  cbind(sealdf, seals_geno) 
## add marker type as names to the data.frame
names(df)[6:66] <- c(paste0("I", 1:13), paste0("N", 1:48))


lm_by_loc <- lapply(1:61, function(x) {
  ## extract data of given marker x
  value <- df[,x + 5]  
  ## run linear models
  res <- summary(lm(as.numeric(df$health) ~ value))   
  conf <- confint(lm(as.numeric(df$health) ~ value))
  f <- res$fstatistic
  pf(f[1], f[2], f[3], lower=FALSE)
  out <- data.frame(beta = res$coefficients[2,1],
                    lcl = conf[2,1],
                    ucl = conf[2,2])
}) %>% 
  do.call("rbind",.) %>% 
  cbind(., data.frame(names = colnames(seals)[seq(8, ncol(seals), 2)] %>% 
                          substring(., first = 1, last = nchar(.) - 2),
                      type = c(rep("Immune", 13),rep("Neutral", 48)),
                      dummy = ""))

## order by effect size
lm_by_loc <- lm_by_loc[with(lm_by_loc, order(type, beta, decreasing = F)),]    
lm_by_loc$num <- 1:61

## create data frame to label effects
names_df <-  data.frame(label = lm_by_loc$names,
                        num = lm_by_loc$num)
``` 

Create a plot to feature each loci and their relevant effect sizes:

```{r}
het_by_loci_plot <- ggplot(lm_by_loc, aes(x = num, y = beta, col = type)) +
  geom_errorbar(aes(ymin = lcl, ymax = ucl),
                width=0.6, alpha=0.7, size = 0.7) +
  geom_point(size = 1) +
  scale_x_continuous(expand = c(0,0), breaks = 1:61, labels = names_df$label) +
  scale_y_continuous(expand = c(0,0)) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  coord_flip(xlim = c(0, 61.5), ylim = c(-0.5,0.6)) +
  scale_color_manual(values=c("#C1403D","#0294A5"), 
                     name = "",
                     breaks = c("Neutral", "Immune"),
                     labels = c("Neutral", "Immune")) +
  theme_classic() +
  xlab("") +
  ylab("Effect size") +
  theme(legend.justification = c(0,1),
        legend.position = c(0,1.05),
        legend.background = element_rect(fill = NA),
        text = element_text(size = 12),
        axis.text.y = element_text(size = 5),
        legend.text = element_text(size = 7),
        panel.border = element_blank(),
        strip.background =element_rect(fill = "white", colour = "white"),
        strip.text = element_text(colour = 'white'),
        plot.margin = grid::unit(c(2,2,2,2), 'mm')) +
  guides(color = guide_legend(
    keywidth = 0.05,
    keyheight = 0.05,
    default.unit = "inch")) +
  facet_wrap(~dummy) +
  ggtitle("c)")
```

``` {r het_by_loci_plot, echo=FALSE}
ggsave(het_by_loci_plot, filename = "het_by_loci_plot.tiff",
     units = "in", dpi = 300, width = 4, height = 6)
```

First we can look at the mean and standard error between neutral and immune loci, then we can look for local effects between the effect sizes of neutral and immune loci, using a Wilcoxon test:
```{r local effects, results='hide'}
mean(lm_by_loc$beta[1:13]) 
sd(lm_by_loc$beta[1:13])/13 

mean(lm_by_loc$beta[14:61]) 
sd(lm_by_loc$beta[14:61])/48 

wilcox.test(lm_by_loc$beta[1:13],lm_by_loc$beta[14:61])
```

***

### Create the final combined figure
To create a combination plot of all figures (as in the manuscript):
```{r create combo_plot, message=FALSE, results='hide', warning=FALSE}
## define layout of the plot
lay <- rbind(c(1,3),
             c(2,3))

## combine figures
combo_plot <- grid.arrange(g2_neut_bs_hist, 
                           het_plot, 
                           het_by_loci_plot, ncol = 3, layout_matrix = lay)
```
``` {r combo_plot, message=FALSE, echo=FALSE, results='hide'}
combo_plot
ggsave(combo_plot, filename = "combo_plot.tiff",
     units = "in", dpi = 600, width = 6, height = 6)
```

***

### Modeling effects of sMLH on bacterial infection status

To test for associations between microsatellite heterozygosity and death from bacterial infection, we constructed several alternative generalized linear mixed-models (GLMMs) incorporating relevant predictor variables and quantified their relative support using AICc weights within a multi-model inference framework.  All of the models had pup survival as a binary response variable (coded as 0 = alive and 1 = dead) and included year as a random effect to statistically control for any variation in survivorship attributable to inter-annual variation.  The following GLMMs were considered:

```{r AIC models, results='hide'}
## define models
models <- list(
  glmer(health ~  1 + (1|year), data = sealdf_year, family = 'binomial'),
  glmer(health ~  All + (1|year), data = sealdf_year, family = 'binomial'),
  glmer(health ~  Immune + (1|year), data = sealdf_year, family = 'binomial'),
  glmer(health ~  Neutral + (1|year), data = sealdf_year, family = 'binomial'), 
  glmer(health ~  1 + birthweight + (1|year), data = sealdf_year, family = 'binomial'),
  glmer(health ~  All + birthweight + (1|year), data = sealdf_year, family = 'binomial'),
  glmer(health ~  Immune + birthweight + (1|year), data = sealdf_year, family = 'binomial'),
  glmer(health ~  Neutral + birthweight + (1|year), data = sealdf_year, family = 'binomial')) 
names(models) <- paste0("m", 1:length(models))

## model selection using delta AICc
kableExtra::kable(AICcmodavg::aictab(models, second.ord = T), booktabs = TRUE,
                  longtable = FALSE, caption = "Model selection")

# Model averaged paramaters, one by one
AICcmodavg::modavg(models, parm = "All")
AICcmodavg::modavg(models, parm = "Immune")
AICcmodavg::modavg(models, parm = "Neutral")
AICcmodavg::modavg(models, parm = "birthweight")
```

These included 'null models' without any genetic effects (models i and v) as well as models that included sMLH combined over all loci or calculated separately for the neutral versus immune loci.  Models v to viii also included pup birth weight (in kg) to incorporate any potential effects of body size on survivorship. All of the models were specified using the glmer function of the package “lme4” with a binomial error structure.[^7] Using the R package `AICcmodavg`, the most parsimonious model was selected based on the delta AICc value, which compares weights as a measure of the likelihood of a particular model.[^8] The best supported model has $\Delta$ AICc = 0 and a difference of two or more units was applied as a criterion for choosing one model over a competing model.[^9]

[^7]: Bates, D., Maechler, AICcmodavg::aictab(models, second.ord = T)M., Bolker, B., & Walker, S. (2014). Fitting linear mixed-effects models using lme4. arXiv preprint arXiv:1406.5823.

[^8]: Mazerolle, M. J., & Mazerolle, M. M. J. (2017). Package ‘AICcmodavg’. R package.

[^9]: Anderson, D. R., & Burnham, K. P. (2002). Avoiding pitfalls when using information-theoretic methods. The Journal of Wildlife Management, 912-918.



***

# Supplementary Data
#### (a) Patterns of allelic richness and cross-amplification

We tested for patterns in allelic richness among markers (i.e. immune vs neutral), developmental source (i.e. designed for Antarctic fur seals, phocids or otariids). Secondly, we evaluate the cross-amplification success of loci in two other species of pinnipeds, namely the Grey seal and Nothern Elephant seal. 

```{r heat_map, message=FALSE, warning=FALSE, results='hide'}
## read and format genotypes
library(readxl)
heatmap_df <- readxl::read_xlsx("data/data.xlsx")[, c(3, 8:ncol(seals))]

## randomly select six individuals per species
heatmap_df <- heatmap_df[c(sample(which(heatmap_df[["Species"]] == "Fur seal"), 
                                  size = 6, replace = F),
           sample(which(heatmap_df[["Species"]] == "Grey seal"), 
                  size = 6, replace = F),
           sample(which(heatmap_df[["Species"]] == "Northern Elephant seal"), 
                  size = 6, replace = F)),]

## extract geno
marker_geno <- apply(heatmap_df[,-1], 2, as.character)

## get loci names
loci_names <- colnames(marker_geno)[seq(1, ncol(marker_geno), 2)] %>% 
  substring(., first = 1, last = nchar(.) - 2)

## define a vector of immune marker names
immune_marker_names <- c("Agi01", "Agi02", "Agi03", "Agi04",
                         "Agi05", "Agi06", "Agi07", "Agi08",
                         "Agi09", "Agi10", "Agi11", "Agt10", "Agt47")

## collapse information for each locus in one column
marker_geno <- lapply(seq(1, ncol(marker_geno), 2), function(x) {
  marker_geno[,x:(x + 1)] %>% 
    apply(., 1, paste0, collapse = "/")
}) %>% 
  do.call("cbind",.) %>% 
  ## rename loci
  set_colnames(x = ., value = paste0("Locus", 1:61))

## set missing data to NA
marker_geno[which(marker_geno == "NA/NA")] <- NA

## convert to GENIND object
genind <- adegenet::df2genind(marker_geno, ploidy = 2, sep = "/", 
                              pop = heatmap_df[["Species"]] %>% as.factor)

## convert to GENPOP
genpop <- adegenet::genind2genpop(genind)
```

```{r}
heatmap_df <- lapply(levels(genpop@loc.fac), function(i) {
  df.temp <- genpop@tab[,which(genpop@loc.fac == i)]   ## fetch data 
  if (is.null(dim(df.temp))) {
    df.temp[df.temp > 0] <- 1
    df.temp[df.temp == 0] <- 0
    
  } else {
    df.temp <-   apply(df.temp, 2, function(x) ifelse(x > 0, 1, 0)) %>%
      ## presence/absence of allele
    rowSums(na.rm = T) ## count alleles
  }
  # return results
  return(data.frame(Species = names(df.temp),
                    Locus = i,
                    Alleles = df.temp))
}) %>% 
  do.call("rbind", .)

## set zero to NA
heatmap_df[["Alleles"]][which(heatmap_df[["Alleles"]] == 0)] <- NA

heatmap_df[["Locus"]] <- factor(heatmap_df[["Locus"]], labels = loci_names)
heatmap_df[["Type"]] <- 'Neutral'
heatmap_df[["Type"]][which(heatmap_df[["Locus"]] %in% immune_marker_names)] <- 'Immune'

## sort by species
heatmap_df[["Species"]] <- factor(heatmap_df[["Species"]],
                          levels = c("Fur seal", "Grey seal", 
                                     "Northern Elephant seal"),
                          labels = c("Antarctic fur seal", "Grey seal", 
                                     "Northern Elephant seal"))

## define colors for marker types
col_key <- ifelse(levels(heatmap_df[["Locus"]]) %in% immune_marker_names, 
                  "#C1403D", "#0294A5")
```

```{r}
plot <- ggplot(data = heatmap_df, aes(x = Species, y = Locus, fill = Alleles)) +
  theme_classic() +
  geom_tile(colour = "Black", size = .75) + 
  scale_fill_viridis_c(name = "Alleles", na.value = "Grey75") +
  scale_x_discrete(expand = c(0,0)) +
  theme(
    plot.margin = margin(t = 5, r = 25, b = 5, l = 15, unit = "mm"),
    legend.position = c(1,1), 
    legend.justification = c(0, 1),
    legend.direction = "vertical",
    legend.margin = margin(0,0,0,5, "mm"),
    axis.text.y = element_text(hjust = 0, colour = col_key, size = 8),
    axis.line.x = element_blank(),
    axis.text.x = element_text(size = 10)) + 
  xlab("Species") +
  ylab("") 

ggsave(plot,
       filename = 'HeatmapLoci.tiff',
       width = 6,
       height = 9,
       units = "in",
       dpi = 300)
```

```{r, echo=FALSE}
plot
```

The heatmap above shows several patterns which are tested statistically next.

```{r, warning=FALSE, message=FALSE}
## get raw data again
genotypes_raw <- readxl::read_xlsx("data/data.xlsx", skip = 1)[, c(3, 8:ncol(seals))]

## extract genotypes
marker_geno <- apply(genotypes_raw[,-1], 2, as.character)

## collapse information for each locus in one column
marker_geno <- lapply(seq(1, ncol(marker_geno), 2), function(x) {
  marker_geno[,x:(x + 1)] %>% 
    apply(., 1, paste0, collapse = "/")
}) %>% 
  do.call("cbind",.) %>% 
  ## rename loci
  set_colnames(x = ., value = paste0("Locus", 1:61))

## set missing data to NA
marker_geno[which(marker_geno == "NA/NA")] <- NA

## create GENIND for Antarctic fur seal alone
genind_afs <- adegenet::df2genind(marker_geno[1:78,], ploidy = 2, sep = "/")
```

```{r, warning=FALSE, results='hide', message=FALSE}
## extract allele numbers for both marker types
immune_afs <- genind@loc.n.all[1:13]
mean(immune_afs)
sd(immune_afs)

neutral_afs <- genind@loc.n.all[14:61]
mean(neutral_afs)
sd(neutral_afs)

## compare marker types
wilcox.test(immune_afs, neutral_afs, paired = F)
```

```{r, warning=FALSE, results='hide', message=FALSE}
## compare neutral markers by origin
neutral_afs <- genind@loc.n.all[14:22]
mean(neutral_afs)
sd(neutral_afs)

neutral_others <- genind@loc.n.all[23:61]
mean(neutral_others)
sd(neutral_others)

## compare by marker
wilcox.test(neutral_afs, neutral_others, paired = F)
```

```{r, warning=FALSE, results='hide', message=FALSE}
## cross-amplification
immune <- dplyr::filter(heatmap_df, Species != "Antarctic fur seal", 
                        Type == "Immune")[["Alleles"]]
immune <- ifelse(is.na(immune), 0, 1) # check if amplified
mean(immune) ## cross-amplification rate

neutral <- dplyr::filter(heatmap_df, Species != "Antarctic fur seal")[27:44, "Alleles"]
neutral <- ifelse(is.na(neutral), 0, 1) # check if amplified 
mean(neutral) ## cross-amplification rate


wilcox.test(neutral, immune, paired = F)
```

#### (b) FDR Correction
After receiving HWE values, apply a false discovery rate correction for a table of p-values.
```{r FDR, message=FALSE, warning=FALSE, results='hide'}
pval <- 
  read.table("data/pvalues.txt",header = F, sep = ",") %>% 
  as.vector() %>% .[[1]]
qobj <- qvalue(pval)
qvalues <- qobj$qvalues
pi0 <- qobj$pi0
lfdr <- qobj$lfdr

summary(qobj)
df <- data.frame(p = qobj$pvalues,
                 q = qobj$qvalues)
## view(df)
```

#### (c) Sensitivity of loci number on estimates of *g*~2~

Here, we repeat the estimation of *g*~2~ for each marker type and for the entire dataset

```{r g2_neut_resampled, eval=FALSE, message=FALSE, warning=FALSE, results='hide'}
g2_neutral_resampled <- pbapply::pblapply(seq(4, 48, 4), function(x) { 
  subs <- lapply(1:100, function(y) {
    rand <- sample(1:48, x, replace = FALSE)
    loci <- neutral_markers[, rand]
    g2 <- g2_microsats(loci, nperm = 0, nboot = 9999, verbose = F)
    df <- data.frame(bs = g2$g2_boot,
                     lcl = g2$CI_boot[[1]],
                     ucl = g2$CI_boot[[2]],
                     g2  = g2$g2,
                     p = g2$p_val)
    return(df[1,])
  }) %>% do.call("rbind", .)
  return(data.frame(g2 = mean(subs$g2),
                    lcl = confidence_interval(subs$g2)[1],
                    ucl = confidence_interval(subs$g2)[2]))
}) %>%  do.call("rbind", .)
g2_neutral_resampled$loci <- seq(4, 48, 4)
```

```{r g2_neut_resampledb, eval=FALSE, message=FALSE, warning=FALSE, results='hide'}
## load saved dataset
load("data/g2_neutral_resampled.RData")

g2_neutral_resampled_plot <-
  ggplot(data = g2_neutral_resampled, aes(x = loci, y = g2)) +
  geom_line() +
  geom_point(size = 1.5) +
  geom_errorbar(aes(ymin = lcl,
                    ymax = ucl),
                width = 0.8, alpha = 0.7, size = 0.8, colour = "black") +
  geom_hline(yintercept = 0, linetype = "dotted") +
  theme_classic() +
  theme(legend.position = "none",
        panel.border = element_blank(),
        strip.background = element_blank(),
        text = element_text(size = 12),
        aspect.ratio = 1,
        axis.title.y = element_text(face = "italic"),
        plot.margin = grid::unit(c(2,2,2,2), 'mm')) +
  xlab("Number of loci") +
  labs(y = expression(italic(g)["2"])) +
  scale_x_continuous(expand = c(0,0), limits = c(0, 50))
```

```{r g2_neut_resampledc, eval=FALSE, message=FALSE, warning=FALSE, results='hide'}
g2_neutral_resampled_plot
```

## Microsatellite development

### Identifying microsatellites within the transcriptome

Short tandem repeats were identified within the Antarctic fur seal transcriptome assembly using the script `misa.pl`. The required initiation file called `misa.ini` is available in the folder `data` and defines the minimum number of five repeats for di-, tri- and tetranucleotide motifs. *As already mentioned before, [MISA](http://pgrc.ipk-gatersleben.de/misa/misa.html) needs to be downloaded by the user*

```bash
# identify microsats
perl misa.pl arc_gaz_transcriptome.fasta
```

The code above generates one ouput file `arc_gaz_transcriptome.fasta.misa` containing a total of `r length( readLines("data/arc_gaz_transcriptome.fasta.misa")) - 1` microsatellites found within the transcriptome. The resulting data table was subsequently reformatted for further filtering steps. 

```{r, eval=T}
## read the data
data <- readLines("data/arc_gaz_transcriptome.fasta.misa")

## create a matrix
microsats_table <- matrix(ncol = length(strsplit(data[1], split = "\t")[[1]]),
                          nrow = length(data))
for (i in 1:length(data)) {
microsats_table[i,1:length(strsplit(data[i], split = "\t")[[1]])] <- 
  strsplit(data[i], split = "\t")[[1]]
}
microsats_table <- as.data.frame(microsats_table[2:nrow(microsats_table),])

## add column names
names(microsats_table) <-
  c('contig.name', 'ssr.no','ssr.type','ssr.seq','ssr.size','ssr.start','ssr.end') 

## read contig length information
contig_length <- read.table("data/transcriptlength.txt", header = T)
names(contig_length) <- c('MatchID', 'contig.length')

## set to character
microsats_table[["contig.name"]] <-
  as.character(microsats_table[["contig.name"]])
contig_length[["MatchID"]] <-
  as.character(contig_length[["MatchID"]])

## correct 'MatchID' for cross-referencing
microsats_table[["MatchID"]] <- NA
for (i in 1:nrow(microsats_table)) {
  # discard chunk following the underscore (e.g. 4708387_length... becomes 4708387)
  microsats_table$MatchID[i] <-
   strsplit(microsats_table$contig.name[i],split = "_")[[1]][1]
  microsats_table$MatchID[i] <-
    strsplit(microsats_table$MatchID[i],split = " ")[[1]][1]
} 

## correct contig_length
for (i in 1:nrow(contig_length)) {
  # remove everything after the underscore, see above
  contig_length$MatchID[i] <-
    strsplit(contig_length$MatchID[i],split = "_")[[1]][1]
}

## merge data frames
microsats_table <-
dplyr::left_join(microsats_table,contig_length, by = "MatchID")

## some data class conversions
microsats_table[["ssr.start"]] <-
  as.numeric(as.character(microsats_table[["ssr.start"]])) 
microsats_table[["ssr.end"]] <-
  as.numeric(as.character(microsats_table[["ssr.end"]]))
microsats_table[["contig.length"]] <-
  as.numeric(as.character(microsats_table[["contig.length"]]))

## export Supplementary Table S1
write.csv(microsats_table, file = "data/Supplementary Table S1.csv", row.names = F)

## summary of microsatellite types
summary(microsats_table$ssr.type)
```

### Filtering Microsatellites

Among the `r nrow(microsats_table)` identified microsatellites, there are some compound microsatellites as well as repeats that do not offer adequate flanking sites for primer design. These were discarded next.

```{r}
## remove compound microsats
microsats_table <- 
  subset(microsats_table, microsats_table[["ssr.type"]] != 'c') 
microsats_table <- 
  subset(microsats_table, microsats_table[["ssr.type"]] != 'c*')

## selection based on flanking sites
microsats_table[["temp"]] <-
  rep(1,nrow(microsats_table)) # flag for removal
for (i in 1:nrow(microsats_table)) {
  # inspect flanking site upstream
if (microsats_table[["ssr.start"]][i] <= 100) {
      microsats_table[["temp"]][i] <- 0
      # inspect flanking site downstream
} else if ((microsats_table[["contig.length"]][i] -
            microsats_table[["ssr.end"]][i]) <= 100) {
    microsats_table[["temp"]][i] <- 0
  }
}  

## Remove flagged microsats
microsats_table <- subset(microsats_table, microsats_table[["temp"]] != 0) 
microsats_table <- microsats_table[,1:9]
```

After the above filtering `r nrow(microsats_table)` microsatellites were retained. Now, we selected microsatellites that are associated to immunity based on [Gene Ontology](http://www.geneontology.org/) Gene annotations.

```{r}
## Keywords including 'immune*'
annotation <- readLines("data/arc_gaz_transcriptome_annotations.txt")
annotations <- matrix(ncol = 18, nrow = length(annotation))

## fill table
for (i in 1:nrow(annotations)) {
annotations[i, 1:length(strsplit(annotation[i], split = "\t")[[1]])] <-
  strsplit(annotation[i], split = "\t")[[1]]
}
annotations <- data.frame(annotations)[-1,]
annotations <- annotations[, c(1,14:18)]
names(annotations) <- c('MatchID','goTerm','cellular.components','biological.processes',
                        'molecular.functions','keywords')

annotations[["MatchID"]] <- as.character(annotations[["MatchID"]])
for (i in 1:nrow(annotations)) {
  annotations[["MatchID"]][i] <-
    strsplit(annotations[["MatchID"]][i],split = "_")[[1]][1]
} 

annotations.extd <- dplyr::left_join(microsats_table, annotations,by = 'MatchID')

immuneTable2 <- data.frame(annotations.extd) %>% # Check for matches with keywords
  dplyr::filter(grepl('immun*', keywords)) 

ImmuneMarker_Keywords <- immuneTable2 # 13 within just keywords
```


```{r, echo=FALSE, eval=FALSE}
## overview
ImmuneMarker_whole_file$keywords
kableExtra::kable(ImmuneMarker_Keywords[,c(1, 14)] , 
                  booktabs = TRUE, longtable = TRUE,
                  caption = "Keywords 'Immune'") %>%
kableExtra::kable_styling(latex_options = c("hold_position", "repeat_header"))
```

For 13 microsatellites, we found a match to the term 'immun*' under the keywords of the GO annotations. To increase the number of suitable microsatellites, we repeated the initial search to all categories of the GO annotations with an extended list of search terms shown below.

```{r}
## define list of keywords
immune <- c('immun*', 'antigen', 'chemokine', 'T cell',
            'MHC', 'Antibody', 'histocompatibility',
            'Interleukin', 'Leucocyte', 'Lymphocyte')

immuneLines <- NULL
for (i in immune) { 
  immuneLines <- c(immuneLines, annotation[grep(i, annotation, ignore.case = T)])
}

immuneTable <- matrix(ncol = 18, nrow = length(immuneLines))

for (i in 1:length(immuneLines)) {
immuneTable[i,1:length(strsplit(immuneLines[i], split = "\t")[[1]])] <-
  strsplit(immuneLines[i], split = "\t")[[1]]
}
immuneTable <- data.frame(immuneTable)[,c(1,10,14:18)] 
names(immuneTable) <-
  c('MatchID','geneID','goTerm','cellular.components','biological.processes',
                        'molecular.functions','keywords') 

immuneTable[["MatchID"]] <-
  as.character(immuneTable[["MatchID"]])
for (i in 1:nrow(immuneTable)) {
  immuneTable[["MatchID"]][i] <-
    strsplit(immuneTable[["MatchID"]][i],split = "_")[[1]][1]
} 

ImmuneMarker_whole_file <-
  unique(dplyr::inner_join(microsats_table, immuneTable, by = "MatchID")) 

## write to file
write.csv2(ImmuneMarker_whole_file, file = "data/immune_microsats_raw.csv", row.names = F)
```

The extended search yielded a total of `r nrow(ImmuneMarker_whole_file)` microsatellites. The entire list is shown below.   

```{r, echo=FALSE}
kableExtra::kable(ImmuneMarker_whole_file[, c(1, 4, 6:7, 10)],
                  booktabs = TRUE, longtable = TRUE,
                  col.names = c('Contig', 'Motif', 'Start', 'End', 'Gene ID'),
                  caption = "Annotated microsatellites") %>%
kableExtra::kable_styling(latex_options = c("hold_position", "repeat_header"))
```

\newpage

### Designing primers

For all of the `r nrow(ImmuneMarker_whole_file)` we developed oligonucleotide primers using the primer design tool [primer3](http://primer3.sourceforge.net/releases.php). In order to use the command line interface, the list of microsatellites should be re-formatted accordingly. 

```{r, eval=T}
## list of microsatellites 
data <- read.csv(file = "data/immune_microsats_raw.csv", sep = ';')[,1:7] 
names(data) <- c('ID','SSR nr.','SSR type','SSR','size','start','end')
data[["ID"]] <- as.character(data[["ID"]])

for (i in 1:nrow(data)) {
  if ((nchar(data[["ID"]][i]) > 20)) {
    data[["ID"]][i] <- 
      paste0(strsplit(data[["ID"]][[i]],split = "_")[[1]][1]," ", 
             strsplit(data[["ID"]][[i]],split = "_")[[1]][2]," ",
             strsplit(data[["ID"]][[i]],split = "_")[[1]][3]," ",
             strsplit(data[["ID"]][[i]],split = "_")[[1]][4],"_",
             strsplit(data[["ID"]][[i]],split = "_")[[1]][5],"_",
             strsplit(data[["ID"]][[i]],split = "_")[[1]][6],"_",
             strsplit(data[["ID"]][[i]],split = "_")[[1]][7])
  }
}
write.table(row.names = FALSE, quote = FALSE,x = data,
            sep = "\t",file = 'data/arc_gaz_transcriptome.fasta.misa2')
```

### Invoke primer3 for primer design

```bash
perl p3_in_fur_seal.pl arc_gaz_transcriptome.fasta.misa2
primer3_core <arc_gaz_transcriptome.fasta.p3in> arc_gaz_transcriptome.fasta.p3out
``` 

### Overview of initially tested microsatellites

The table below summarises the results of testing 96 primers on 12 Antarctic fur seal indiviudals. See the manuscript for further details. 

\newpage
\blandscape

```{r, echo=FALSE}
## overview
results <- read.csv("data/Supplementary Table S2.csv")[,-c(2, 5:7)]
results$Forward.primer_.5..3._._M13..21..tail <- stringr::str_remove(
  string = results$Forward.primer_.5..3._._M13..21..tail,
  pattern = "TGTAAAACGACGGCCAGT")
kableExtra::kable(results, booktabs = TRUE, 
                  longtable = TRUE, caption = "Overview microsatellite testing. Primers used in the present study are named Agi01-Agi11", 
                  col.names = c("Contig", "Gene ID", "Motif", 
                                'Marker ID', 
                                "Forward primer 5'-3'",
                                "Reverse primer 5'-3'",
                                "PCR result")) %>%
kableExtra::kable_styling(latex_options = c("hold_position", "repeat_header"))
```

\elandscape
***

```{r, echo = FALSE}
sessioninfo::session_info()
```