#load packages library(vegan) library(adegenet) library(adespatial) library(spdep) setwd("G:/My Drive/BobcatData/") Bobcats<-read.csv("Bobcat_Genotypes_All_Edited.csv", header = TRUE) Bobcats <- Bobcats[1:102,1:19] row.names(Bobcats) <- Bobcats$..CatID Bobcats$..CatID <- NULL ####################sPCA, conducted in adegenet###################### #adegenet takes genepop files, it MUST be space delimited with a .gen extension# Bobcat<-read.genepop("Bobcat102.gen", ncode=3L, quiet=FALSE) Bobcat #Read XY coordinates, they must all be unique, use the jitter function# XY<-read.table(file="XY_102test.txt", header=T) # Matrix needs diff name than coordinates# xy.jitter<-matrix(ncol=2, nrow=102) #The XY[,1] is adding noise to the first column (i.e., X coordinates) and XY[,2] is the second column/Y coordinates# xy.jitter[,1]<-jitter(XY[,1], a=0.01) xy.jitter[,2]<-jitter(XY[,2], a=0.01) #sPCA function, nfposi is the number of positive axes to retain. Take 2 positive and one negative axis# mySpca<-spca(Bobcat, xy=xy.jitter, nfposi=2, nfnega=1) #Big figure# plot(mySpca) #Plot of each axis vs. spatial autocorrelation (see guide for this interpretation)# screeplot(mySpca) summary(mySpca) View(mySpca$li[,1]) pc1.mctest <- moran.mc(mySpca$li[,1], 999) plot(pc1.mctest) #Colorplot to map principal components onto geographic space (not used in manuscript)# colorplot(mySpca,cex=3,main="colorplot of mySpca, first global score") #Scores for resultant sPCA, these are the response variables in the RDA. Can copy and paste them from R into a text document, or write a table in R# mySpca$li View(mySpca$li) View(mySpca$tab) summary.spca ###Which markers contributed to each Axis#### myLoadings <- mySpca$c1[,1]^2 boxplot(myLoadings~mySpca$tab.fac, las=3, ylab="Contribution", xlab="Marker", main="Contributions by markers \nto the first global score", col="grey") ####################### DAPC analysis, conduted in adegenet################################ #find.clusters first transforms the data using PCA, asking the user to specify the number of retained PCs grp <- find.clusters(Bobcat, max.n.clust=30) names(grp) grp grp$size table(pop(Bobcat), grp$grp) #cross-validation to determine number of principal components to use #k-means: use as many PCs as need, DAPC: want to minimize PCs, overfitting bad set.seed(999) Bobcatx <- xvalDapc(tab(Bobcat, NA.method = "mean"), grp$grp) Bobcatx[-1] #dapc, no prior grouping dapc1 <- dapc(Bobcat, grp$grp) dapc1 scatter(dapc1) dapc1$var.contr myCol <- c("blue","yellow","red","orange","darkblue","purple") scatter(dapc1,1,1, col=myCol, bg="white", scree.da=FALSE, legend=TRUE, solid=.8) contrib <- loadingplot(dapc1$var.contr, axis=2,thres=.07, lab.jitter=1) #Interpreting group memberships class(dapc1$posterior) dim(dapc1$posterior) round(head(dapc1$posterior),3) round(dapc1$posterior) summary(dapc1) assignplot(dapc1, subset= 1:50) compoplot(dapc1, subset=1:50, posi="bottomright", txt.leg=paste("Cluster", 1:3), lab="", n.col=1, xlab="individuals", col=funky(3)) #admixed individuals temp <- which(apply(dapc1$posterior,1, function(e) all(e<0.9))) temp compoplot(dapc1, subset=temp, posi="bottomright", txt.leg=paste("Cluster", 1:3), n.col=2,col=funky(3)) #####################Redundancy Analysis, conducted in vegan############################### #Bobcat Input Files: Dispersal and Home Range# HR<-read.table(file="HR_EnvData_102cats.txt", header=T) Disp<-read.table(file="Disp_EnvData_102cats.txt", header=T) #HR RDA analysis for HR. Each column is a variable. In this case, the order is: 2: X coordinate, 3: Y coordinate, 4: LULC_URban, 5: LULC_Agriculture, 6: LULC_HerbaceousRangeland, 7: Mixed Rangeland, 8: Stream Density, 9: sPCA Axis 1 Scores, 10: sPCA Axis 2 Scores# try<-rda(HR[,9]+HR[,10]~HR[,4]+HR[,5]+HR[,6]+HR[,7]+HR[,8]+Condition(HR[,2]+HR[,3])) #Setting up model selection for stepwise (ordistep). I use a null model mod0 and the full model mod1 to give the function the search scope. It then adds each variable combination, checks the AIC and rejects if it is not better. It then gives the best model plus a long output# mod0<-rda(HR[,9]+HR[,10] ~ 1) mod1<-try sel<-ordistep(mod0, scope = formula(mod1), perm.max = 1000) plot(HR[,6],HR[,9])##sPCA Axis1 plot(HR[,6],HR[,10]) ##sPCA Axis 2 plot(HR[,5],HR[,9]) plot(HR[,5],HR[,10]) plot(mod1) #Analysis variance of final selection model# RDAc<-anova.cca(sel, by="terms") RDAc #Rsquared Adjusted# RsquareAdj(sel) #Dispersal analysis for HR. Same as before, but order is adjusted due to the addition of Mesquite Juniper Bush. 2: X coordinate, 3: Y coordinate, 4: LULC_Urban, 5: LULC_Agriculture, 6: LULC_HerbaceousRange, 7: Mesquite_Juniper, 8: LULC_MixedRangeland, 9: Stream Density, 10: sPCA Axis 1 Scores, 11: sPCA Axis 2 Scores# try<-rda(Disp[,10]+Disp[,11]~Disp[,4]+Disp[,5]+Disp[,6]+Disp[,7]+Disp[,8]+Disp[,9]+Condition(Disp[,2]+Disp[,3])) mod0<-rda(Disp[,10]+Disp[,11] ~ 1) mod1<-try sel<-ordistep(mod0, scope = formula(mod1), perm.max = 1000) plot(Disp[,8],Disp[,10]) plot(Disp[,8],Disp[,11]) plot(Disp[,5], Disp[,10]) plot(Disp[,5],Disp[,11]) RDAc<-anova.cca(sel, by="terms") RDAc RsquareAdj(sel)