Plotting smoothed surface diagrams of allele frequencies obtained from AlleleShift on a baseline map via ggmap
Roeland Kindt
23/7/2021
- First version: 23-JUL-2021
- Updated: 18-JAN-2022
1 Packages needed
library(AlleleShift)
library(BiodiversityR) # also loads vegan
library(ggmap)
library(dplyr)
library(ggrepel)
library(patchwork)2 Introduction
The AlleleShift package (Kindt 2021) predicts changes in allele frequencies from baseline to changed climates. The package includes the shift.surf.ggplot function that generates and plots smoothed regression surfaces, internally via vegan::ordisurf and BiodiversityR::ordisurfgrid.long (see Kindt 2021 Figure 5).
Here I show how smoothed regression surfaces can be plotted on a baseline map obtained via the ggmap package.
3 Step 1: Obtain the geospatial data
3.1 Load the locations of the populations
data(Poptri.loc)
head(Poptri.loc)## Pop Population Lat Long orig.index Latitude.index
## P1 P1 Chilliwack 49.39234 -121.7522 1 10
## P2 P2 Columbia 46.10763 -122.9206 2 3
## P3 P3 Dean 52.81705 -126.7579 3 18
## P4 P4 Harrison 50.27599 -122.8218 4 15
## P5 P5 Homathko 51.12257 -124.8980 5 16
## P6 P6 Kitimat 54.14684 -128.5952 6 193.2 Get a stamen map
Changing the zoom level will change the resolution and the size of the final graph.
bbox_study <- c(left=min(Poptri.loc$Long) - 0.9,
bottom=min(Poptri.loc$Lat) - 0.9,
right=max(Poptri.loc$Long) + 0.9,
top=max(Poptri.loc$Lat) + 0.9)
stamen_study <- get_stamenmap(bbox_study, zoom=6, maptype="terrain-background")## Source : http://tile.stamen.com/terrain-background/6/8/20.png## Source : http://tile.stamen.com/terrain-background/6/9/20.png## Source : http://tile.stamen.com/terrain-background/6/10/20.png## Source : http://tile.stamen.com/terrain-background/6/8/21.png## Source : http://tile.stamen.com/terrain-background/6/9/21.png## Source : http://tile.stamen.com/terrain-background/6/10/21.png## Source : http://tile.stamen.com/terrain-background/6/8/22.png## Source : http://tile.stamen.com/terrain-background/6/9/22.png## Source : http://tile.stamen.com/terrain-background/6/10/22.png## Source : http://tile.stamen.com/terrain-background/6/8/23.png## Source : http://tile.stamen.com/terrain-background/6/9/23.png## Source : http://tile.stamen.com/terrain-background/6/10/23.png## Source : http://tile.stamen.com/terrain-background/6/8/24.png## Source : http://tile.stamen.com/terrain-background/6/9/24.png## Source : http://tile.stamen.com/terrain-background/6/10/24.png3.3 Check the map
map_study <- ggmap(stamen_study) +
geom_point(data=Poptri.loc,
aes(x=Long, y=Lat),
colour="red")
map_study
4 Step 2: Obtain data on changes in allele frequencies
4.1 Data on frequencies for each population
AlleleShift includes baseline and future allele frequency data that were generated by prediction methods available in the package.
data(Poptri.freq.baseline)
data(Poptri.freq.future)
head(Poptri.freq.baseline)## Pop Pop.label Pop.index N Allele Allele.freq A B
## 1 Chilliwack 1 1 76 X01_10838495.1 0.10526316 8 68
## 2 Columbia 2 2 360 X01_10838495.1 0.23055556 83 277
## 3 Dean 3 3 14 X01_10838495.1 0.21428571 3 11
## 4 Harrison 4 4 44 X01_10838495.1 0.15909091 7 37
## 5 Homathko 5 5 18 X01_10838495.1 0.05555556 1 17
## 6 Kitimat 6 6 16 X01_10838495.1 0.25000000 4 12
## Ap Bp N.e1 Freq.e1 Freq.e2 LCL UCL increasing
## 1 8.950722 67.04928 76 0.1177727 0.1144837 0.09803780 0.1309296 TRUE
## 2 88.915206 271.08479 360 0.2469867 0.2079033 0.18915823 0.2266484 FALSE
## 3 1.816999 12.18300 14 0.1297857 0.1290953 0.11110733 0.1470833 FALSE
## 4 9.875440 34.12456 44 0.2244418 0.1898990 0.17229199 0.2075060 TRUE
## 5 3.289572 14.71043 18 0.1827540 0.2008387 0.17374709 0.2279304 TRUE
## 6 1.694003 14.30600 16 0.1058752 0.1050042 0.09154483 0.1184635 FALSE
## Long Lat Latitude.index
## 1 -121.7522 49.39234 10
## 2 -122.9206 46.10763 3
## 3 -126.7579 52.81705 18
## 4 -122.8218 50.27599 15
## 5 -124.8980 51.12257 16
## 6 -128.5952 54.14684 194.2 Fit regression surfaces via the ordisurf function from ordisurf
Regression surfaces are fitted separately for baseline or future frequencies for selected alleles.
First get the names of the different alleles.
Allele.names <- unique(as.character(Poptri.freq.baseline$Allele))
Allele.names## [1] "X01_10838495.1" "X01_16628872.1" "X01_23799134.1" "X01_38900650.1"
## [5] "X01_41284594.1"Obtain data for the regression surfaces via vegan::ordisurf and BiodiversityR::ordisurfgrid.long. A similar procedure is used internally in the shift.surf.ggplot function of AlleleShift.
# Baseline frequencies for allele #5
Poptri.B1 <- Poptri.freq.baseline[Poptri.freq.baseline$Allele == Allele.names[1], ]
ordiplot.B1 <- ordiplot(Poptri.B1[, c("Long", "Lat")])## species scores not available
surf.B1 <- ordisurfgrid.long(ordisurf(ordiplot.B1,
y=Poptri.B1$Allele.freq))
surf.B1 <- surf.B1[complete.cases(surf.B1), ]
head(surf.B1)## x y z
## 31 -128.5952 54.14684 0.2203447
## 62 -128.3188 54.14684 0.2160397
## 91 -128.0424 53.14904 0.1937580
## 92 -128.0424 53.64794 0.2033720
## 121 -127.7661 52.65014 0.1779051
## 122 -127.7661 53.14904 0.1886474# Baseline frequencies for allele #3
Poptri.B3 <- Poptri.freq.baseline[Poptri.freq.baseline$Allele == Allele.names[3], ]
ordiplot.B3 <- ordiplot(Poptri.B3[, c("Long", "Lat")])## species scores not available
surf.B3 <- ordisurfgrid.long(ordisurf(ordiplot.B3,
y=Poptri.B3$Allele.freq))
surf.B3 <- surf.B3[complete.cases(surf.B3), ]
# Future frequencies for allele #1
Poptri.F1 <- Poptri.freq.future[Poptri.freq.future$Allele == Allele.names[1], ]
ordiplot.F1 <- ordiplot(Poptri.F1[, c("Long", "Lat")])## species scores not available
surf.F1 <- ordisurfgrid.long(ordisurf(ordiplot.F1,
y=Poptri.F1$Freq.e2))
surf.F1 <- surf.F1[complete.cases(surf.F1), ]
# Future frequencies for allele #3
Poptri.F3 <- Poptri.freq.future[Poptri.freq.future$Allele == Allele.names[3], ]
ordiplot.F3 <- ordiplot(Poptri.F3[, c("Long", "Lat")])## species scores not available
surf.F3 <- ordisurfgrid.long(ordisurf(ordiplot.F3,
y=Poptri.F3$Freq.e2))
surf.F3 <- surf.F3[complete.cases(surf.F3), ]5 Step 3: Use geom_raster to plot information on the regression surface
map_study.B1 <- ggmap(stamen_study) +
geom_raster(data=surf.B1,
aes(x=x, y=y, fill=z),
alpha=0.8) +
geom_point(data=Poptri.loc,
aes(x=Long, y=Lat),
colour="black", size=4, shape=21, fill="grey50") +
geom_label_repel(data=Poptri.loc,
aes(x=Long, y=Lat, label=Pop),
size=2, alpha=0.7) +
scale_fill_viridis_b() +
labs(fill="Frequency") +
theme(axis.title = element_blank(),
legend.position=c(0.2, 0.2)) +
coord_fixed() ## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.B1
# Include information on direction of changes
Poptri.loc2 <- dplyr::left_join(Poptri.loc,
Poptri.F1[, c("Pop", "increasing")],
by=c("Population" = "Pop"))
map_study.F1 <- ggmap(stamen_study) +
geom_raster(data=surf.F1,
aes(x=x, y=y, fill=z),
alpha=0.8) +
geom_point(data=Poptri.loc2,
aes(x=Long, y=Lat),
shape=1, size=4, colour="black") +
geom_point(data=Poptri.loc2,
aes(x=Long, y=Lat, colour=increasing),
shape=16, size=4) +
scale_fill_viridis_b() +
scale_colour_manual(values=c("firebrick3", "chartreuse4"),
labels=c("decrease", "increase")) +
labs(fill="Frequency", colour="Change") +
theme(axis.title = element_blank(),
legend.position=c(0.2, 0.3)) +
coord_fixed() ## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.F1
map_study.B3 <- ggmap(stamen_study) +
geom_raster(data=surf.B3,
aes(x=x, y=y, fill=z),
alpha=0.8) +
geom_point(data=Poptri.loc,
aes(x=Long, y=Lat),
colour="black", size=4, shape=21, fill="grey50") +
scale_fill_viridis_b() +
labs(fill="Frequency") +
theme(axis.title = element_blank(),
legend.position=c(0.2, 0.2)) +
coord_fixed() ## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.B3
# Include information on direction of changes
Poptri.loc2 <- dplyr::left_join(Poptri.loc,
Poptri.F3[, c("Pop", "increasing")],
by=c("Population" = "Pop"))
map_study.F3 <- ggmap(stamen_study) +
geom_raster(data=surf.F3,
aes(x=x, y=y, fill=z),
alpha=0.8) +
geom_point(data=Poptri.loc2,
aes(x=Long, y=Lat),
size=4, shape=1, colour="black") +
geom_point(data=Poptri.loc2,
aes(x=Long, y=Lat, colour=increasing),
size=4, shape=16) +
scale_fill_viridis_b() +
scale_colour_manual(values=c("firebrick3", "chartreuse4"),
labels=c("decrease", "increase")) +
labs(fill="Frequency", colour="Change") +
theme(axis.title = element_blank(),
legend.position=c(0.2, 0.3)) +
coord_fixed() ## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.F3
6 Step 4: Generate a high resolution figure via patchwork and ragg
map_study <- (map_study.B1 + map_study.F1) /
(map_study.B3 + map_study.F3) +
plot_annotation(tag_levels="A") &
theme(plot.tag = element_text(size = 20))
# plot(map_study)
filesave <- paste0(getwd(), "//Allele frequency map.png")
ragg::agg_png(filesave, width = 180, height = 250, units = "mm", res = 300, scaling=0.9)
plot(map_study)
dev.off()## png
## 2Show the figure that was created
knitr::include_graphics(filesave)
7 UPDATE 1: Using geom_contour_filled
With the coordinate reference system of the baseline map, it is better to use gglot2’s geom_contour_filled or geom_contour (See also here).
7.1 Modify step 3
map_study.B1 <- ggmap(stamen_study) +
geom_contour_filled(data=surf.B1,
aes(x=x, y=y, z=z),
alpha=0.8) +
geom_point(data=Poptri.loc,
aes(x=Long, y=Lat),
colour="black", size=4, shape=21, fill="grey50") +
geom_label_repel(data=Poptri.loc,
aes(x=Long, y=Lat, label=Pop),
size=2, alpha=0.7) +
scale_fill_viridis_d() +
labs(fill="Frequency") +
theme(axis.title = element_blank(),
axis.text = element_text(size=7),
legend.position="right") +
coord_sf(crs=4326, expand=FALSE) ## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.B1
# Include information on direction of changes
Poptri.loc2 <- dplyr::left_join(Poptri.loc,
Poptri.F1[, c("Pop", "increasing")],
by=c("Population" = "Pop"))
map_study.F1 <- ggmap(stamen_study) +
geom_contour_filled(data=surf.F1,
aes(x=x, y=y, z=z),
alpha=0.8) +
geom_point(data=Poptri.loc2,
aes(x=Long, y=Lat),
shape=1, size=4, colour="black") +
geom_point(data=Poptri.loc2,
aes(x=Long, y=Lat, colour=increasing),
shape=16, size=4) +
scale_fill_viridis_d() +
scale_colour_manual(values=c("firebrick3", "chartreuse4"),
labels=c("decrease", "increase")) +
labs(fill="Frequency", colour="Change") +
theme(axis.title = element_blank(),
axis.text = element_text(size=7),
legend.position="right") +
coord_sf(crs=4326, expand=FALSE) ## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.F1
map_study.B3 <- ggmap(stamen_study) +
geom_contour_filled(data=surf.B3,
aes(x=x, y=y, z=z),
alpha=0.8) +
geom_point(data=Poptri.loc,
aes(x=Long, y=Lat),
colour="black", size=4, shape=21, fill="grey50") +
scale_fill_viridis_d() +
labs(fill="Frequency") +
theme(axis.title = element_blank(),
axis.text = element_text(size=7),
legend.position="right") +
coord_sf(crs=4326, expand=FALSE) ## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.B3
# Include information on direction of changes
Poptri.loc2 <- dplyr::left_join(Poptri.loc,
Poptri.F3[, c("Pop", "increasing")],
by=c("Population" = "Pop"))
map_study.F3 <- ggmap(stamen_study) +
geom_contour_filled(data=surf.F3,
aes(x=x, y=y, z=z),
alpha=0.8) +
geom_point(data=Poptri.loc2,
aes(x=Long, y=Lat),
size=4, shape=1, colour="black") +
geom_point(data=Poptri.loc2,
aes(x=Long, y=Lat, colour=increasing),
size=4, shape=16) +
scale_fill_viridis_d() +
scale_colour_manual(values=c("firebrick3", "chartreuse4"),
labels=c("decrease", "increase")) +
labs(fill="Frequency", colour="Change") +
theme(axis.title = element_blank(),
axis.text = element_text(size=7),
legend.position="right") +
coord_sf(crs=4326, expand=FALSE)## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.F3
7.2 Repeat step 4
map_study <- (map_study.B1 + map_study.F1) /
(map_study.B3 + map_study.F3) +
plot_annotation(tag_levels="A") &
theme(plot.tag = element_text(size = 20))
# plot(map_study)
filesave2 <- paste0(getwd(), "//Allele frequency map contour filled.png")
ragg::agg_png(filesave2, width = 180, height = 250, units = "mm", res = 300, scaling=0.9)
plot(map_study)
dev.off()## png
## 2Show the figure that was created
knitr::include_graphics(filesave2)
8 UPDATE 2: Using geom_contour
Instead ofgeom_contour_filled, now we use geom_contour.
8.1 Modify step 3
map_study.B1 <- ggmap(stamen_study) +
geom_contour(data=surf.B1,
aes(x=x, y=y, z=z, colour=after_stat(level)),
size=2, alpha=0.8) +
geom_point(data=Poptri.loc,
aes(x=Long, y=Lat),
colour="black", size=4, shape=21, fill="grey50") +
geom_label_repel(data=Poptri.loc,
aes(x=Long, y=Lat, label=Pop),
size=2, alpha=0.7) +
scale_colour_viridis_b() +
labs(colour="Frequency") +
theme(axis.title = element_blank(),
axis.text = element_text(size=7),
legend.position="right") +
coord_sf(crs=4326, expand=FALSE) ## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.B1
# Include information on direction of changes
Poptri.loc2 <- dplyr::left_join(Poptri.loc,
Poptri.F1[, c("Pop", "increasing")],
by=c("Population" = "Pop"))
map_study.F1 <- ggmap(stamen_study) +
geom_contour(data=surf.F1,
aes(x=x, y=y, z=z, colour=after_stat(level)),
alpha=0.8, size=2) +
geom_point(data=Poptri.loc2,
aes(x=Long, y=Lat, shape=increasing),
size=4, colour="black", fill="white") +
geom_point(data=subset(Poptri.loc2, increasing==FALSE),
aes(x=Long, y=Lat),
shape=25, size=4, colour="black", fill="firebrick3", show.legend=FALSE) +
geom_point(data=subset(Poptri.loc2, increasing==TRUE),
aes(x=Long, y=Lat),
shape=24, size=4, colour="black", fill="chartreuse4", show.legend=FALSE) +
scale_colour_viridis_b() +
scale_shape_manual(values=c(25, 24),
labels=c("decrease", "increase")) +
labs(colour="Frequency", shape="Change") +
theme(axis.title = element_blank(),
axis.text = element_text(size=7),
legend.position="right") +
coord_sf(crs=4326, expand=FALSE) ## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.F1
map_study.B3 <- ggmap(stamen_study) +
geom_contour(data=surf.B3,
aes(x=x, y=y, z=z, colour=after_stat(level)),
size=2, alpha=0.8) +
geom_point(data=Poptri.loc,
aes(x=Long, y=Lat),
colour="black", size=4, shape=21, fill="grey50") +
scale_colour_viridis_b() +
labs(colour="Frequency") +
theme(axis.title = element_blank(),
axis.text = element_text(size=7),
legend.position="right") +
coord_sf(crs=4326, expand=FALSE) ## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.B3
# Include information on direction of changes
Poptri.loc2 <- dplyr::left_join(Poptri.loc,
Poptri.F3[, c("Pop", "increasing")],
by=c("Population" = "Pop"))
map_study.F3 <- ggmap(stamen_study) +
geom_contour(data=surf.F3,
aes(x=x, y=y, z=z, colour=after_stat(level)),
alpha=0.8, size=2) +
geom_point(data=Poptri.loc2,
aes(x=Long, y=Lat, shape=increasing),
size=4, colour="black", fill="white") +
geom_point(data=subset(Poptri.loc2, increasing==FALSE),
aes(x=Long, y=Lat),
shape=25, size=4, colour="black", fill="firebrick3", show.legend=FALSE) +
geom_point(data=subset(Poptri.loc2, increasing==TRUE),
aes(x=Long, y=Lat),
shape=24, size=4, colour="black", fill="chartreuse4", show.legend=FALSE) +
scale_colour_viridis_b() +
scale_shape_manual(values=c(25, 24),
labels=c("decrease", "increase")) +
labs(colour="Frequency", shape="Change") +
theme(axis.title = element_blank(),
axis.text = element_text(size=7),
legend.position="right") +
coord_sf(crs=4326, expand=FALSE)## Coordinate system already present. Adding new coordinate system, which will replace the existing one.map_study.F3
8.2 Repeat step 4
map_study <- (map_study.B1 + map_study.F1) /
(map_study.B3 + map_study.F3) +
plot_annotation(tag_levels="A") &
theme(plot.tag = element_text(size = 20))
# plot(map_study)
filesave3 <- paste0(getwd(), "//Allele frequency map contour.png")
ragg::agg_png(filesave3, width = 180, height = 250, units = "mm", res = 300, scaling=0.9)
plot(map_study)
dev.off()## png
## 2Show the figure that was created
knitr::include_graphics(filesave3)
9 Session Information
sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United Kingdom.1252
## [2] LC_CTYPE=English_United Kingdom.1252
## [3] LC_MONETARY=English_United Kingdom.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United Kingdom.1252
##
## attached base packages:
## [1] tcltk stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] patchwork_1.1.1 ggrepel_0.8.2 dplyr_1.0.2
## [4] ggmap_3.0.0 ggplot2_3.3.3 BiodiversityR_2.14-1
## [7] vegan_2.5-6 lattice_0.20-41 permute_0.9-5
## [10] AlleleShift_1.0-2
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.1.7 Hmisc_4.4-0
## [4] systemfonts_0.3.2 plyr_1.8.6 igraph_1.2.6
## [7] sp_1.4-6 splines_4.0.2 digest_0.6.25
## [10] htmltools_0.5.1.1 gdata_2.18.0 relimp_1.0-5
## [13] magrittr_1.5 checkmate_2.0.0 cluster_2.1.0
## [16] openxlsx_4.1.5 tcltk2_1.2-11 gmodels_2.18.1
## [19] sandwich_2.5-1 prettyunits_1.1.1 jpeg_0.1-8.1
## [22] colorspace_1.4-1 mitools_2.4 textshaping_0.2.1
## [25] haven_2.3.1 xfun_0.15 crayon_1.3.4
## [28] lme4_1.1-23 survival_3.1-12 zoo_1.8-8
## [31] ape_5.4 glue_1.4.1 gtable_0.3.0
## [34] seqinr_4.2-4 car_3.0-8 abind_1.4-5
## [37] scales_1.1.1 DBI_1.1.0 Rcpp_1.0.7
## [40] isoband_0.2.1 viridisLite_0.3.0 xtable_1.8-4
## [43] progress_1.2.2 spData_0.3.8 htmlTable_2.0.0
## [46] units_0.6-7 foreign_0.8-80 spdep_1.1-5
## [49] Formula_1.2-3 survey_4.0 htmlwidgets_1.5.1
## [52] httr_1.4.2 RColorBrewer_1.1-2 acepack_1.4.1
## [55] ellipsis_0.3.1 farver_2.0.3 pkgconfig_2.0.3
## [58] nnet_7.3-14 deldir_0.2-3 labeling_0.3
## [61] tidyselect_1.1.0 rlang_0.4.8 reshape2_1.4.4
## [64] later_1.1.0.1 munsell_0.5.0 cellranger_1.1.0
## [67] tools_4.0.2 generics_0.1.0 ade4_1.7-16
## [70] evaluate_0.14 stringr_1.4.0 fastmap_1.0.1
## [73] ragg_0.4.0 yaml_2.2.1 knitr_1.28
## [76] zip_2.0.4 purrr_0.3.4 RgoogleMaps_1.4.5.3
## [79] nlme_3.1-148 mime_0.9 adegenet_2.1.3
## [82] compiler_4.0.2 rstudioapi_0.11 curl_4.3
## [85] png_0.1-7 e1071_1.7-3 tibble_3.0.1
## [88] statmod_1.4.34 stringi_1.4.6 forcats_0.5.0
## [91] Matrix_1.2-18 classInt_0.4-3 nloptr_1.2.2.1
## [94] vctrs_0.3.4 effects_4.1-4 RcmdrMisc_2.7-2
## [97] pillar_1.4.4 LearnBayes_2.15.1 lifecycle_0.2.0
## [100] data.table_1.12.8 bitops_1.0-6 raster_3.5-11
## [103] httpuv_1.5.4 R6_2.4.1 latticeExtra_0.6-29
## [106] promises_1.1.1 KernSmooth_2.23-17 gridExtra_2.3
## [109] rio_0.5.16 codetools_0.2-16 boot_1.3-25
## [112] MASS_7.3-51.6 gtools_3.8.2 rjson_0.2.20
## [115] withr_2.2.0 nortest_1.0-4 Rcmdr_2.7-2
## [118] mgcv_1.8-31 expm_0.999-5 parallel_4.0.2
## [121] hms_0.5.3 terra_1.4-22 grid_4.0.2
## [124] rpart_4.1-15 tidyr_1.1.2 coda_0.19-3
## [127] class_7.3-17 minqa_1.2.4 rmarkdown_2.3
## [130] carData_3.0-4 sf_0.9-6 shiny_1.4.0.2
## [133] base64enc_0.1-3