library("adegenet") library("vegan") library("gdistance") library("rgdal") library("raster") library("PopGenReport") library("ade4") library("hierfstat") #load coord latlong<-read.csv ("loc_pop.csv",header = T) geoloc<-as.matrix(latlong[,2:3]) geoloc2<-as.matrix(read.table("Coord_corrected_for_IBR-road.txt",header=T)) ####Extract climate variables, with only bio1 as example bio1<-extract(raster("./Climate/wc2.0_bio_30s_01.tif"), geoloc)# Annual Mean Temperature ####Combine all extracted variables envdata<-cbind(bio1,bio2,bio3,bio4,bio5,bio6,bio7,bio8,bio9,bio10,bio11,bio12,bio13,bio14, bion,bio16,bio17,bio18,bio19) ####load resistance surface layers #IBR-climate cost <-raster("./resistance surface/IBR-climate.tif") trcost<- geoCorrection(transition(cost,mean,8,symm=F),type="c") #IBR-landcover cost2 <-raster("./resistance surface/IBR-landcover.tif") trcost2<- geoCorrection(transition(cost2,mean,8,symm=F),type="c") #IBR-road cost3 <-raster("./resistance surface/IBR-road.tif") trcost3<- geoCorrection(transition(cost3,mean,8,symm=F),type="c") ####Distance calculation library("gdistance") #cost distance costdis <-costDistance(trcost,geoloc) costdis2<-costDistance(trcost2,geoloc) costdis3<-costDistance(trcost3,geoloc2) #Environmental distance pcaenv<-dudi.pca(envdata, center=TRUE, scannf = F,nf=2,scale=T) envdis<- dist(data.frame(pcaenv$li)[1:2],method = "euclidean") #Geographic distance geodis<- as.dist(pointDistance(SpatialPoints(geoloc), lonlat=T,allpairs = T)) #fst/1-fst genepop<-read.genepop("gennotime.gen",ncode = 3L) fst<-pairwise.fst(genepop, res.type = "dist") gendis<-fst/(1-fst) #nei's D neid<-dist.genpop(genind2genpop(genepop),method=1) #########Test,only fst/(1-fst) as example######## ##MMRR library("PopGenReport") IBDMMRR<-lgrMMRR(gen.mat=as.matrix(gendis),cost.mats=list(as.matrix(geodis)),nperm=9999) IBEMMRR<-lgrMMRR(gen.mat=as.matrix(gendis),cost.mats=list(as.matrix(envdis)),nperm=9999) IBR1MMRR<-lgrMMRR(gen.mat=as.matrix(gendis),cost.mats=list(as.matrix(costdis)),nperm=9999) IBR2MMRR<-lgrMMRR(gen.mat=as.matrix(gendis),cost.mats=list(as.matrix(costdis2)),nperm=9999) IBR3MMRR<-lgrMMRR(gen.mat=as.matrix(gendis),cost.mats=list(as.matrix(costdis3)),nperm=9999) IBDMMRR IBEMMRR IBR1MMRR IBR2MMRR IBR3MMRR ##partial Mantel TEST library('vegan') PIBDtest<- mantel.partial(gendis,geodis,envdis,method="pearson", permutations=9999) PIBEtest<- mantel.partial(gendis,envdis,geodis,method="pearson", permutations=9999) PIBRtest1<- mantel.partial(gendis,costdis,geodis,method="pearson", permutations=9999) PIBRtest2<- mantel.partial(gendis,costdis2,geodis,method="pearson", permutations=9999) PIBRtest3<- mantel.partial(gendis,costdis3,geodis,method="pearson", permutations=9999) c(PIBDtest$statistic,PIBDtest$signif) c(PIBEtest$statistic,PIBEtest$signif) c(PIBRtest1$statistic,PIBRtest1$signif) c(PIBRtest2$statistic,PIBRtest2$signif) c(PIBRtest3$statistic,PIBRtest3$signif) ###Map least-cost path n=nrow(geoloc) pdf("IBR-climate.pdf") plot(cost,xlim=c(74,90.9),ylim=c(41,49),xaxs="i",yaxs="i", maxpixels=1e7,horizontal=T) title(main="(a) IBR-climate") plot(SpatialPoints(geoloc),add=T,pch=20,col="red") i <- 1 while (i <= n){ j <- i+1 while (j <= n) { plot(shortestPath(trcost, geoloc[i,], geoloc[j,], output="SpatialLines"),col="#666666", add=TRUE) j <- j+1} i <- i+1} dev.off() pdf("IBR-landcover.pdf") plot(cost2,xlim=c(74,90.9),ylim=c(41,49),xaxs="i",yaxs="i", maxpixels=1e7,horizontal=T) title(main="(b) IBR-landcover") plot(SpatialPoints(geoloc),add=T,pch=20,col="red") i <- 1 while (i <= n){ j <- i+1 while (j <= n) { plot(shortestPath(trcost2, geoloc[i,], geoloc[j,], output="SpatialLines"),col="#666666", add=TRUE) j <- j+1} i <- i+1} dev.off() pdf("IBR-road.pdf") plot(cost3,xlim=c(74,90.9),ylim=c(41,49),xaxs="i",yaxs="i", legend=T, maxpixels=1e7,horizontal=T) title(main="(c) IBR-road") plot(SpatialPoints(geoloc2),add=T,pch=20,col="red") i <- 1 while (i <= n){ j <- i+1 while (j <= n) { plot(shortestPath(trcost3, geoloc2[i,], geoloc2[j,], output="SpatialLines"),col="#666666", add=TRUE) j <- j+1} i <- i+1} dev.off() #######GDM######### library("gdm") #Readers should carefully read the GDM input format requirements #Geographic distance pcostdis<-as.matrix(read.csv("./gdm/geodis.csv",header = T)) #Climate data envdata2<-read.csv("envdata.csv",header = T) envdata3<-dudi.pca(envdata2[,2:20],scannf = F, nf = 2)$li preenv<-cbind(envdata2$site,envdata3,envdata2$longitude,envdata2$latitude) names(preenv)<-c("site","Axis1","Axis2","envdata2.longitude","envdata2.latitude") #three kinds of cost distance prere1<-as.matrix(read.csv("./gdm/resis1.csv",header = T)) prere2<-as.matrix(read.csv("./gdm/costdis2.csv",header = T)) prere3<-as.matrix(read.csv("./gdm/costdis3.csv",header = T)) #genetic data gdmgendis<-as.matrix(read.csv("./gdm/gendis.csv")) gdmneid<-as.matrix(read.csv("./gdm/neid.csv")) #fake data prefake<-as.matrix(read.csv("fake.csv",header = T)) #GDM geodis gdm.geotab<-subset(formatsitepair(bioData=gdmgendis,bioFormat=3,siteColumn = "site",XColumn ="envdata2.longitude",YColumn="envdata2.latitude", predData = prefake, distPreds = list(pcostdis)),select = -c(s1.fake,s2.fake)) #GDM climate gdm.envtab<- formatsitepair(bioData=gdmgendis,bioFormat=3,siteColumn = "site",XColumn ="envdata2.longitude",YColumn="envdata2.latitude", predData = preenv) #gdm costdis1 gdm.re1tab<-subset(formatsitepair(bioData=gdmgendis,bioFormat=3,siteColumn = "site",XColumn ="envdata2.longitude",YColumn="envdata2.latitude", predData = prefake, distPreds = list(prere1)),select = -c(s1.fake,s2.fake)) #gdm costdis2 gdm.re2tab<-subset(formatsitepair(bioData=gdmgendis,bioFormat=3,siteColumn = "site",XColumn ="envdata2.longitude",YColumn="envdata2.latitude", predData = prefake, distPreds = list(prere2)),select = -c(s1.fake,s2.fake)) #gdm costdis3 gdm.re3tab<-subset(formatsitepair(bioData=gdmgendis,bioFormat=3,siteColumn = "site",XColumn ="envdata2.longitude",YColumn="envdata2.latitude", predData = prefake, distPreds = list(prere3)),select = -c(s1.fake,s2.fake)) #gdm all gdm.alltab<- formatsitepair(bioData=gdmgendis,bioFormat=3,siteColumn = "site",XColumn ="envdata2.longitude",YColumn="envdata2.latitude", predData = preenv, distPreds = list(pcostdis,prere1,prere2,prere3)) #gdm geo+env gdm.geoenvtab<- formatsitepair(bioData=gdmgendis,bioFormat=3,siteColumn = "site",XColumn ="envdata2.longitude",YColumn="envdata2.latitude", predData = preenv,distPreds = list(pcostdis)) #gdm env+re1 gdm.envre1tab<- formatsitepair(bioData=gdmgendis,bioFormat=3,siteColumn = "site",XColumn ="envdata2.longitude",YColumn="envdata2.latitude", predData = preenv, distPreds = list(prere1)) #gdm env+re2 gdm.envre2tab<- formatsitepair(bioData=gdmgendis,bioFormat=3,siteColumn = "site",XColumn ="envdata2.longitude",YColumn="envdata2.latitude", predData = preenv, distPreds = list(prere2)) #gdm env+re3 gdm.envre3tab<- formatsitepair(bioData=gdmgendis,bioFormat=3,siteColumn = "site",XColumn ="envdata2.longitude",YColumn="envdata2.latitude", predData = preenv, distPreds = list(prere3)) ###GDM FIT gdm.geo<-gdm(gdm.geotab) gdm.env<-gdm(gdm.envtab) gdm.re1<-gdm(gdm.re1tab) gdm.re2<-gdm(gdm.re2tab) gdm.re3<-gdm(gdm.re3tab) gdm.all<-gdm(gdm.alltab) gdm.geo_env<-gdm(gdm.geoenvtab) gdm.envre1<-gdm(gdm.envre1tab) gdm.envre2<-gdm(gdm.envre2tab) gdm.envre3<-gdm(gdm.envre3tab) ##GDM SUMMARY summary(gdm.geo) summary(gdm.env) summary(gdm.re1) summary(gdm.re2) summary(gdm.re3) summary(gdm.all) summary(gdm.geo_env) summary(gdm.envre1) summary(gdm.envre2) summary(gdm.envre3) #######RDA###### #Individual coord. indloc<-as.matrix(read.csv("loc_ind.csv",header = T)[,2:3]) #Extract individual-level climatic variables, with only bio1 as example ibio1<-extract(raster("./Climate/wc2.0_bio_30s_01.tif"), indloc) #Combine climatic variable indenv<-cbind(ibio1,ibio2,ibio3,ibio4,ibio5,ibio6,ibio7,ibio8, ibio9,ibio10,ibio11,ibio12,ibio13,ibio14,ibion,ibio16,ibio17, ibio18,ibio19) #Individual least-cost distance indres <- scores(pcnm(as.matrix(read.csv("ind-LCD.csv",header=T,row.names = 1)))) keep<-round(length(which(indres$values>0))/2) indres<-scores(indres)[,1:keep] #Convert individual geographic distance to PCNM indgeo <- pcnm(as.matrix(as.dist(pointDistance(SpatialPoints(indloc),lonlat=T,allpairs = T)))) keep<-round(length(which(indgeo$values>0))/2) indgeo<-scores(indgeo)[,1:keep] #vif-test ie <- subset(indenv,select = -c(ibio3,ibio4,ibio18,ibio10,ibio11,ibio17,ibio16, ibio14,ibio13,ibio5,ibio7,ibio19,ibio9)) rda_test<-dbrda(idis ~ie, scale=T) vif.cca(rda_test) #rda fitting rda_env <- dbrda(idis ~ie,scale=T) rda_res2 <- dbrda(idis ~indres2,scale=T) rda_env_res2 <- dbrda(idis ~ie + Condition(indres2),scale=T) rda_res2_env <- dbrda(idis ~indres2 + Condition(ie),scale=T) rda_envres2 <- dbrda(idis ~ie + indres2 ,scale=T) summary(rda_env)$constr.chi/summary(rda_env)$tot.chi summary(rda_res2)$constr.chi/summary(rda_res2)$tot.chi summary(rda_env_res2)$constr.chi/summary(rda_env_res2)$tot.chi summary(rda_res2_env)$constr.chi/summary(rda_res2_env)$tot.chi RsquareAdj(rda_env) RsquareAdj(rda_res2) RsquareAdj(rda_env_res2) RsquareAdj(rda_res2_env) anova(rda_env) anova(rda_res2) anova(rda_env_res2) anova(rda_res2_env) ####Draw RDA result, with only envionment fitting as example #Read the prepered Q-values matrix from STRUCTURE result to color the RDA points. indcolor<-as.matrix(read.csv("loc_ind.csv",header = T)[,4]) #Draw figure pdf("rda_env.pdf") plot(rda_env,choices = c(1, 2), type="n",xlab="RDA1",ylab="RDA2") points(rda_env, display= "sites", pch=20, col=rgb(color(indcolor),alpha=190,max=255)) text(rda_env, scaling=2, display="bp", col="black", cex=1.2,labels=c("Bio1","Bio2","Bio6","Bio8","Bio12","Bion")) dev.off()