TAGfile <- read.csv("RADtags_cydmelhec_fstset2_haplo.csv") hec <- c(3:6) agl <- c(7:10) ama <- c(11:14) timE <- c(15:18) timP <- c(19:22) cord <- c(23:26) hue <- c(27:30) melC <- c(31:34) chio <- c(35:38) cyd <- c(39:42) wey <- c(43:46) Ecu <- c(47:50) melG <- c(51:54) melP <- c(55:58) ros <- c(59:62) sspnames <- c("hec", "agl", "ama", "timE", "timP", "cord", "heu", "melC", "chio", "cyd", "wey", "Ecu", "melG", "melP", "ros") sspcols <- cbind(hec, agl, ama, timE, timP, cord, hue, melC, chio, cyd, wey, Ecu, melG, melP, ros) popnames <- vector() pop <- 0 for (z in 1:14) { for (j in 1:(15-z)) { k <- j+z pop <- pop + 1 popnames[pop] <- paste("fst_",sspnames[z],"_",sspnames[k], sep ="") } } fst.calc <- function(ssp1,ssp2){ ssp1_A_freq <- (( ifelse (ssp1[,1]=="A", 1, 0)+ ifelse (ssp1[,2]=="A", 1, 0)+ ifelse (ssp1[,3]=="A", 1, 0)+ ifelse (ssp1[,4]=="A", 1, 0))/4) ssp1_C_freq <- (( ifelse (ssp1[,1]=="C", 1, 0)+ ifelse (ssp1[,2]=="C", 1, 0)+ ifelse (ssp1[,3]=="C", 1, 0)+ ifelse (ssp1[,4]=="C", 1, 0))/4) ssp1_G_freq <- (( ifelse (ssp1[,1]=="G", 1, 0)+ ifelse (ssp1[,2]=="G", 1, 0)+ ifelse (ssp1[,3]=="G", 1, 0)+ ifelse (ssp1[,4]=="G", 1, 0))/4) ssp1_T_freq <- (( ifelse (ssp1[,1]=="T", 1, 0)+ ifelse (ssp1[,2]=="T", 1, 0)+ ifelse (ssp1[,3]=="T", 1, 0)+ ifelse (ssp1[,4]=="T", 1, 0))/4) ssp1_het <- 1-(ssp1_A_freq^2 + ssp1_G_freq^2 + ssp1_C_freq^2 + ssp1_T_freq^2 ) ssp2_A_freq <- (( ifelse (ssp2[,1]=="A", 1, 0)+ ifelse (ssp2[,2]=="A", 1, 0)+ ifelse (ssp2[,3]=="A", 1, 0)+ ifelse (ssp2[,4]=="A", 1, 0))/4) ssp2_C_freq <- (( ifelse (ssp2[,1]=="C", 1, 0)+ ifelse (ssp2[,2]=="C", 1, 0)+ ifelse (ssp2[,3]=="C", 1, 0)+ ifelse (ssp2[,4]=="C", 1, 0))/4) ssp2_G_freq <- (( ifelse (ssp2[,1]=="G", 1, 0)+ ifelse (ssp2[,2]=="G", 1, 0)+ ifelse (ssp2[,3]=="G", 1, 0)+ ifelse (ssp2[,4]=="G", 1, 0))/4) ssp2_T_freq <- (( ifelse (ssp2[,1]=="T", 1, 0)+ ifelse (ssp2[,2]=="T", 1, 0)+ ifelse (ssp2[,3]=="T", 1, 0)+ ifelse (ssp2[,4]=="T", 1, 0))/4) ssp2_het <- 1-(ssp2_A_freq^2 + ssp2_G_freq^2 + ssp2_C_freq^2 + ssp2_T_freq^2 ) pop_A_freq <- (( ssp1_A_freq + ssp2_A_freq )/2) pop_G_freq <- (( ssp1_G_freq + ssp2_G_freq )/2) pop_C_freq <- (( ssp1_C_freq + ssp2_C_freq )/2) pop_T_freq <- (( ssp1_T_freq + ssp2_T_freq )/2) Exp_het <- 1-( pop_A_freq^2 + pop_G_freq^2 + pop_C_freq^2 + pop_T_freq^2 ) ssp1Het <- mean(ssp1_het) ssp2Het <- mean(ssp2_het) ExpHetM <- mean(Exp_het) fst <- (ExpHetM - ((ssp1Het + ssp2Het)/2))/ExpHetM return(fst) } i <- 1 fst.table <- matrix(ncol=105, nrow=length(TAGfile[,1]), dimnames = list(c(),c(popnames))) for (i in 1:(length(TAGfile[,1]))) { infile <- paste("RADgeno",i,".txt", sep ="") tag <- read.table(infile, header = T, as.is = T) pop <- 0 for (z in 1:14) { for (j in 1:(15-z)) { k <- j+z pop <- pop + 1 ssp1 <- tag[,sspcols[,z]] ssp2 <- tag[,sspcols[,k]] fst.table[i,pop] <- fst.calc(ssp1,ssp2) } } } Resfile <- cbind(TAGfile, fst.table) write.csv(Resfile, "RAD_Fst_cydmelhec_allpairs.csv")