---
title: 'Script to assess genetic differentiation of polar cod using nine microsatellite loci'
author: "Sarah Maes"
date: "09/04/2019"
output:
  html_document:
    df_print: paged
  pdf_document: default
---

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

#### This script performs analyses as described in: ####
#### Maes, S.M., Christiansen, H., Mark, F.C., Lucassen, M., Van de Putte, A.P., Volckaert, F.A.M., Flores, H. (2021) High gene flow in polar cod (Boreogadus saida) from West-Svalbard and the Eurasian Basin. Journal of Fish Biology ####

We assume here that the relevant packages are already installed. If this not the case, install them using install.packages() or another function, as you prefer. These are the major packages used here:

```{r libraries, include = FALSE}
packages <- c('adegenet', 'diveRsity', 'pegas', 'poppr')
lapply(packages, require, ch = T)
```

If you use any of these, don't forget to cite them:

```{r pkgs, echo = FALSE}
lapply(packages, citation)
```

In addition we use the package 'here' for relative file paths. We now load the microsatellite genotypes and convert a Genepop file to a genind object.

```{r load}
library(here)
library(adegenet)
pcod <- read.genepop("genotypesBS.gen", 3)
pcod
```

Add population names and check the data (number of individuals, group sizes, number of alleles, missing data, observed and expected heterozygosity).

```{r popnames}
pcod@type
levels(pcod$pop) <- c('Nansen Basin', 'Amundsen Basin', 'Hornsund',
                      'Kongsfjorden', 'Billefjorden', 'Sophia Basin')
summary(pcod)
```

Each locus is analysed independently to avoid bias introduced by a single locus. The number of alleles, expected and observed heterozygosity, fixation index FST with 95% confidence intervals and inbreeding coefficient FIS with 95% confidence intervals at each locus and per putative population is determined using diVersity package v1.9.90 (Keenan et al. 2013). You may need to install the package HWxtest for this function to work. Use the code below, if necessary:
install.packages('devtools')
devtools::install_github('wrengels/HWxtest', subdir = 'pkg')

Now the per locus indices can be calculated.

```{r bylocus}
library(diveRsity)
basicStats(infile = 'genotypesBS.gen', outfile = 'diVersity_basicStats_output',
           fis_ci = TRUE, ar_ci = FALSE, fis_boots = 1000, ar_boots = NULL,
           mc_reps = 999, rarefaction = TRUE, ar_alpha = 0.05,
           fis_alpha = 0.05)
```

For each locus, departures from Hardy-Weinberg proportions are assessed after 10,000 permutations with package pegas v0.10. Departures from Hardy-Weinberg for each putative population are assessed in the diveRsity package using the basicstats function.

```{r hwtest}
library(pegas)
hw.test(pcod, B=10000)
```

Nei's GST, Jost's D and Hedrick's G''ST for each locus with 95% confidence intervals are also calculated with the diVersity package after bootstrapping (10,000 permutations). Population differentiation is measured using pairwise FST comparisons according to Weir & Cockerham (1984) with the same package and 95% confidence intervals after bootstrapping (10,000 permutations). 

```{r gst}
library(diveRsity)
diffCalc(infile = 'genotypesBS.gen', outfile = 'diVersity_diffCalc_output', fst = TRUE, bs_locus = TRUE, pairwise = TRUE, bs_pairwise = TRUE, boots = 10000)

```

The number of private alleles are estimated per putative population with the poppr package v2.8.1.

```{r par}
library(poppr)
private_alleles(pcod, form = alleles ~ ., report = "table",
                level = "population", count.alleles = TRUE, drop = FALSE)
```

Allelic richness after rarefaction and private allelic richness after rarefaction were calculated with HP-Rare software v1.1 outside the R environment.

A Discriminant Analysis of Principal Components (DAPC) is conducted with package adegenet v2.0.0 with the aim of assessing genetic differentiation and clustering patterns. The analysis here is based on 15 principal components and two discriminant functions. 

```{r dapc}
pcod.dapc <- dapc(pcod, n.da = 2, n.pca = 15)
pcod.dapc$ind.coord
scatter(pcod.dapc, ratio.pca = .2, scree.pca = TRUE, ratio.da = .2, 
        col = c("#009E73","#E69F00", "#0072B2","#56B4E9",  "#D55E00", "#CC79A7" ),
        xax = 1, yax = 2, solid = .6, pch = 17:22, lwd = 2, axesell = FALSE,
        inset.solid = 0.75, clabel = FALSE, cstar = 0, cell = 0, cex = 1.8,
        leg = T, cleg = 0.75, posi.pca = "topleft", posi.leg = 
        'topright')
```

The Bayesian clustering aproach to assign individuals to populations is performed with STRUCTURE v2.3.4 and not in R. 

We test for isolation-by-distance (IBD) using package adegenet. IBD is tested using the Mantel test between a matrix of geographic and genetic distances. We first convert the genind data back to a genpop format.


```{r ibd}
library(adegenet)
library(ade4)
IBD.coord <- read.table("IBD_coord.txt", header = TRUE)
IBD_data <- genind2genpop(pcod)
Dgen <- dist(IBD_data)
Dgeo <- dist(IBD.coord)
mantel.randtest(Dgen, Dgeo)
``` 

We assess linkage disequilibrium in GENEPOP v4.6, not in R.

We test for statistical power of the nine microsatellite loci with POWSIM v4.1. For this analysis, we need allele frequencies. These are calculated with the R package PopGenReport v3.0.4. 

```{r}
library(PopGenReport)
library(knitr)
allele.dist(pcod)
```

