10_analyze_novembre2008_full.Rmd

---
title: "Genes and Geography in Europe"
output:
  html_document: default
  pdf_document: default
---

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

In this exercise, the aim is to reproduce parts of the analysis of one of the seminal papers in human population genetics, ["Genes Mirror Geography within Europe"](https://www.nature.com/articles/nature07331), Novembre et al 2018 Nature. We need to load two R packages for this: 

```{r env, message=FALSE, warning=FALSE}
library(tidyverse)
library(vegan)
library(maps)
```

<br />
The dataset for reproducing the main figures is available from the [Novembre lab Github repository](https://github.com/NovembreLab/Novembre_etal_2008_misc). Here, we will use a minimally reformatted versions of the main PCA results table. There are two expected outcomes of this exercise:

**Create a replication of the PCA plot in Figure 1, starting from the PCA results table**

To achieve this, you should carry out the following steps:

* Read the [dataset](https://www.dropbox.com/s/xloe4h44rokni49/novembre_2008_pca.tsv?dl=0) containing PCA coordinate and geographic locations, as well as the [color palette table](https://www.dropbox.com/s/l8avwtbb5h3sx55/novembre_2008_colours.tsv?dl=0) into R

* Rotate the original PCA coordinates to maximize the correlation between PCs and geographic locations. This can be achieved using the `procrustes()` function of the `vegan` R package. The function takes two matrix objects as input, the first one the target matrix (i.e. the geographic coordinates), and the second one the matrix to be rotated (i.e. the PCA result). More details on how to use the function can be found in the help page by typing `?procrustes`

* Visualize the resulting rotated PCA using `ggplot2`. Start by trying to replicate the overall PCA plot of samples (using `geom_text()` with the given population labels), with the color scale in the given palette table (using `scale_colour_manual()`). Add population median positions using summary tables produced by `group_by()` and `summarise()`. 

* Adding rotated axes and their labels requires a bit more processing, with one possible solution being rotating points at the ends of the original axes into the new coordinates by using the `predict()` function on the procrustes result object. 

<br />
**Predict the geographic location of samples from a chosee country using their PCA position**

To achieve this, you can follow the linear regression approach from Novembre et al:

* Split your data into a test set containing only samples from the country of interest, and a training set of the remaining individuals.

* Create two linear regression models, one for longitude and one for latitude as outcomes. For each model, use both rotated PC1 and PC2, as well as their interaction as predictors. The syntax for such a model in R is `lm(response ~ predictor1 * predictor2, data = data)`.

* Obtain the predicted longitude and latitude for the test samples using the `predict()` function on both models. 

* Plot the actual and predicted locations of the test samples on a map. You can find a simple example of how to plot spatial data in `ggplot2` towards the end of the chapter "[Data visualization](https://r4ds.had.co.nz/data-visualisation.html)" of the "R for data science" book.

<br />
The following sections of this document show an example solution for these tasks

## Preparing the dataset

First we'll set up the environment. We need the `tidyverse` packages, as well as the `vegan` package to perfrorm the procrustes analysis

```{r import}
pca <- read_tsv("novembre_2008_pca.tsv")

colors <- read_tsv("novembre_2008_colours.tsv",col_names = c("country", "color"))
```

<br />
The variable `pca` contains the results of the principal component analysis of 1,387 European individuals from Novembre et al 2008. Let's have a look at the first few lines of the dataset:

```{r data1, echo=FALSE}
knitr::kable(
  head(pca), booktabs = TRUE
)
```

<br />
The color palette used for the countries is stored as a table in the variable `colors`:

```{r data2, echo=FALSE}
knitr::kable(
  head(colors), booktabs = TRUE
)
```

## Rotating the PCA

In order to find the rotation for the PCA that aligns best with the geographic location of the samples, we will use the `procrustes()` function from the package `vegan`. We first need to convert both the principal components and geographic locations of the samples into matrix format.

```{r procrustes1}
loc <- pca %>%
    select(longitude, latitude) %>%
    as.matrix()

pca1 <- pca %>%
    select(PC1, PC2) %>%
    as.matrix()
```
   
<br />
Now we carry out the procrustes analysis to find the optimal rotation

```{r procrustes2}
rot <- procrustes(loc, pca1)
```

<br />
We can examine the result by using base R `plot` function on the resulting object

```{r procrustes3}
plot(rot)
```

<br />
The rotated PC coordinates can be extracted from the `Yrot` element of the result object. Here we add two new columns for the rotated PCs to our `pca` object

```{r procrustes4}
pca <- pca %>%
    mutate(PC1_rot = rot$Yrot[,1], 
           PC2_rot = rot$Yrot[,2])
```

## Recreating Figure 1 of Novembre et al 2008

With the rotated PCA results calculated, everything is in place for creating the figure. First we set up the color scale for the countries

```{r viz1}
colorscale <- colors$color
names(colorscale) <- colors$country
```

<br />
Next, we calculate the median PC values for each country

```{r viz2}
pca_summary <- pca %>%
    group_by(country, alabels) %>%
    summarise_at(c("PC1_rot", "PC2_rot"), median)
```

<br />
To add rotated axes lines, we take the coordinates of the four points at the end of the axis ranges, and rotate them into the new coordinate space using the `predict` function on the procrustes output object

```{r viz3}
lim_pc1 <- extendrange(pca$PC1)
lim_pc2 <- extendrange(pca$PC2)

axis <- rbind(c(lim_pc1[1], 0), c(lim_pc1[2], 0), c(0, lim_pc2[1]), c(0, lim_pc2[2]))
axis_rot <- predict(rot, axis, truemean = FALSE)
axis_rot_df  <- tibble(x = axis_rot[c(1, 3),1], 
                       xend = axis_rot[c(2, 4),1], 
                       y = axis_rot[c(1, 3),2], 
                       yend = axis_rot[c(2, 4), 2])
```

<br />
Similarly, we calculate rotated coordinates for two axis labels, in a position shifted towards the origin from the axis end point

```{r viz4, message=FALSE, warning=FALSE}
labels <- rbind(c(lim_pc1[2] * 0.85, 0), c(0, lim_pc2[2] * 0.85))
labels_rot <- predict(rot, labels, truemean = FALSE) %>%
    as_tibble() %>%
    mutate(labels = c("PC1", "PC2"))
colnames(labels_rot)[1:2] <- c("PC1_rot", "PC2_rot")
```

<br />
Finally, we need to calculate the rotation angle from the rotation matrix for the axis label text

```{r viz5}
theta1 <- acos(rot$rotation[1,1]) * 180 / pi
```
Our rotation angle is `r round(120-theta1, 2)`, closely replicating the -16 degrees from Novembre et al.

<br />
With all objects in place we can now build the final plot

```{r viz6}
p <- ggplot(pca, aes(x = PC1_rot, 
                     y = PC2_rot))
p +
    geom_segment(aes(x = x, 
                     xend = xend, 
                     y = y, 
                     yend = yend), 
                 size = 0.25, 
                 arrow = arrow(ends = "both", 
                               angle = 10, 
                               type = "closed", 
                               length = unit(0.02, "npc")), 
                 data = axis_rot_df) +
    geom_point(size = 10, 
               colour = "white", 
               data = labels_rot) +
    geom_text(aes(label = labels), 
              angle = c(theta1 + 180, theta1 + 270), 
              data = labels_rot) +
    geom_text(aes(label = alabels, 
                  colour = country), 
              alpha = 0.8, 
              size = 2) +
    geom_point(aes(colour = country), 
               size = 7, 
               data = pca_summary) +
    geom_text(aes(label = alabels), 
              size = 3, 
              data = pca_summary) +
    scale_color_manual(name = "Country", 
                       values = colorscale) +
    coord_equal() +
    theme_void() +
    theme(legend.position = "none")
```


## Predicting location of origin from the PCA position

In this last section, we will follow the paper and build a predictor for an individual's geographic location from its PCA position, using linear regression. As test case, we take all individuals of a country of choice, remove them from the dataset and use the remaining individuals to build the model. Then, we use the built model to predict the location of the test individuals. 

<br />
First we split our dataset into training and test data
```{r pred1}
test_country <- "Germany"
d_train <- filter(pca, country != test_country)
d_test <- filter(pca, country == test_country)
```

<br />
Next, we carry out two linear regressions, modelling both longitude and latitude as functions of the rotated PCs and including an interaction term
```{r pred2}
lm_long <- lm(longitude ~ PC1_rot * PC2_rot, data = d_train)
lm_lat <- lm(latitude ~ PC1_rot * PC2_rot, data = d_train)
```

<br />
We can examine the fit of the linear models using the `summary()` function
```{r pred2a}
summary(lm_long) 
summary(lm_lat) 
```
Both models result in significant fits with the largest estimated effect sizes for the rotated PC1 / PC2 coordinates on longitude / latitude respectively.

<br />
Now we use the `predict()` function to infer the predicted geographic location for the test samples, and reshape the results into a format useful for plotting 
```{r pred3}
d1 <- select(d_test, sampleId, longitude, latitude) %>%
    mutate(longitude_pr = predict(lm_long, d_test), 
           latitude_pr = predict(lm_lat, d_test), 
           longitude_orig = longitude, 
           latitude_orig = latitude) %>%
    select(-longitude:-latitude) %>%
    pivot_longer(-sampleId) %>%
    separate(name, into = c("variable", "group")) %>%
    pivot_wider(names_from = variable, 
                values_from = value)
```

<br />
Finally, we can plot the inferred versus actual sample locations using the basic map plotting functionality in `ggplot2`
```{r pred4}
w <- map_data("world")

lim_long <- extendrange(pca$longitude)
lim_lat <- extendrange(pca$latitude)

pal <- c("black", "red")
names(pal) <- c("orig", "pr")
sh <- c(16, 4)
names(sh) <- c("orig", "pr")

p <- ggplot(w)
p +
    geom_polygon(aes(long, 
                     lat, 
                     group = group), 
                 fill = "white", 
                 color = "grey", 
                 size = 0.25) +
    geom_point(aes(x = longitude, 
                   y = latitude, 
                   colour = group, 
                   shape = group), 
               size = 2, 
               data = d1) +
    scale_color_manual(values = pal) +
    scale_shape_manual(values = sh) +
    coord_quickmap(xlim = lim_long, ylim = lim_lat)
```