R Notebook

#Load excel sheet
#excel sheet should contain longitue and latitude
setwd("~/Desktop/Genomics Paper/4.21")
library(readxl)


tissue <- read_excel("4.21.tissue.scale.w.genomics.xlsx")

#PACKAGES REQUIRED

#if (!require(tidyverse)) install.packages('tidyverse')
#if (!require(sf)) install.packages('sf')
#if (!require(leaflet)) install.packages('leaflet')
#if (!require(mapview)) install.packages('mapview')

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.0     ✓ dplyr   1.0.5
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(sf)

## Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1

library(leaflet)
library(mapview)

#load whole map using st_as_sf function
whole.map = st_as_sf(tissue, coords = c("Longitude", "Latitude"),  crs=4326)


#class(whole.map)

#should be:
#[1] "sf"         "tbl_df"     "tbl"        "data.frame"

#quick view of map with no color 'filter' of 'Species' column on excel sheet
#mapview(whole.map, zcol = "Species")


#this allows visualization on the desired column name on the map, but no filtering 

#Subdivide map based on genetic species

#Montane vole~montanus 
Genetic.Microtus.montanus <- whole.map  %>% 
  filter(GENETIC.GROUP == "Montane")

#Woodland vole~pinetorium
Genetic.Microtus.pinetorum <- whole.map  %>% 
  filter(GENETIC.GROUP == "Woodland")

#Meadow vole clusters ~ Pennsylvanicus 
Genetic.Meadow.cluster.1  <- whole.map  %>% 
  filter(GENETIC.GROUP == "Meadow cluster 1")

Genetic.Meadow.cluster.2  <- whole.map  %>% 
  filter(GENETIC.GROUP == "Meadow cluster 2")

Genetic.Meadow.cluster.1.2.Hybrid  <- whole.map  %>% 
  filter(GENETIC.GROUP == "Meadow cluster 1.2.Hybrid")

#unknown group
Unknown <- whole.map  %>% 
  filter(GENETIC.GROUP == "Unknown")

#Map out genetic species 

genetic.map= 
mapview(Genetic.Microtus.montanus, cex =16, col.regions = "purple") +  mapview(Genetic.Microtus.pinetorum, cex =16, col.regions = "red") + mapview(Unknown, cex =16, col.regions = "black")+ mapview(Genetic.Meadow.cluster.1, cex =16, col.regions = "olivedrab1")+ mapview(Genetic.Meadow.cluster.1.2.Hybrid, cex =16, col.regions = "seagreen")+ mapview(Genetic.Meadow.cluster.2, cex =16, col.regions = "sienna1")

genetic.map