README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  warning = FALSE,
  message = FALSE
)
```

# mapmixture

<!-- badges: start -->
[![R-CMD-check](https://github.com/Tom-Jenkins/mapmixture/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/Tom-Jenkins/mapmixture/actions/workflows/R-CMD-check.yaml)
[![cran](https://www.r-pkg.org/badges/version/mapmixture)](https://CRAN.R-project.org/package=mapmixture)
[![download](https://cranlogs.r-pkg.org/badges/grand-total/mapmixture)](https://CRAN.R-project.org/package=mapmixture)
<!-- [![Codecov test coverage](https://codecov.io/gh/Tom-Jenkins/mapmixture/branch/main/graph/badge.svg)](https://app.codecov.io/gh/Tom-Jenkins/mapmixture?branch=main) -->
<!-- badges: end -->

`mapmixture` is an R package and Shiny app that enables users to visualise admixture as pie charts on a projected map. It also allows users to visualise admixture as traditional structure barplots or facet barplots.

## Installation

`mapmixture` requires [R](https://www.r-project.org/) (>= 4.2) to be installed on your system. Click [here](https://cran.r-project.org/bin/windows/base/) to download the latest version of R for Windows.

Install the latest stable release from CRAN:
``` r
install.packages("mapmixture")
```

Install the latest development version from GitHub:

``` r
# install.packages("devtools")
devtools::install_github("Tom-Jenkins/mapmixture")
```

## Reference
```r
mapmixture()         # main function
structure_plot()     # plot traditional structure or facet barplot
scatter_plot()       # plot PCA or DAPC results
launch_mapmixture()  # launch mapmixture Shiny app
```

## Citation
Jenkins TL (2024). [mapmixture]{.smallcaps}: an R package and web app for spatial visualisation of admixture and population structure. *Molecular Ecology Resources*, 24: e13943. DOI: 10.1111/1755-0998.13943.

## Examples

#### Basic usage of mapmixture

<details markdown="1">
<summary>Code</summary>
```{r dpi=300}
# Load package
library(mapmixture)

# Read in admixture file format 1
file <- system.file("extdata", "admixture1.csv", package = "mapmixture")
admixture1 <- read.csv(file)

# Read in coordinates file
file <- system.file("extdata", "coordinates.csv", package = "mapmixture")
coordinates <- read.csv(file)

# Run mapmixture
map1 <- mapmixture(admixture1, coordinates, crs = 3035)
# map1
```
</details>
```{r dpi=300, echo=FALSE}
map1
```


#### Customised usage of mapmixture with a high resolution map

<details markdown="1">
<summary>Code</summary>
```{r dpi=300}
# Load packages
library(mapmixture)
library(rnaturalearthhires)

# Install rnaturalearthhires package using:
# install.packages("rnaturalearthhires", repos = "https://ropensci.r-universe.dev", type = "source")

# Read in admixture file format 1
file <- system.file("extdata", "admixture1.csv", package = "mapmixture")
admixture1 <- read.csv(file)

# Read in coordinates file
file <- system.file("extdata", "coordinates.csv", package = "mapmixture")
coordinates <- read.csv(file)

# Run mapmixture
map2 <- mapmixture(
  admixture_df = admixture1,
  coords_df = coordinates,
  cluster_cols = c("#f1a340","#998ec3"),
  cluster_names = c("Group A","Group B"),
  crs = 3035,
  basemap = rnaturalearthhires::countries10[, c("geometry")],
  boundary = c(xmin=-15, xmax=16, ymin=40, ymax=62),
  pie_size = 1,
  pie_border = 0.3,
  pie_border_col = "white",
  pie_opacity = 1,
  land_colour = "#d9d9d9",
  sea_colour = "#deebf7",
  expand = TRUE,
  arrow = TRUE,
  arrow_size = 1.5,
  arrow_position = "bl",
  scalebar = TRUE,
  scalebar_size = 1.5,
  scalebar_position = "tl",
  plot_title = "Admixture Map",
  plot_title_size = 12,
  axis_title_size = 10,
  axis_text_size = 8
)
# map2
```
</details>
```{r dpi=300, echo=FALSE}
map2
```

#### Admixture map with single coloured circles

<details markdown="1">
<summary>Code</summary>
```{r dpi=300}
# Load package
library(mapmixture)

# Read in admixture file format 3
file <- system.file("extdata", "admixture3.csv", package = "mapmixture")
admixture3 <- read.csv(file)

# Read in coordinates file
file <- system.file("extdata", "coordinates.csv", package = "mapmixture")
coordinates <- read.csv(file)

# Run mapmixture
map3 <- mapmixture(admixture3, coordinates, crs = 3035)
# map3
```
</details>
```{r dpi=300, echo=FALSE}
map3
```


#### Add additional geoms or theme options to mapmixture ggplot object

<details markdown="1">
<summary>Code</summary>
```{r dpi=300}
# Load packages
library(mapmixture)
library(ggplot2)

# Read in admixture file format 1
file <- system.file("extdata", "admixture1.csv", package = "mapmixture")
admixture1 <- read.csv(file)

# Read in coordinates file
file <- system.file("extdata", "coordinates.csv", package = "mapmixture")
coordinates <- read.csv(file)

# Run mapmixture
map4 <- mapmixture(
  admixture_df = admixture1,
  coords_df = coordinates,
  cluster_cols = c("#f1a340","#998ec3"),
  cluster_names = c("Ancestry 1","Ancestry 2"),
  crs = 4326,
  boundary = c(xmin=-15, xmax=16, ymin=40, ymax=62),
  pie_size = 1,
)+
  # Add additional label to the map
  annotate("label",
    x = -10,
    y = 46.5,
    label = "Atlantic Ocean",
    size = 3,
  )+
  # Add additional text to the map
  annotate("text",
    x = 2.5,
    y = 57,
    label = "North Sea",
    size = 3,
  )+
  # Adjust ggplot theme options
  theme(
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8),
  )+
  # Adjust the size of the legend keys
  guides(fill = guide_legend(override.aes = list(size = 5, alpha = 1)))
# map4
```
</details>
```{r dpi=300, echo=FALSE}
map4
```


#### Combine admixture map and barplot ggplot objects into a single figure

<details markdown="1">
<summary>Code</summary>
```{r dpi=300}
# Load packages
library(mapmixture)
library(ggplot2)
library(gridExtra)

# Read in admixture file format 1
file <- system.file("extdata", "admixture1.csv", package = "mapmixture")
admixture1 <- read.csv(file)

# Read in coordinates file
file <- system.file("extdata", "coordinates.csv", package = "mapmixture")
coordinates <- read.csv(file)

# Run mapmixture
map5 <- mapmixture(
  admixture_df = admixture1,
  coords_df = coordinates,
  cluster_cols = c("#f1a340","#998ec3"),
  cluster_names = c("Ancestry 1","Ancestry 2"),
  crs = 4326,
  boundary = c(xmin=-20, xmax=20, ymin=40, ymax=62),
  pie_size = 1.3,
)+
  # Adjust theme options
  theme(
    legend.position = "top",
    plot.margin = margin(l = 10, r = 10),
  )+
  # Adjust the size of the legend keys
  guides(fill = guide_legend(override.aes = list(size = 5, alpha = 1)))

# Traditional structure barplot
structure_barplot <- structure_plot(
  admixture_df = admixture1,
  type = "structure",
  cluster_cols = c("#f1a340","#998ec3"),
  site_dividers = TRUE,
  divider_width = 0.4,
  site_order = c(
    "Vigo","Ile de Re","Isles of Scilly","Mullet Peninsula",
    "Shetland","Cromer","Helgoland","Flodevigen","Lysekil","Bergen"
  ),
  labels = "site",
  flip_axis = FALSE,
  site_ticks_size = -0.05,
  site_labels_y = -0.35,
  site_labels_size = 2.2
)+
  # Adjust theme options
  theme(
    axis.title.y = element_text(size = 8, hjust = 1),
    axis.text.y = element_text(size = 5),
  )

# Arrange plots
# grid.arrange(map5, structure_barplot, nrow = 2, heights = c(4,1))
```
</details>
```{r dpi=300, echo=FALSE, }
grid.arrange(map5, structure_barplot, nrow = 2, heights = c(4,1))
```

<details markdown="1">
<summary>Code</summary>
```{r dpi=300}
# Load packages
library(mapmixture)
library(ggplot2)
library(gridExtra)

# Read in admixture file format 1
file <- system.file("extdata", "admixture1.csv", package = "mapmixture")
admixture1 <- read.csv(file)

# Read in coordinates file
file <- system.file("extdata", "coordinates.csv", package = "mapmixture")
coordinates <- read.csv(file)

# Run mapmixture
map6 <- mapmixture(
  admixture_df = admixture1,
  coords_df = coordinates,
  cluster_cols = c("#f1a340","#998ec3"),
  cluster_names = c("Ancestry 1","Ancestry 2"),
  crs = 4326,
  boundary = c(xmin=-20, xmax=20, ymin=40, ymax=62),
  pie_size = 1.3,
)+
  # Adjust theme options
  theme(
    legend.position = "top",
    plot.margin = margin(l = 10, r = 10),
  )+
  # Adjust the size of the legend keys
  guides(fill = guide_legend(override.aes = list(size = 5, alpha = 1)))

# Facet structure barplot
facet_barplot <- structure_plot(admixture1,
  type = "facet",
  cluster_cols = c("#f1a340","#998ec3"),
  facet_col = 2,
  ylabel = "Admixture proportions",
)+
  theme(
    axis.title.y = element_text(size = 10),
    axis.text.y = element_text(size = 5),
    strip.text = element_text(size = 6, vjust = 1, margin = margin(t=1.5, r=0, b=1.5, l=0)),
  )

# Arrange plots
# grid.arrange(map6, facet_barplot, ncol = 2, widths = c(3,2))
```
</details>
```{r dpi=300, echo=FALSE, }
grid.arrange(map6, facet_barplot, ncol = 2, widths = c(3,2))
```


#### Use a raster as the basemap 

The raster (TIFF) used in the example below was downloaded from Natural Earth [here](https://www.naturalearthdata.com/downloads/50m-raster-data/50m-natural-earth-1/). You need to install the [terra](https://github.com/rspatial/terra) package to use this feature. Currently, the `basemap` argument accepts a `SpatRaster` or a `sf` object.

<details markdown="1">
<summary>Code</summary>
```{r dpi=300}
# Load packages
library(mapmixture)
library(terra)

# Create SpatRaster object
earth <- terra::rast("../NE1_50M_SR_W/NE1_50M_SR_W.tif")

# Read in admixture file format 1
file <- system.file("extdata", "admixture1.csv", package = "mapmixture")
admixture1 <- read.csv(file)

# Read in coordinates file
file <- system.file("extdata", "coordinates.csv", package = "mapmixture")
coordinates <- read.csv(file)

# Run mapmixture
map7 <- mapmixture(admixture1, coordinates, crs = 3035, basemap = earth)
# map7
```
</details>
```{r dpi=300, echo=FALSE}
map7
```


#### Add pie charts to an existing map

The vector data (shapefile) used in the example below was downloaded from the Natural England Open Data Geoportal [here](https://naturalengland-defra.opendata.arcgis.com/datasets/Defra::marine-conservation-zones-england/explore?location=53.749917%2C-5.534585%2C6.27).

<details markdown="1">
<summary>Code</summary>
```{r dpi=300}
# Load packages
library(mapmixture)
library(rnaturalearthhires)
library(ggplot2)
library(dplyr)
library(sf)

# Read in admixture file format 1
file <- system.file("extdata", "admixture1.csv", package = "mapmixture")
admixture1 <- read.csv(file)

# Read in coordinates file
file <- system.file("extdata", "coordinates.csv", package = "mapmixture")
coordinates <- read.csv(file)

# Parameters
crs <- 3035
boundary <- c(xmin=-11, xmax=13, ymin=50, ymax=60) |> transform_bbox(bbox = _, crs)

# Read in world countries from Natural Earth and transform to CRS
world <- rnaturalearthhires::countries10[, c("geometry")]
world <- st_transform(world, crs = crs)

# Read in Marine Conservation Zones shapefile
# Extract polygons for Western Channel, Offshore Brighton and Swallow Sand
# Transform to CRS
mczs <- st_read("../Marine_Conservation_Zones_England/Marine_Conservation_Zones___Natural_England_and_JNCC.shp", quiet = TRUE) |>
  dplyr::filter(.data = _, MCZ_NAME %in% c("Western Channel", "Offshore Brighton", "Swallow Sand")) |>
  st_transform(x = _, crs = crs)

# Run mapmixture helper functions to prepare admixture and coordinates data
admixture_df <- standardise_data(admixture1, type = "admixture") |> transform_admix_data(data = _)
coords_df <- standardise_data(coordinates, type = "coordinates")
admix_coords <- merge_coords_data(coords_df, admixture_df) |> transform_df_coords(df = _, crs = crs)

# Plot map and add pie charts
map8 <- ggplot()+
  geom_sf(data = world, colour = "black", fill = "#d9d9d9", size = 0.1)+
  geom_sf(data = mczs, aes(fill = "MCZs"), linewidth = 0.3)+
  scale_fill_manual(values = c("yellow"))+
  coord_sf(
    xlim = c(boundary[["xmin"]], boundary[["xmax"]]),
    ylim = c(boundary[["ymin"]], boundary[["ymax"]])
  )+
  add_pie_charts(admix_coords,
    admix_columns = 4:ncol(admix_coords),
    lat_column = "lat",
    lon_column = "lon",
    pie_colours = c("green","blue"),
    border = 0.3,
    opacity = 1,
    pie_size = 0.8
  )+
  theme(
    legend.title = element_blank(),
  )
# map8
```
</details>
```{r dpi=300, echo=FALSE}
map8
```


#### Scatter plot of PCA or DAPC results from genotypes

<details markdown="1">
<summary>Code</summary>
```{r dpi=300}
# Load packages
library(mapmixture)
library(ggplot2)
library(adegenet)
library(RColorBrewer)
library(gridExtra)

# Load example genotypes
data("dapcIllus")
geno = dapcIllus$a

# Change population labels
popNames(geno) = c("Pop1","Pop2","Pop3","Pop4","Pop5","Pop6")

# Region names
region_names <- rep(c("Region1", "Region2"), each = 300)

# Define colour palette
cols = brewer.pal(nPop(geno), "RdYlBu")

# Perform PCA
pca1 = dudi.pca(geno, scannf = FALSE, nf = 3)

# Percent of genetic variance explained by each axis
percent = round(pca1$eig/sum(pca1$eig)*100, digits = 1)

# Scatter plot with centroids and segments
scatter1 <- scatter_plot(
  dataframe = pca1$li,
  group_ids = geno$pop,
  type = "points",
  axes = c(1,2),
  percent = percent,
  colours = cols,
  point_size = 2,
  point_type = 21,
  centroid_size = 2,
  stroke = 0.1,
  plot_title = "PCA coloured by group_ids"
)+
  theme(
    legend.position = "none",
    axis.title = element_text(size = 8),
    axis.text = element_text(size = 6),
    plot.title = element_text(size = 10),
  )

# Same as scatter1 but no segments and axis 1 and 3 are shown
scatter2 <- scatter_plot(
  dataframe = pca1$li,
  group_ids = geno$pop,
  type = "points",
  axes = c(1,3),
  percent = percent,
  colours = cols,
  point_size = 2,
  point_type = 21,
  centroids = TRUE,
  centroid_size = 2,
  segments = FALSE,
  stroke = 0.1,
  plot_title = "PCA no segments and axis 1 and 3 shown"
)+
  theme(
    legend.position = "none",
    axis.title = element_text(size = 8),
    axis.text = element_text(size = 6),
    plot.title = element_text(size = 10),
  )

# Same as scatter1 but coloured by region
scatter3 <- scatter_plot(
  dataframe = pca1$li,
  group_ids = geno$pop,
  other_group = region_names,
  type = "points",
  axes = c(1,2),
  percent = percent,
  colours = cols,
  point_size = 2,
  point_type = 21,
  centroid_size = 2,
  stroke = 0.1,
  plot_title = "PCA coloured by other_group"
)+
  theme(
    legend.position = "none",
    axis.title = element_text(size = 8),
    axis.text = element_text(size = 6),
    plot.title = element_text(size = 10),
  )

# Scatter plot with labels instead of points
scatter4 <- scatter_plot(
  dataframe = pca1$li,
  group_ids = geno$pop,
  type = "labels",
  labels = rownames(pca1$li),
  colours = cols,
  size = 2,
  label.size = 0.10,
  label.padding = unit(0.10, "lines"),
  plot_title = "PCA using labels instead of points"
)+
  theme(
    legend.position = "none",
    axis.title = element_text(size = 8),
    axis.text = element_text(size = 6),
    plot.title = element_text(size = 10),
  )

# Arrange plots
# grid.arrange(scatter1, scatter2, scatter3, scatter4)
```
</details>
```{r dpi=300, echo=FALSE}
grid.arrange(scatter1, scatter2, scatter3, scatter4)
```


## Launch interactive Shiny app

```{r eval=FALSE}
# Load package
library(mapmixture)

# Launch Shiny app
launch_mapmixture()

# Tested with the following package versions:
# shiny v1.8.0 (important)
# shinyFeedback v0.4.0
# shinyjs v2.1.0
# shinyWidgets 0.8.4
# bslib 0.7.0
# colourpicker 1.3.0
# htmltools v0.5.8.1
# waiter 0.2.5
```

## Link to online Shiny app
https://tomjenkins.shinyapps.io/mapmixture/

## Format

```{r}
# Load package
library(mapmixture)

# Admixture Format 1
file <- system.file("extdata", "admixture1.csv", package = "mapmixture")
admixture1 <- read.csv(file)
head(admixture1)

# Admixture Format 2
file <- system.file("extdata", "admixture2.csv", package = "mapmixture")
admixture2 <- read.csv(file)
admixture2

# Admixture Format 3
file <- system.file("extdata", "admixture3.csv", package = "mapmixture")
admixture3 <- read.csv(file)
admixture3

# Coordinates
file <- system.file("extdata", "coordinates.csv", package = "mapmixture")
coordinates <- read.csv(file)
coordinates
```

## Related software

- [pophelper](https://github.com/royfrancis/pophelper)
- [pophelperShiny](https://github.com/royfrancis/pophelperShiny)