load packages
#load packages
library(vcfR)
##
## ***** *** vcfR *** *****
## This is vcfR 1.12.0
## browseVignettes('vcfR') # Documentation
## citation('vcfR') # Citation
## ***** ***** ***** *****
library(ggplot2)
library(adegenet)
## Loading required package: ade4
##
## /// adegenet 2.1.5 is loaded ////////////
##
## > overview: '?adegenet'
## > tutorials/doc/questions: 'adegenetWeb()'
## > bug reports/feature requests: adegenetIssues()
library(SNPfiltR)
## This is SNPfiltR v.1.0.0
##
## Detailed usage information is available at: devonderaad.github.io/SNPfiltR/
##
## If you use SNPfiltR in your published work, please cite the following papers:
##
## DeRaad, D.A. (2022), SNPfiltR: an R package for interactive and reproducible SNP filtering. Molecular Ecology Resources, 00, 1-15. http://doi.org/10.1111/1755-0998.13618
##
## Knaus, Brian J., and Niklaus J. Grunwald. 2017. VCFR: a package to manipulate and visualize variant call format data in R. Molecular Ecology Resources, 17.1:44-53. http://doi.org/10.1111/1755-0998.12549
library(StAMPP)
## Loading required package: pegas
## Loading required package: ape
## Registered S3 method overwritten by 'pegas':
## method from
## print.amova ade4
##
## Attaching package: 'pegas'
## The following object is masked from 'package:ape':
##
## mst
## The following object is masked from 'package:ade4':
##
## amova
## The following objects are masked from 'package:vcfR':
##
## getINFO, write.vcf
library(viridis)
## Loading required package: viridisLite
library(rgeoboundaries)
library(elevatr)
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following objects are masked from 'package:ape':
##
## rotate, zoom
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.4.0, PROJ 8.1.1; sf_use_s2() is TRUE
library(ggpubr)
##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:raster':
##
## rotate
## The following object is masked from 'package:ape':
##
## rotate
read in data
#read in your sample info file
pops<-read.csv("~/Desktop/phil.dicaeum/dicaeum.retained.sampling.csv")
#it should look roughly like this:
head(pops)
## ID Institition Tissue Species Taxa State
## 1 D_hypoleucum_18070 KU 18070 hypoleucum Luzon Camarines Norte
## 2 D_hypoleucum_18191 KU 18191 hypoleucum Zamboanga
## 3 D_hypoleucum_14037 KU 14037 hypoleucum Mindanao Dinagat
## 4 D_hypoleucum_14079 KU 14079 hypoleucum Mindanao Dinagat
## 5 D_hypoleucum_14065 KU 14065 hypoleucum Mindanao Dinagat
## 6 D_hypoleucum_14120 KU 14120 hypoleucum Mindanao Eastern Samar
## Specificlocality
## 1 Luzon Island; Municipality Labo; Barangay Tulay na Lupa; Mt. labo
## 2 Mindanao Island; City of Zamboanga; Brgy. Pasonanca; Pasonanca Nat. Park
## 3 Dinagat Island; Muncipality of Loreto; Brgy. Santiago
## 4 Dinagat Island; Muncipality of Loreto; Brgy. Santiago
## 5 Dinagat Island; Muncipality of Loreto; Brgy. Santiago
## 6 Samar Island; Muncipality of Taft; Brgy. San Rafael
## Latitude Longitude Altitude
## 1 14.039000 122.7860 580
## 2 7.018483 122.0273 750
## 3 10.343680 125.6181 200
## 4 10.343680 125.6181 200
## 5 10.343680 125.6181 190
## 6 11.828000 125.2760 200
#remove outgroup
pops<-pops[c(1:54,57:60),]
table(pops$Latitude)
##
## 6.109 7.018483 7.042 7.1078 7.164 7.848 8 8.473
## 5 3 3 2 2 5 3 1
## 8.721 8.7325 9.075 9.7056 10.34368 10.405 11.828 13.7644
## 1 1 5 1 4 1 3 1
## 14.039 15.64362 15.6851 17 18.219 18.438
## 2 1 3 2 5 4
#make dataframe with only the 8 unique sampling locs
unq<-pops[!duplicated(pops$Latitude),]
#download elevation data for the Philippines
phil_bound <- geoboundaries("Philippines")
elevation_data <- get_elev_raster(locations = phil_bound, z = 6, clip = "locations")
## Registered S3 method overwritten by 'httr':
## method from
## print.cache_info hoardr
## Mosaicing & Projecting
## Clipping DEM to locations
## Note: Elevation units are in meters.
elevation_data <- as.data.frame(elevation_data, xy = TRUE)
colnames(elevation_data)[3] <- "elevation"
# remove rows of data frame with one or more NA's,using complete.cases
elevation_data <- elevation_data[complete.cases(elevation_data), ]
# remove rows of data frame where elevation is lower than 0, to distinguish land from sea
elevation_data <- elevation_data[elevation_data$elevation >0,]
#set how big you want boundaries around points to be
g=1
#make subset plotting area corresponding to our sampling
#subset elevation data based on latitude and longitude limits
elevation_data<-elevation_data[elevation_data$x >= min(unq$Longitude)-g & elevation_data$x <= max(unq$Longitude)+g &
elevation_data$y >= min(unq$Latitude)-g & elevation_data$y <= max(unq$Latitude)+g,]
#subset map outline based on limits
phil_bound <- st_crop(phil_bound, xmin = min(unq$Longitude)-g, xmax = max(unq$Longitude)+g,
ymin = min(unq$Latitude)-g, ymax = max(unq$Latitude)+g)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
#plot the map with ggplot
ggplot() +
geom_raster(data = elevation_data, aes(x = x, y = y, fill = elevation)) +
#geom_sf(data = phil_bound, color = "black", fill = NA, cex=.3) +
geom_point(data = unq, aes(x = Longitude, y = Latitude, color=Taxa), size=3, alpha =1, show.legend=FALSE) +
scale_color_manual(values=c("black","grey","white"))+
geom_point(data = unq, aes(x = Longitude, y = Latitude), shape = 21,size=3, colour = "black")+
#geom_text(data = unq, aes(x = long, y = lat, label = sample.loc), size = 6, color="white")+
#scale_fill_gradient(low = "white", high = "black")+
scale_fill_gradientn(colours = terrain.colors(50))+
labs(x = "Longitude", y = "Latitude", fill = "Elevation\n(meters)")+
theme_classic()+
theme(legend.position = c(0.75, 0.75), legend.justification = c(0.01, 0.01),
legend.background = element_blank())
## Warning: Raster pixels are placed at uneven vertical intervals and will be
## shifted. Consider using geom_tile() instead.
make dots sized according to sample size
#add sample size for each point into the dataframe
unq$sample.size<-as.vector(table(pops$Latitude)[order(match(names(table(pops$Latitude)),unq$Latitude))])
#plot the map with ggplot
ggplot() +
geom_raster(data = elevation_data, aes(x = x, y = y, fill = elevation)) +
#geom_sf(data = phil_bound, color = "black", fill = NA, cex=.3) +
geom_point(data = unq, aes(x = Longitude, y = Latitude, color=Taxa, size=sample.size), alpha=1, show.legend=FALSE) +
scale_color_manual(values=c("black","grey","white"))+
geom_point(data = unq, aes(x = Longitude, y = Latitude, size=sample.size), shape=21, colour="black")+
#geom_text(data = unq, aes(x = long, y = lat, label = sample.loc), size = 6, color="white")+
#scale_fill_gradient(low = "white", high = "black")+
scale_size(range = c(3,6))+
scale_fill_gradientn(colours = terrain.colors(50))+
labs(x = "Longitude", y = "Latitude", fill = "Elevation\n(meters)", size="sample size")+
theme_classic()+
theme(legend.position = c(0.77, 0.6), legend.justification = c(0.01, 0.01),
legend.background = element_blank())
## Warning: Raster pixels are placed at uneven vertical intervals and will be
## shifted. Consider using geom_tile() instead.
library(ggview)
#ggview(units="in", width=5.3, height=11)
#ggsave(file="~/Desktop/phil.dicaeum/dicaeum.sampling.map.pdf", units="in",width=5.3,height=11) #saves g
Add species range outline in photoshop
#see final product
knitr::include_graphics("/Users/devder/Desktop/phil.dicaeum/dicaeum.fixed.sampling.png")
