Almost every data set we will work on has a spatial context. Until relatively recently, this was largely ignored and at most given as a 2-dimensional projection of sampling locations on a black and white map in our manuscripts. Here we explore methodologies to provide more intuitive approaches for plotting the spatial location of our sampling sites.
First, we need to load in the proper libraries. Each of the libraries are needed for:
We will start by loading in the default data set for gstudio and extracting the coordinates of our populations.
# load in the libraries and grab the data
library(gstudio)
data(arapat)
# pull out the coordinates - Population is the default stratum
coords <- strata_coordinates(arapat)
coords$Stratum <- as.character( coords$Stratum )
coords[1:5,] Stratum Longitude Latitude
1 88 -114.2935 29.32541
11 9 -113.9449 29.01457
20 84 -113.6679 28.96651
29 175 -113.4897 28.72796
36 177 -113.9914 28.66056We can easily turn this coordinate dataset into an interactive map using the leaflet package.
library(leaflet)
coords |>
leaflet() |>
addTiles() |>
addMarkers( ~Longitude, ~Latitude, popup = ~Stratum)It is also possible to plot the data in a spatial context where the size of our markers represents some biological property, like allele frequencies or genetic diversity. Here, we’ll plot the locations with circle sizes proportional to the allele frequency of the LTRS locus or the effective allele size.
# Break out allele frequencies for LTRS by population
f <- frequencies(arapat,loci="LTRS",stratum="Population")
# lets only use the allele 01
f <- f[ f$Allele=="01",]
# here are what the data look like
f[1:10,] Stratum Locus Allele Frequency
1 101 LTRS 01 0.2777778
3 102 LTRS 01 0.3125000
5 12 LTRS 01 0.8000000
7 153 LTRS 01 0.1500000
9 156 LTRS 01 0.0000000
11 157 LTRS 01 0.0000000
13 159 LTRS 01 0.8888889
15 160 LTRS 01 0.8000000
17 161 LTRS 01 0.9500000
19 162 LTRS 01 0.9500000Now we merge these frequencies into our coordinates:
Stratum Longitude Latitude Frequency
1 101 -110.5744 27.90509 0.2777778
2 102 -109.1263 26.38014 0.3125000
3 12 -112.6655 27.18232 0.8000000
4 153 -110.4624 24.13389 0.1500000
5 156 -109.9890 24.04380 0.0000000We can also extract the effective allele size (Ae) and merge it:
# grab effect for LTRS locus
ae <- genetic_diversity(arapat, stratum = "Population", mode="Ae")
ae <- ae[ ae$Locus=="LTRS",]
# merge in the data
coords <- merge( coords, ae[,c(1,3)])
coords[1:5,] Stratum Longitude Latitude Frequency Ae
1 101 -110.5744 27.90509 0.2777778 1.670103
2 102 -109.1263 26.38014 0.3125000 1.753425
3 12 -112.6655 27.18232 0.8000000 1.470588
4 153 -110.4624 24.13389 0.1500000 1.342282
5 156 -109.9890 24.04380 0.0000000 1.000000We can now plot our coordinates interactively, with circle sizes scaled by the frequency of allele 01 or the effective allele size (Ae).
coords |>
leaflet() |>
addTiles() |>
addCircleMarkers( ~Longitude, ~Latitude, popup = ~Stratum, radius = ~Frequency * 15, color = "blue", fillOpacity = 0.6)You can change the underlying base map by passing on other Tiles to the construction. A list of available tile providers can be found here.
coords |>
leaflet() |>
addProviderTiles("Esri.WorldTopoMap") |>
addMarkers( ~Longitude, ~Latitude, popup=~Stratum )Or you can use dark-themed tiles and add custom-sized points scaled to represent effective allele size:
coords$Size <- log(coords$Ae) * 10 + 1
coords |>
leaflet() |>
addProviderTiles("CartoDB.DarkMatter") |>
addCircleMarkers( ~Longitude, ~Latitude, popup=~Stratum, radius=~Size, color="red", fillOpacity=0.75) In the information that pops out of the marker, you can embed a HTML table that contains richer detailed data, such as several diversity measures (\(H_O\), \(H_E\), and \(A_E\)).
x <- genetic_diversity(arapat,stratum = "Population",mode="Ho")
x <- x[ x$Locus == "LTRS",]
coords <- merge( coords, x[,c(1,3)])
x <- genetic_diversity(arapat,stratum = "Population",mode="He")
x <- x[ x$Locus == "LTRS",]
coords <- merge( coords, x[,c(1,3)])Let’s make these data into a table format in HTML:
coords$label <- paste("<h2>",coords$Stratum,"</h2><table><tr><th>Parameter</th><th>Value</th></tr><tr><td>H<sub>O</sub></td><td>", format(coords$Ho,digits = 3),"</td></tr><tr><td>H<sub>E</sub></td><td>",format(coords$He, digits=3), "</td></tr><tr><td>A<sub>E</sub></td><td>",format(coords$Ae,digits=3),"</td></tr></table>",sep="")
# plot it
coords |>
leaflet() |>
addTiles() |>
addCircleMarkers( ~Longitude, ~Latitude, popup=~label, radius=~Size, color="red", fillOpacity=0.75) A Sankey Diagram is a kind of flow diagram used to indicate allocation of one set of objects into two other sets. This is a complicated set of interactions to display. In the following example, I’ll use the populations from the A. attenuatus data set and then map each individual onto its STRUCTURE cluster and mtDNA clade. This plot needs a data.frame for nodes and one for edges created in a particular way.
library(networkD3)
pops <- as.character(unique(arapat$Population))
species <- as.character(unique(arapat$Species))
cluster <- as.character(unique(arapat$Cluster))
sankey.nodes <- data.frame( name = c( pops, species, cluster),
stringsAsFactors = FALSE)
# takes vector of values and finds out the number of times
# each combination is found in the pair.
get_sets <- function( source, target ) {
t <- table(source, target)
m <- which( t > 0, arr.ind = TRUE)
suppressWarnings( df <- data.frame( ind=m ) )
df$source <- rownames(t)[df$ind.source]
df$target <- colnames(t)[df$ind.target]
df$value <- sapply(seq(1,dim(df)[1]),
function(x) {
return(t[df$ind.source[x],df$ind.target[x]])
})
return( df[,3:5] )
}
# For finding the number of the name instead of the name
nameIndex <- function(vals,names){
ret <- sapply( vals, function(x) {which(x==names)-1}, USE.NAMES = FALSE)
return(ret)
}
sankey.edges <- rbind( get_sets(arapat$Species,arapat$Population),
get_sets(arapat$Population, arapat$Cluster) )
sankey.edges$targetNum <- nameIndex(sankey.edges$target, sankey.nodes$name)
sankey.edges$sourceNum <- nameIndex(sankey.edges$source, sankey.nodes$name)What it does produce is a very cool and dynamic plot. The colored boxes represent the groups–Species & STRUCTURE Cluster on the outside, population on the inside. Connecting lines indicate individuals who are in each of these partially overlapping groups. For example, if you grab the population Mat, you will see it has samples in it who are identified as either mtDNA type Peninsula or Cape as well as belonging to STRUCTURE groups CBP-C and SCBP-A. If you have difficult data such as these, this kind of plot may be very helpful.
sankeyNetwork( Links=sankey.edges, Nodes=sankey.nodes,
Source = "sourceNum", Target = "targetNum",
Value = "value", NodeID = "name",
width = 700, fontSize = 12, nodeWidth = 30)