#R script by Sam Banks to calculate landscape-wide summary stats of geentic diversity and associations with landscape fire variables across all simulation outputs #Can also use selected bits of the Gen_analyses_2 function (not as a defined function) to get per-cell genetic and fire data for individual scenarios for plotting #setwd("foldername") #set foldername to folder containing all simulation output files containing individual genotypes filelist<-list.files("foldername", pattern=".csv") #again, fix up folder name to allow reading of all .csv file names in relevant folder #Try one file as example #x<-"2000_TP01_WR01_MD01_SP01_SV02_RP04.csv" #just used a single replicate to test function #remove(x) #Define function to get basic genetic diversity stats Gen_analyses_2 <- function(x) { scenario<-x data <- read.csv(x, sep="\t",header=TRUE) #Name pops data$Pop<-paste("Cell",data$LocID,sep="_") #Create aggregated 3x3 cell populations within observation window data$locx3<-floor((data$locx-14)/3) data$locy3<-floor((data$locy-161)/3) data$Pop3<-paste("Cell",data$locx3,data$locy3,sep="_") #data frame of Pop3 info Agg3Table<-aggregate(data[, 9:10], list(Pop3=data$Pop3), mean) Agg3popsize<-aggregate(data[, 1], list(Pop3=data$Pop3), length) Agg3tb<-data.frame(Agg3Table) Agg3ps<-data.frame(Agg3popsize) Agg3df<-cbind(Agg3tb,Agg3ps$x) #####Allele Freqs for selected locus Agg3A31Acount<-aggregate(data[,41], list(Pop3=data$Pop3), sum) Agg3df<-merge(Agg3df,Agg3A31Acount,by="Pop3") names(Agg3df)[names(Agg3df)=="x"] <- "Sel_allele_count" names(Agg3df)[names(Agg3df)=="Agg3ps$x"] <- "Popsize" Agg3df$sel_alf<-(Agg3df$Sel_allele_count/(Agg3df$Popsize*2)) Agg3df$sel_alf2<-1-((Agg3df$Sel_allele_count/(Agg3df$Popsize*2))) #Can plot this to show maps of selected allele freqs for each scenario Agg3df$He_sel<-(1-((Agg3df$sel_alf)^2)-(((1-Agg3df$sel_alf))^2)) ######### #Hs_neutral Agg3neutralcounts<-aggregate(data[,11:40], list(Pop3=data$Pop3), sum) Agg3neutralcounts<-Agg3neutralcounts[,2:31] neutralf1<-Agg3neutralcounts/(Agg3df$Popsize*2) neutralf2<-1-neutralf1 Agg3neutral_He<-1-(neutralf1^2)-(neutralf2^2) #Can plot this to show maps of neutral Hs for each scenario meanHe<-rowMeans(Agg3neutral_He) Agg3df<-cbind(Agg3df,meanHe) ########## #Neutral FST and HT across entire landscape #HT first Alf_loc<-colSums(data[,11:40])/(2*(length(data[,11]))) HT_loc<-1-(Alf_loc^2)-((1-Alf_loc)^2) meanHT<-mean(HT_loc) #Mean HS meanHS<-colMeans(Agg3neutral_He) FST<-(meanHT-meanHS)/meanHT meanFST<-mean(FST) #Linear models of per-cell genetic diversity against mean IFI and habitat selection variable (scaled time since fire) #1 linear model of "early successional" selected allele freq against habitat selection variable lm_Ph_Hab<-lm(Agg3df$sel_alf~Agg3df$Hab) R2_Ph_Hab<-summary(lm_Ph_Hab)$r.squared SLP_Ph_Hab<-summary(lm_Ph_Hab)$coefficients[2,1] cor_Ph_Hab<-cor(Agg3df$sel_alf,Agg3df$Hab) #2 linear model of "late successional" selected allele freq against habitat selection variable lm_Ph2_Hab<-lm(Agg3df$sel_alf2~Agg3df$IDI) R2_Ph2_Hab<-summary(lm_Ph2_Hab)$r.squared SLP_Ph2_Hab<-summary(lm_Ph2_Hab)$coefficients[2,1] cor_Ph2_Hab<-cor(Agg3df$sel_alf2,Agg3df$Hab) #3 linear model of "early successional" selected allele freq against mean IFI lm_Ph_IDI<-lm(Agg3df$sel_alf~Agg3df$IDI) R2_Ph_IDI<-summary(lm_Ph_IDI)$r.squared SLP_Ph_IDI<-summary(lm_Ph_IDI)$coefficients[2,1] cor_Ph_IDI<-cor(Agg3df$sel_alf,Agg3df$IDI) #4 linear model of "late successional" selected allele freq against mean IFI lm_Ph2_IDI<-lm(Agg3df$sel_alf2~Agg3df$IDI) R2_Ph2_IDI<-summary(lm_Ph2_IDI)$r.squared SLP_Ph2_IDI<-summary(lm_Ph2_IDI)$coefficients[2,1] cor_Ph2_IDI<-cor(Agg3df$sel_alf2,Agg3df$IDI) #5 linear model of neutral locus Hs against habitat selection variable lm_He_Hab<-lm(Agg3df$meanHe~Agg3df$Hab) R2_He_Hab<-summary(lm_He_Hab)$r.squared SLP_He_Hab<-summary(lm_He_Hab)$coefficients[2,1] cor_He_Hab<-cor(Agg3df$meanHe,Agg3df$Hab) #6 linear model of neutral locus Hs against mean IFI lm_He_IDI<-lm(Agg3df$meanHe~Agg3df$IDI) R2_He_IDI<-summary(lm_Ph_IDI)$r.squared SLP_He_IDI<-summary(lm_Ph_IDI)$coefficients[2,1] cor_He_IDI<-cor(Agg3df$meanHe,Agg3df$IDI) out2<-as.data.frame(cbind(scenario,R2_Ph_Hab, SLP_Ph_Hab, cor_Ph_Hab, R2_Ph2_Hab, SLP_Ph2_Hab, cor_Ph2_Hab, R2_Ph_IDI, SLP_Ph_IDI, cor_Ph_IDI, R2_Ph2_IDI, SLP_Ph2_IDI, cor_Ph2_IDI,R2_He_Hab,SLP_He_Hab,cor_He_Hab, R2_He_IDI,SLP_He_IDI,cor_He_IDI, meanHT, meanFST)) return(out2) } ######################################################################## #Apply it to folder of all simulation output genotypes out2 <- lapply(filelist, Gen_analyses_2) outtable2<-as.data.frame(do.call(rbind, out2)) out1 write.table(outtable2,"Gen_summaries.csv",sep=",") ######################################################################## #Worked off the summary output file for subsequent analyses to avoid repeating the above calcs #Added Fisher Z-transformed correlation coefficients to the table for further analysis require(stringr) regtable<-read.csv("Gen_summaries.csv", heade=TRUE) regtable$scenario<-str_replace_all(regtable$scenario, "_", " ") scens<-str_split_fixed(regtable$scenario, " ", 7) regtable<-cbind(scens,regtable) nms<- c("Yr","Topo", "Weather", "Dispersal", "Selection", "Severity","Replicate", "scenario","R2_Ph_Hab","SLP_Ph_Hab","cor_Ph_Hab", "R2_Ph2_Hab", "SLP_Ph2_Hab", "cor_Ph2_Hab", "R2_Ph_IDI", "SLP_Ph_IDI","cor_Ph_IDI","R2_Ph2_IDI", "SLP_Ph2_IDI", "cor_Ph2_IDI","R2_He_Hab","SLP_He_Hab", "cor_He_Hab","R2_He_IDI", "SLP_He_IDI", "cor_He_IDI", "meanHT", "meanFST") require(data.table) setnames(regtable,nms) # ARCTANH transformed correlation coeffs for analysis regtable$Z_cor_Ph_Hab<-atanh(regtable$cor_Ph_Hab) regtable$Z_cor_Ph_IDI<-atanh(regtable$cor_Ph_IDI) regtable$Z_cor_Ph2_Hab<-atanh(regtable$cor_Ph2_Hab) regtable$Z_cor_Ph2_IDI<-atanh(regtable$cor_Ph2_IDI) regtable$Z_cor_HeN_Hab<-atanh(regtable$cor_He_Hab) regtable$Z_cor_HeN_IDI<-atanh(regtable$cor_He_IDI) write.table(regtable,"gen_v_env3.csv",sep=",") #Used this summary table for ANOVA and subsequent landscape-level plotting and analysis ###################################################################