# Code to assess change in the effect of adding # larvae & embryos to adult samples ## Load packages library(adegenet) library(hierfstat) library(reshape) library(ecodist) library(plyr) # Read in data files dist <- read.csv("Pond_Dist.csv")[,3] a <- read.genepop("GenePop/adult_nosibs.gen", ncode = 3) e <- read.genepop("GenePop/embryo_nosibs.gen", ncode = 3) l <- read.genepop("GenePop/larvae_nosibs.gen", ncode = 3) a <- genind2df(a) e <- genind2df(e) l <- genind2df(l) for(i in 1:ncol(a)){ a[,i] <- as.numeric(a[,i]) e[,i] <- as.numeric(e[,i]) l[,i] <- as.numeric(l[,i]) } e.l <- rbind(e,l) # Calculate Fst set.seed <- 12345 dat1 <- df2genind(a[,2:16],pop=as.factor(a$pop), ncode = 3) Fst.adult <- pairwise.fst(dat1, res.type = "matrix") Fst.adult <- lower(Fst.adult) (Mantelr.adult_Fst <- mantel(Fst.adult~dist, nperm = 100000)) # Caclualte Dc dat <- genind2genpop(dat1,other.action = "NA", quiet = T) dat <- as(dat,'genpop') Dc <- lower(dist.genpop(dat,method = 2,upper = T)) (Mantelr.adult_Dc <- mantel(Dc~dist,nperm=100000)) iters <- 1000 msat.sub <- c(5, 10, 15) LENGTH <- length(seq(from = 0, to = 1,by = .05)) MixedList <- vector(mode = 'list', length = LENGTH) names(MixedList) <- seq(from = 0, to = 1,by = .05) Ho <- Ar <- vector(mode = 'list', length = iters) MixedMantelp_Fst <- vector(mode = 'list', length = iters) MixedMantelr_Fst <- vector(mode = 'list', length = iters) MixedMantelp_Dc <- vector(mode = 'list', length = iters) MixedMantelr_Dc <- vector(mode = 'list', length = iters) FinalMantelr_Fst <- FinalMatenlp_Fst <- FinalMantelr_Dc <- FinalMantelp_Dc <- vector(mode = 'list', length = length(msat.sub)) progress_bar <- create_progress_bar("text") progress_bar$init(LENGTH*iters) COUNTER <- 0 for(h in 1:length(msat.sub)) { for(i in 1:LENGTH) { i.prop <- seq(from = 0, to = 1,by = .05)[i] COUNTER <- COUNTER + 1 for(j in 1:iters) { # Subsample msats msats <- sample(2:16, msat.sub[h], replace = F) ## Subsample n.mix <- round(18 * i.prop) a.sub <- subsampind(dat=a,sampsize=18-n.mix)[,c(1,msats)] e.l.sub <- subsampind(dat=e.l,sampsize=n.mix)[,c(1,msats)] dat <- rbind(a.sub,e.l.sub) ## Calculate allelic richness Ar[j] <- mean(allelic.richness(data=dat)$Ar) ## Create genind object dat <- df2genind(dat[,-1],pop=as.factor(dat$pop), ncode = 3) ## Calculate Heterozygosity Ho[j] <- mean(summary(dat, verbose=F)$Hobs) # Calculate Fst Fst <- pairwise.fst(dat, res.type = "matrix") Fst <- lower(Fst) RESULTS <- mantel(Fst~dist,data=dat,nperm=10000)[c(1,2)] ## Mantel MixedMantelr_Fst[j] <- RESULTS[1] MixedMantelp_Fst[j] <- RESULTS[2] # Make genind object dat1 <- df2genind(dat[,2:16],pop=as.factor(dat$pop), ncode = 3) dat <- genind2genpop(dat1,other.action = "NA", quiet = T) dat <- as(dat,'genpop') Dc <- lower(dist.genpop(dat,method = 2,upper = T)) RESULTS <- mantel(Dc~dist,data=dat,nperm=10000)[c(1,2)] ## Mantel MixedMantelr_Dc[j] <- RESULTS[1] MixedMantelp_Dc[j] <- RESULTS[2] progress_bar$step() } # CLose iteration loop Fst.r <- unlist(MixedMantelr_Fst) Fst.p <- unlist(MixedMantelp_Fst) Dc.r <- unlist(MixedMantelr_Dc) Dc.p <- unlist(MixedMantelp_Dc) Ho2 <- unlist(Ho) Ar2 <- unlist(Ar) MixedList[[COUNTER]] <- data.frame(Fst.r = Fst.r, Fst.p = Fst.p, Dc.r = Dc.r, Dc.p = Dc.p, Ar = Ar2, Ho = Ho2, msat = rep(msat.sub[h], iters), ratio = rep(i.prop, iters)) } # Close mixture sample loop } # Close msat sample loop composite.df <- rbind.fill(MixedList) detach("package:plyr", unload=TRUE) library(dplyr) Fst.r_summary <- composite.df %>% group_by(msat, ratio) %>% summarise(n=n(), mn=mean(Fst.r), sd=sd(Fst.r)) %>% mutate(se=sd/sqrt(n), LCI=mn+qnorm(0.025)*se, UCI=mn+qnorm(0.975)*se) %>% as.data.frame() Fst.p_summary <- composite.df %>% group_by(msat, ratio) %>% summarise(n=n(), mn=mean(Fst.p), sd=sd(Fst.p)) %>% mutate(se=sd/sqrt(n), LCI=mn+qnorm(0.025)*se, UCI=mn+qnorm(0.975)*se) %>% as.data.frame() Dc.r_summary <- composite.df %>% group_by(msat, ratio) %>% summarise(n=n(), mn=mean(Dc.r), sd=sd(Dc.r)) %>% mutate(se=sd/sqrt(n), LCI=mn+qnorm(0.025)*se, UCI=mn+qnorm(0.975)*se) %>% as.data.frame() Dc.p_summary <- composite.df %>% group_by(msat, ratio) %>% summarise(n=n(), mn=mean(Dc.p), sd=sd(Dc.p)) %>% mutate(se=sd/sqrt(n), LCI=mn+qnorm(0.025)*se, UCI=mn+qnorm(0.975)*se) %>% as.data.frame() Ar_summary <- composite.df %>% group_by(msat, ratio) %>% summarise(n=n(), mn=mean(Ar), sd=sd(Ar)) %>% mutate(se=sd/sqrt(n), LCI=mn+qnorm(0.025)*se, UCI=mn+qnorm(0.975)*se) %>% as.data.frame() Ho_summary <- composite.df %>% group_by(msat, ratio) %>% summarise(n=n(), mn=mean(Ho), sd=sd(Ho)) %>% mutate(se=sd/sqrt(n), LCI=mn+qnorm(0.025)*se, UCI=mn+qnorm(0.975)*se) %>% as.data.frame()