#Trying changepoint install.packages("changepoint") library(changepoint) mydata = read.csv("SNPwindows.csv") #cpts() meanMBIC = cpt.mean(mydata$fracrare, method = 'PELT', penalty = 'MBIC') meanAIC = cpt.mean(mydata$fracrare, method = 'PELT', penalty = 'AIC') meanBIC = cpt.mean(mydata$fracrare, method = 'PELT', penalty = 'BIC') bothMBIC = cpt.meanvar(mydata$fracrare, method = 'PELT', penalty = 'MBIC') bothAIC = cpt.meanvar(mydata$fracrare, method = 'PELT', penalty = 'AIC') bothBIC = cpt.meanvar(mydata$fracrare, method = 'PELT', penalty = 'BIC') length(cpts(meanMBIC)) length(cpts(meanAIC)) length(cpts(meanBIC)) length(cpts(bothMBIC)) length(cpts(bothAIC)) length(cpts(bothBIC)) write.csv(cpts(bothMBIC), "MBIC.csv") write.csv(cpts(bothAIC), "AIC.csv") write.csv(cpts(bothBIC), "BIC_500snpwindows.csv") write.csv(param.est(bothMBIC), "MBICmeans.csv") write.csv(param.est(bothAIC), "AICmeans.csv") write.csv(param.est(bothBIC), "BICmeans_500snpwindows.csv") BIC = read.csv("BIC.csv") BICmeans = read.csv("BICmeans.csv") plot(BIC$lengths, BICmeans$mean) lines(BICmeans$mean[1:100]) #500 SNP windows mydata = read.delim("500snp.simple.txt") mydata = read.csv("GH_new.csv") head(mydata) try1 = lm(AvgLogMass ~ LoadIn5New + LoadOut5New + Set, data = mydata) try2 = lm(AvgLogMass ~ LoadIn7.5New + LoadOut7.5New + Set, data = mydata) try3 = lm(AvgLogMass ~ LoadIn10New + LoadOut10New + Set, data = mydata) try4 = lm(AvgLogMass ~ NewTot + Set, data = mydata) summary(try1) summary(try2) summary(try3) summary(try4) anova(try4) summary(try5) summary(try6) summary(try7) try8 = lm(AvgLogMass_1 ~ LoadIn7.5New_1, data = mydata) summary(try8) anova(try8) try2 = lm(AvgLogMass ~ LoadIn5New + Set, data = mydata) try4 = lm(AvgLogMass_1 ~ LoadIn7.5New_1 + Set_1, data = mydata) try6 = lm(AvgLogMass ~ LoadIn10New + Set, data = mydata) try3 = lm(AvgLogMass ~ LoadOut5New + Set, data = mydata) try5 = lm(AvgLogMass ~ LoadOut7.5New + Set, data = mydata) try7 = lm(AvgLogMass ~ LoadOut10New + Set, data = mydata) try1 = lm(AvgLogMass ~ NewTot + Set, data = mydata) #separated out GHwZ = lm(AvgLogMass ~ LoadIn5 + Set, data = mydata) GHwoutZ = lm(AvgLogMass ~ LoadIn5, data = mydata) summary(GHwZ) anova(GHwZ) summary(GHwoutZ) anova(GHwoutZ) #W/out Z12 GHwZ = lm(AvgLogMass_1 ~ LoadIn5_1 + Set_1, data = mydata) GHwoutZ = lm(AvgLogMass_1 ~ LoadIn5_1, data = mydata) summary(GHwZ) anova(GHwZ) summary(GHwoutZ) anova(GHwoutZ) mydata = read.csv("Field_new.csv") head(mydata) library(rsq) surv1 <- glm(Surv ~ NewTot, binomial, data = mydata) surv2 <- glm(Surv ~ LoadIn5New, binomial, data = mydata) surv3 <- glm(Surv ~ LoadOut5New, binomial, data = mydata) surv4 <- glm(Surv ~ LoadIn7.5New, binomial, data = mydata) surv5 <- glm(Surv ~ LoadOut7.5New, binomial, data = mydata) surv6 <- glm(Surv ~ LoadIn10New, binomial, data = mydata) surv7 <- glm(Surv ~ LoadOut10New, binomial, data = mydata) SS1 <- glm(SS ~ NewTot, poisson, data = mydata) SS2 <- glm(SS ~ LoadIn5New, poisson, data = mydata) SS3 <- glm(SS ~ LoadOut5New, poisson, data = mydata) SS4 <- glm(SS ~ LoadIn7.5New, poisson, data = mydata) SS5 <- glm(SS ~ LoadOut7.5New, poisson, data = mydata) SS6 <- glm(SS ~ LoadIn10New, poisson, data = mydata) SS7 <- glm(SS ~ LoadOut10New, poisson, data = mydata) summary(surv1) rsq(surv1) summary(surv2) rsq(surv2) summary(surv3) rsq(surv3) summary(surv4) rsq(surv4) summary(surv5) rsq(surv5) summary(surv6) rsq(surv6) summary(surv7) rsq(surv7) summary(SS1) rsq(SS1) summary(SS2) rsq(SS2) summary(SS3) rsq(SS3) summary(SS4) rsq(SS4) summary(SS5) rsq(SS5) summary(SS6) rsq(SS6) summary(SS7) rsq(SS7) install.packages("lme4") library(lme4) fielddata <- read.csv("Field_new.csv") head(fielddata) glmm_5 <- glmer(Surv ~ LoadIn5New + (1|LineA), data = mydata, family = binomial) glmm_tot <- glmer(Surv ~ NewTot + (1|LineA), data = mydata, family = binomial) glmm_old <- glmer(Surv ~ fracminor + (1|LineA), data = mydata, family = binomial) summary(glmm_5) anova(glmm_5) summary(glmm_tot) anova(glmm_tot) summary(glmm_old) anova(glmm_old) install.packages("aster") library(aster) fieldfit <- read.csv('forAster_New.csv') pred <- c(0,1,2) fam <- c(1,1,3) reaster_fracminor <- reaster(resp ~ varb + fit:fracminor, random = list(line = ~ 0 + fit : LineA), pred, fam, varb, id, root, data = fieldfit) reaster_tot <- reaster(resp ~ varb + fit:tot, random = list(line = ~ 0 + fit : LineA), pred, fam, varb, id, root, data = fieldfit) reaster_loadin5 <- reaster(resp ~ varb + fit:LoadIn5New, random = list(line = ~ 0 + fit : LineA), pred, fam, varb, id, root, data = fieldfit) summary(reaster_fracminor) summary(reaster_tot) summary(reaster_loadin5) aster_fracminor <- aster(resp ~ varb + fit:fracminor, pred, fam, varb, id, root, data = fieldfit) aster_loadin5 <- aster(resp ~ varb + fit:LoadIn5New, pred, fam, varb, id, root, data = fieldfit) summary(aster_fracminor) summary(aster_loadin5) library(ggplot2) #graphing for just survival ##Graphing for Figure 2 mydata = read.csv("Field_new.csv") #individual data #The following simple models were used to get intercept and effect estimates for graphing summary(glm(Surv ~ LoadIn5New, data = mydata, family = binomial)) summary(glm(SS ~ LoadIn5New, data = mydata, family = poisson)) Fielddata = read.csv("FieldAvg_new.csv") #line averages Curves = read.csv("Curves_new.csv") #math ggplot() + geom_point(data = Fielddata, aes(x=LoadIn5New, y=Surv), size = 2, color = scales::hue_pal()(2)[2])+ geom_line(data = Curves, aes(x= X, y = NewSurv), size = 1.2, color = scales::hue_pal()(2)[1]) + theme_bw(base_size= 15) ggplot() + geom_point(data = Fielddata, aes(x=LoadIn5New, y=SS), size = 2, color = scales::hue_pal()(2)[2])+ geom_line(data = Curves, aes(x= X, y = NewSS), size = 1.2, color = scales::hue_pal()(2)[1]) + theme_bw(base_size= 15) ggplot() + geom_point(data = Fielddata, aes(x=fracrare, y=Surv), size = 2, color = scales::hue_pal()(2)[2])+ geom_line(data = Curves, aes(x= X, y = FracrareSurv), size = 1.2, color = scales::hue_pal()(2)[1]) + theme_bw(base_size= 15) ggplot() + geom_point(data = Fielddata, aes(x=fracrare, y=SS), size = 2, color = scales::hue_pal()(2)[2])+ geom_line(data = Curves, aes(x= X, y = FracrareSS), size = 1.2, color = scales::hue_pal()(2)[1]) + theme_bw(base_size= 15) mydata = read.delim("Scaff2_256_624_645.txt") scaff2_test = cpt.meanvar(mydata$fracrare, method = 'PELT', penalty = 'BIC') plot(scaff2_test) data_256 = read.delim("Scaff2_256.txt") scaff2_256 = cpt.meanvar(data_256$fracrare, method = 'PELT', penalty = 'BIC') plot(scaff2_256) data_624 = read.delim("Scaff2_624.txt") scaff2_624 = cpt.meanvar(data_624$fracrare, method = 'PELT', penalty = 'BIC') plot(scaff2_624) data_645 = read.delim("Scaff2_645.txt") scaff2_645 = cpt.meanvar(data_645$fracrare, method = 'PELT', penalty = 'BIC') plot(scaff2_645) cpts(scaff2_645) param.est(scaff2_645) plot(scaff2_256, xlab = "Window", ylab = "Fraction rare allele")