rm(list = ls()) setwd("C:/Users/els419/Dropbox/Ongoing Projects/Cockles_2013_final/RELAXED COND ANALYSIS 2017/RDA") #PC setwd("~/Dropbox/Ongoing Projects/Cockles_2013_final/RELAXED COND ANALYSIS 2017/RDA") #MAC ### Download libraries library(vegan) library(reshape2) library(FactoMineR) library(SoDA) library(dplyr) library(adespatial) library(psych) library(adegenet) library(ggvegan) library(ggord) library(ggplot2) ######################################## ############# AEM ANALYSIS ############# ######################################## ###read nodes-to-edges-matrix### #nodes40 <- read.csv("connectivity 40km_June2019/settlement_radius_40km/edges-to-nodes40.csv", header = TRUE, row.names = 1) #nodes40.m = as.matrix(nodes) nodes <- read.csv("connectivity_newsimulations_Feb2019/edges-to-nodes2.csv", header = TRUE, row.names = 1) nodes.m = as.matrix(nodes) ###read larval dispersal as weights### weights <- read.csv("connectivity_newsimulations_Feb2019/edges_weights2.csv", header = TRUE, row.names = 1) #weights40 <- read.csv("connectivity 40km_June2019/settlement_radius_40km/edges_weights40.csv", header = TRUE, row.names = 1) weights.m = as.matrix(weights$average) weight <- as.vector(weights.m) #normalised weights? #x <- weight #normalized = (x-min(x))/(max(x)-min(x)) #normalized #unweighted AEM vectors #cockles.aem.uw <- aem(binary.mat = nodes.m, rm.link0 = FALSE) #cockles.aem.uw # weighted AEM vectors construction cockles.aem <- aem(binary.mat = nodes.m, weight = weight, rm.link0 = FALSE) cockles.aem AEM <- as.data.frame(cockles.aem$vectors) colnames(AEM) <- c("AEM1", "AEM2", "AEM3", "AEM4") rownames(AEM) <- c("DEE", "RW", "DY", "DA", "BU", "BB", "FS" ) AEM ############################################## ########## TEMPERATURE DATA ################## ############################################## # surface (sst) and bottom (sbt) temperatures env <- read.csv("connectivity_newsimulations_Feb2019/environment.csv", header = TRUE, row.names = 1) ENV <- env[, -c(1,2)] ######################################## ########## GEOGRAPHICAL DATA ########### dbMEMs ######################################## ### Import geo object (coordinates) geo <- env[,1:2] geo$CAR <- geoXY(geo$lat, geo$lon) # Transform geographic coordinates into cartesian coordinates euclidian_distances <- dist(geo$CAR, method="euclidean") # Calculate euclidian distances based on geographic coordinates pcnm_values <- pcnm(euclidian_distances, thresh= NULL) # Perform PCNM function on euclidian distances geo <- cbind(geo, pcnm_values$vectors) ### Keep only cartesian coordinates and mem vectors GEO <- geo[,-c(1:2)] rownames(GEO) <- c("DEE", "RW", "DY", "DA", "BU", "BB", "FS" ) GEO ####################################### ########## GENOMIC DATA ########### ####################################### datamaf <- read.table('MAF_per_population.txt', header=TRUE, sep = '\t') datamaf ncol(datamaf) maf.out <- datamaf[, c(96, 97, 248, 291,351,548,635,797,811,927,953,1027,1301,1622)] maf.neu <- datamaf[, -c(96, 97, 248, 291,351,548,635,797,811,927,953,1027,1301,1622)] ###### do a dRDA using Fst ######################### #allpfst <- as.dist(as.matrix(read.csv('../../RELAXED COND ANALYSIS 2017/cockles_1725/all_pfst.csv', header = TRUE, row.names = 1)), upper=TRUE) #allpfst #neupfst <- as.dist(as.matrix(read.csv('../../RELAXED COND ANALYSIS 2017/cockles_1725/neu_pfst.csv', header = TRUE, row.names = 1)), upper=TRUE) #neupfst #testing gen distances with AEM - Do PCA first #neufst.pca <- prcomp(neupfst) #summary(neufst.pca) #neu.rda=neufst.pca$x[, 1:4] #outpfst <- as.dist(as.matrix(read.csv('../../RELAXED COND ANALYSIS 2017/cockles_1725/out_pfst.csv', header = TRUE, row.names = 1)), upper=TRUE) #outpfst ####################################################### ########################### DECOSTAND GENO <- decostand(datamaf, "hel") NEU <- decostand(maf.neu, "hel") OUT <- decostand(maf.out, "hel") #GENOFST <- decostand(allpfst, "hel") #NEUFST <- decostand(neupfst, "hel") #OUTFST <- decostand(outpfst, "hel") #rownmaes have to match #match rownames with fst matrix c("BB", "DA", "DEE", "DY", "FS", "RW" , "BU") ENV <- ENV[match(rownames(NEU), rownames(ENV)),] GEO <- GEO[match(rownames(NEU), rownames(GEO)), ] AEM <- AEM[match(rownames(NEU), rownames(AEM)), ] #now check identical(rownames(OUT), rownames(ENV)) identical(rownames(NEU), rownames(AEM)) identical(rownames(OUT), rownames(GEO)) #identical(rownames(GENOFST), rownames(ENV)) ######################################################### #################################### ### Do a RDA for ENV #################################### env.neu.rda <- rda(NEU ~ ., data=ENV, scale = T) set.seed(10); mod.neu <- rda(NEU ~ 1, data = ENV) #set up empty model env.neu.step <- ordistep(mod.neu, scope = formula(env.neu.rda), direction = 'both', permutations = how(nperm=10000)) env.neu.step$anova #July meanSBTG P=0.002 (global) RsquareAdj(env.neu.step) #adjr2=0.1 screeplot(env.neu.rda) ################################################################ env.out.rda <- rda(OUT ~ ., data=ENV, scale = T) set.seed(634); mod.out <- rda(OUT ~ 1, data = ENV) #set up empty model env.out.step <- ordistep(mod.out, scope = formula(env.out.rda), direction = 'both', permutations = how(nperm=10000)) env.out.step$anova #Apr minSST P=0.025; Sept minSST P=0.006 RsquareAdj(env.out.step) #adjr2=0.71 screeplot(env.out.rda) ################################################################ ################################################################ ################################################################ #################################### ### Do a RDA for AEMs #################################### aem.neu.rda <- rda(NEU ~ ., data=AEM, scale = T) set.seed(10); mod.neu.2 <- rda(NEU ~ 1, data = AEM, scale = T) aem.neu.step <- ordistep(mod.neu.2, scope = formula(aem.neu.rda), direction = 'both', permutations = how(nperm=10000)) aem.neu.step$anova # AEM2 P=0.047 RsquareAdj(aem.neu.step) #adjr2=0.04 screeplot(aem.neu.rda) #plotting biplot RDA AEM for neutral markers ggord(aem.neu.step, arrows = TRUE, geom = "text", legend = "none") plot(aem.neu.step, scaling= 3, display='sites') ######################################################################## aem.out.rda <- rda(OUT ~ ., data=AEM, scale = T) set.seed(30); mod.out.2 <- rda(OUT ~ 1, data = AEM, scale = T) aem.out.step <- ordistep(mod.out.2, scope = formula(aem.out.rda), direction = 'both', permutations = how(nperm=10000)) aem.out.step$anova #NS RsquareAdj(aem.out.step) screeplot(aem.out.rda) ############################################################################## ############################################################################## ### Do a RDA for GEO geo.neu.rda <- rda(NEU ~ ., data = GEO) set.seed(10); mod3.neu <- rda(NEU ~ 1, data = GEO, scale =T) geo.neu.step <- ordistep(mod3.neu, scope = formula(geo.neu.rda), direction = 'both', permutations = how(nperm = 10000)) geo.neu.step$anova #PCNM1 P=0.0006 RsquareAdj(geo.neu.step) #adjr2=0.07 screeplot(geo.rda) ############################################################################## geo.out.rda <- rda(OUT ~ ., data = GEO) set.seed(10); mod3.out <- rda(OUT ~ 1, data = GEO, scale =T) geo.out.step <- ordistep(mod3.out, scope = formula(geo.out.rda), direction = 'both', permutations = how(nperm = 10000)) geo.out.step$anova #CAR RsquareAdj(geo.out.step) screeplot(geo.out.rda) ############################################################################## ############################################################################## # final RDAs data <- as.data.frame(cbind(GEO$PCNM1, ENV$July_meanSBT, ENV$April_minSST, ENV$Sept_minSST, AEM$AEM2)) colnames(data) <- c("PCNM1", "JUL_meanSBTG", "Apr_minSST", "Sept_minSST", "AEM2") dataneu <- data[, c(1,2,5)] dataout <- data[, 3:4] #global RDA NEU global.rda.neu <- rda(NEU ~ AEM2+JUL_meanSBTG+PCNM1, data=dataneu) anova(global.rda.neu, step=10000) RsquareAdj(global.rda.neu) #r2a=0.08 p=0.1 ggord(global.rda.neu) + theme(aspect.ratio=1) #square graph summary(global.rda.neu) #global RDA OUT global.rda.out <- rda(OUT ~., data=dataout) anova(global.rda.out,step=10000) #P=0.003 RsquareAdj(global.rda.out) #adjr2=0.71 ggord(global.rda.out) + theme(aspect.ratio=1) summary(global.rda.out) ############################################################################## #NEU PARTIAL RDA, condition ENV partial.neu.rda <- rda(NEU ~ ENV$July_meanSBT + Condition(GEO$PCNM1+ AEM$AEM2 ), scale = T) anova(partial.neu.rda, permutations = 10000) RsquareAdj(partial.neu.rda) screeplot(partial.neu.rda) ######## PARTIAL RDA ######## #NEU of ENV controlling for spatial partialspa.neu.rda <- rda(NEU ~ data$JUL_meanSBT + data$PCNM1, scale = T) anova(partialspa.neu.rda, permutations = 1000) RsquareAdj(partialspa.neu.rda) #NEU of Spatial controlling for ENV partialenv.neu.rda <- rda(NEU ~ data$PCNM1 + data$JUL_meanSBT, scale = T) anova(partialenv.neu.rda, permutations = 1000) RsquareAdj(partialenv.neu.rda) ############################################################################## # plot leg <- c('Bannow Bay', 'Gann', 'Dee', 'Dyfi', 'Flaxfort Strand', 'Red Wharf', 'Burry Inlet') col <- c('forestgreen', 'red1', 'steelblue1', 'orange', 'yellowgreen', 'steelblue4', 'darkred') pops <- c('BAN', 'GAN', 'DEE', 'DYF', 'FLX', 'RWB', 'BUR') library(tidyverse) #### NEU ##################################### load.rda <- scores(global.rda.neu, choices=c(1:2), display="species") sum.neu <- summary(global.rda.neu) neu.sco <- as.data.frame(scores(global.rda.neu)$sites) neu.sco neu.sco <- cbind(neu.sco, pops) neu.sco neudata <- data.frame(sum.neu$biplot[,1:2]) p1 <- ggplot() + geom_hline(yintercept = 0, linetype = 'dashed') + geom_vline(xintercept = 0, linetype = 'dashed') + theme_bw() + xlab('RDA1 (47.63%)') + ylab('RDA2 (28.68%)') + geom_point(data = neu.sco, aes(x=RDA1, y=RDA2, fill = as.factor(pops)), colour = 'black', bg = col, shape = 21, size = 4) + coord_fixed() neuarrows <- attributes(scores(global.rda.neu))$const neuplot <- plot(global.rda.neu) narr <- attributes(neuplot$biplot)$arrow.mul #arr <- ordiArrowMul(scores(final.neu.rda, choices = 1:2, display = 'bp')) #same as above neu.final.plot <- p1 + geom_segment(data=neudata, aes(x=0, xend=RDA1*narr, y=0, yend=RDA2*narr), color="darkgrey", arrow=arrow(length=unit(0.25, "cm"))) + geom_text(aes(x=-0.62, y=-0.10, label = 'July mean SBTG'), colour = 'darkgrey') + geom_text(aes(x=0.7, y=-0.1, label = 'PCNM1'), colour = 'darkgrey') + geom_text(aes(x=-0.62, y=-0.42, label = 'AEM2'), colour = 'darkgrey') neu.final.plot <- neu.final.plot + theme(aspect.ratio = 1) ############################################## #### OUT ##################################### sum.out <- summary(global.rda.out) out.sco <- as.data.frame(scores(global.rda.out)$sites) out.sco out.sco <- cbind(out.sco, pops) out.sco outdata <- data.frame(sum.out$biplot[,1:2]) p3 <- ggplot() + geom_hline(yintercept = 0, linetype = 'dashed') + geom_vline(xintercept = 0, linetype = 'dashed') + theme_bw() + xlab('RDA1 (59.74%)') + ylab('RDA2 (40.26%)') + geom_point(data = out.sco, aes(x=RDA1, y=RDA2, fill = as.factor(pops)), colour = 'black', bg = col, shape = 21, size = 4) + coord_fixed() df3 <- data.frame(sum.out$biplot[,1:2]) #arrows <- attributes(scores(final.neu.rda))$const outplot <- plot(global.rda.out) arr <- attributes(outplot$biplot)$arrow.mul #arr <- ordiArrowMul(scores(final.neu.rda, choices = 1:2, display = 'bp')) #same as above out.final.plot <- p3 + geom_segment(data=df3, aes(x=0, xend=RDA1*arr, y=0, yend=RDA2*arr), color="darkgrey", arrow=arrow(length=unit(0.25, "cm"))) + geom_text(aes(x=-1.5, y=-1.1, label = 'April Min SST'), colour = 'darkgrey') + geom_text(aes(x=0.8, y=-1.88, label = 'September Min SST'), colour = 'darkgrey') + theme(aspect.ratio=1) out.final.plot <- out.final.plot + theme(aspect.ratio=1) ggord(global.rda.out) ######################################### library(gridExtra) grid.arrange(neu.final.plot, out.final.plot , ncol=2)