Modified from hybridViability.Rmd (Oct 2017)
Test for trade-offs between male and female function in F2 hybrids of two sympatric moth-pollinated plant species, Schiedea kaalae and S. hookeri.
setwd("~/MyDocs/MEGA/UCI/Schiedea/Analysis/HybridSeeds")
pol <- read.table("f2pollenEH.csv", header=T, sep="\t", colClasses=c(date="Date"))
pol <- pol[!is.na(pol$V1) & pol$dadfullcross!="9999",]
pol$v <- rowSums(pol[,grepl("V", colnames(pol))])/10 # average viable counts
pol$i <- rowSums(pol[,grepl("I", colnames(pol))])/10 # average inviable counts
pol$vpf <- (pol$v * 444.4/20) * 10 # number viable per flower (5 flrs)
pol$ipf <- (pol$i * 444.4/20) * 10 # number inviable per flower (5 flrs)
pol$vp <- pol$vpf / (pol$vpf + pol$ipf) # proportion viable
pol$tpf <- pol$vpf + pol$ipf
pol$t <- 10*(pol$v + pol$i)
polm2 <- pol[grepl("x",pol$momfullcross),]
polm2.split <- cbind(polm2$crossid, as.data.frame(matrix(unlist(matrix(lapply(matrix(unlist(strsplit(as.character(polm2$momfullcross), " x ", fixed=T)), nrow=2), function(x) { spl <- unlist(strsplit(unlist(strsplit(x, " ", fixed=T)), "-", fixed=T)); c(spl[1],gsub("-NA","",paste(spl[2],spl[3],spl[4], sep="-"))) }), ncol=2)), ncol=4, byrow=T)))
colnames(polm2.split) <- c("crossid","mommompop","mommomid","momdadpop","momdadid")
pold2 <- pol[grepl("x",pol$dadfullcross),]
pold2.split <- cbind(pold2$crossid, as.data.frame(matrix(unlist(matrix(lapply(matrix(unlist(strsplit(as.character(pold2$dadfullcross), " x ", fixed=T)), nrow=2), function(x) { spl <- unlist(strsplit(unlist(strsplit(x, " ", fixed=T)), "-", fixed=T)); c(spl[1],gsub("-NA","",paste(spl[2],spl[3],spl[4], sep="-"))) }), ncol=2)), ncol=4, byrow=T)))
colnames(pold2.split) <- c("crossid","dadmompop","dadmomid","daddadpop","daddadid")
polm1 <- pol[!grepl("x",pol$momfullcross),]
polm1.split <- cbind(polm1$crossid, as.data.frame(matrix(unlist(matrix(lapply(as.character(polm1$momfullcross), function(x) { spl <- unlist(strsplit(unlist(strsplit(x, " ", fixed=T)), "-", fixed=T)); c(spl[1],gsub("-NA","",paste(spl[2],spl[3],spl[4], sep="-"))) }), ncol=1)), ncol=2, byrow=T)))
colnames(polm1.split) <- c("crossid","mompop","momid")
pold1 <- pol[!grepl("x",pol$dadfullcross),]
pold1.split <- cbind(pold1$crossid, as.data.frame(matrix(unlist(matrix(lapply(as.character(pold1$dadfullcross), function(x) { spl <- unlist(strsplit(unlist(strsplit(x, " ", fixed=T)), "-", fixed=T)); c(spl[1],gsub("-NA","",paste(spl[2],spl[3],spl[4], sep="-"))) }), ncol=1)), ncol=2, byrow=T)))
colnames(pold1.split) <- c("crossid","dadpop","dadid")
pol <- merge(pol, polm2.split, by="crossid", all.x=T)
pol <- merge(pol, pold2.split, by="crossid", all.x=T)
pol <- merge(pol, polm1.split, by="crossid", all.x=T)
pol <- merge(pol, pold1.split, by="crossid", all.x=T)
popcols <- c("mompop","dadpop","mommompop","momdadpop","dadmompop","daddadpop")
pol[popcols] <- lapply(pol[popcols], function(x) { sapply(x, mapvalues, from = c("794","866","899","879","892","904","3587"), to = c("WK","WK","WK","879WKG","892WKG","904WPG","3587WP"), warn_missing=F) })
kaalpops <- c("904WPG","3587WP", "892WKG")
pop2sp <- function(x) { factor(ifelse(x %in% kaalpops, "kaal", ifelse(is.na(x),NA,"hook")))}
pol$mommomspecies <- pop2sp(pol$mommompop)
pol$momdadspecies <- pop2sp(pol$momdadpop)
pol$dadmomspecies <- pop2sp(pol$dadmompop)
pol$daddadspecies <- pop2sp(pol$daddadpop)
pol$momcross <- factor(toupper(paste0(substr(pol$mommomspecies,0,1), substr(pol$momdadspecies,0,1))))
pol$dadcross <- factor(toupper(paste0(substr(pol$dadmomspecies,0,1), substr(pol$daddadspecies,0,1))))
pol$momcross <- factor(ifelse(pol$momcross=="NANA", toupper(substr(as.character(pop2sp(pol$mompop)),0,1)), as.character(pol$momcross)))
pol$dadcross <- factor(ifelse(pol$dadcross=="NANA", toupper(substr(as.character(pop2sp(pol$dadpop)),0,1)), as.character(pol$dadcross)))
# is.na(pol$dadcross) <- is.na(pol$momcross) <- pol$crosstype=="control"
pol$dadcross <- factor(pol$dadcross, c("H", "HK", "KH", "K"))
pol$momcross <- factor(pol$momcross, c("H", "HK", "KH", "K"))
pol$crosstype <- with(pol, factor(paste(momcross, dadcross, sep="x")))
pol$species <- factor(ifelse(pol$mompop %in% kaalpops, "kaal", "hook"))
pol$dadsp <- factor(ifelse(pol$dadpop %in% kaalpops, "kaal", "hook"))
pol <- within(pol, smompop <- as.factor(paste(species,mompop,sep="")))
pol <- within(pol, mompid <- as.factor(paste(mompop,momid,sep=".")))
pol <- within(pol, dadpid <- as.factor(paste(dadpop,dadid,sep=".")))
pol$mompop <- as.character(pol$mompop)
pol$mompop[is.na(pol$mompop)] <- as.character(pol$momcross[is.na(pol$mompop)])
pol$mompop <- as.factor(pol$mompop)
pol$dadpop <- as.character(pol$dadpop)
pol$dadpop[is.na(pol$dadpop)] <- as.character(pol$dadcross[is.na(pol$dadpop)])
pol$dadpop <- as.factor(pol$dadpop)
setwd("~/MyDocs/MEGA/UCI/Schiedea/Analysis/HybridSeeds")
sds <- read.table("f2seedsEH.csv", header=T, sep="\t", colClasses=c(date_emasculation="Date", date_pollination="Date"))
msds <- sds %>%
group_by(crossid,pollen_donor) %>%
dplyr::summarize(viable_seeds = mean(viable_seeds), inviable_seeds = mean(inviable_seeds)) %>%
mutate(pviable_seeds = viable_seeds/(viable_seeds+inviable_seeds))
sdspol <- merge(msds,pol,by="crossid")
str(sdspol)
'data.frame': 43 obs. of 59 variables:
$ crossid : Factor w/ 45 levels "101-1","111-1",..: 1 2 3 4 5 6 7 8 8 10 ...
$ pollen_donor : Factor w/ 2 levels "H","K": 1 2 1 2 1 1 2 2 1 2 ...
$ viable_seeds : num 12.29 1.2 1.67 1.75 7.33 ...
$ inviable_seeds: num 0.571 1.8 1.333 0 1.667 ...
$ pviable_seeds : num 0.956 0.4 0.556 1 0.815 ...
$ date : Date, format: "2018-07-23" "2018-08-22" ...
$ momfullcross : Factor w/ 25 levels "3587-14","3587-7 x 879 10-1 ID 1",..: 14 15 16 16 17 18 18 19 19 20 ...
$ dadfullcross : Factor w/ 29 levels "3587-7 x 879 G-2 ID 4",..: 4 9 15 20 11 21 21 19 19 18 ...
$ V1 : int 32 26 19 32 37 5 14 38 38 15 ...
$ I1 : int 9 7 1 10 19 25 14 18 18 12 ...
$ V2 : int 75 40 25 38 27 3 18 37 37 13 ...
$ I2 : int 12 5 5 17 19 25 15 26 26 10 ...
$ V3 : int 67 37 22 31 25 8 10 44 44 21 ...
$ I3 : int 7 9 0 8 11 18 13 29 29 17 ...
$ V4 : int 46 32 32 38 81 7 14 44 44 13 ...
$ I4 : int 11 6 4 6 53 26 10 22 22 16 ...
$ V5 : int 66 29 14 39 29 3 4 30 30 15 ...
$ I5 : int 7 3 9 6 20 16 6 19 19 11 ...
$ V6 : int 74 33 15 33 21 9 11 46 46 17 ...
$ I6 : int 8 8 2 7 11 21 9 21 21 17 ...
$ V7 : int 62 26 12 38 20 6 12 38 38 18 ...
$ I7 : int 8 9 3 6 10 13 10 27 27 17 ...
$ V8 : int 52 31 26 30 30 4 11 37 37 16 ...
$ I8 : int 11 10 5 10 20 14 7 10 10 13 ...
$ V9 : int 61 30 23 32 28 8 12 49 49 16 ...
$ I9 : int 10 5 6 6 12 20 6 17 17 16 ...
$ V10 : int 53 29 22 48 24 8 9 35 35 19 ...
$ I10 : int 6 7 2 9 24 21 7 16 16 18 ...
$ v : num 58.8 31.3 21 35.9 32.2 6.1 11.5 39.8 39.8 16.3 ...
$ i : num 8.9 6.9 3.7 8.5 19.9 19.9 9.7 20.5 20.5 14.7 ...
$ vpf : num 13065 6955 4666 7977 7155 ...
$ ipf : num 1978 1533 822 1889 4422 ...
$ vp : num 0.869 0.819 0.85 0.809 0.618 ...
$ tpf : num 15043 8488 5488 9866 11577 ...
$ t : num 677 382 247 444 521 260 212 603 603 310 ...
$ mommompop : Factor w/ 4 levels "3587WP","879WKG",..: 2 NA 3 3 3 3 3 3 3 NA ...
$ mommomid : Factor w/ 8 levels "1","10-1","2-2",..: 8 NA 1 1 4 4 4 5 5 NA ...
$ momdadpop : Factor w/ 3 levels "879WKG","892WKG",..: 3 NA 1 1 1 1 1 1 1 NA ...
$ momdadid : Factor w/ 10 levels "1-ID-1","10-1-ID",..: 7 NA 4 4 2 4 4 9 9 NA ...
$ dadmompop : Factor w/ 4 levels "3587WP","879WKG",..: NA 2 NA NA NA 3 3 3 3 2 ...
$ dadmomid : Factor w/ 10 levels "1","10","10-1",..: NA 3 NA NA NA 2 2 1 1 10 ...
$ daddadpop : Factor w/ 3 levels "879WKG","892WKG",..: NA 3 NA NA NA 1 1 1 1 3 ...
$ daddadid : Factor w/ 15 levels "10-1-ID","10-ID-7",..: NA 12 NA NA NA 1 1 3 3 11 ...
$ mompop : Factor w/ 6 levels "3587WP","879WKG",..: 5 3 6 6 6 6 6 6 6 4 ...
$ momid : Factor w/ 10 levels "1","10-1","10-1-C1",..: NA 1 NA NA NA NA NA NA NA 6 ...
$ dadpop : Factor w/ 7 levels "3587WP","879WKG",..: 7 5 2 3 2 6 6 6 6 5 ...
$ dadid : Factor w/ 10 levels "10-1","10-4/3/13",..: 4 NA 10 2 6 NA NA NA NA NA ...
$ mommomspecies : Factor w/ 2 levels "hook","kaal": 1 NA 2 2 2 2 2 2 2 NA ...
$ momdadspecies : Factor w/ 2 levels "hook","kaal": 2 NA 1 1 1 1 1 1 1 NA ...
$ dadmomspecies : Factor w/ 2 levels "hook","kaal": NA 1 NA NA NA 2 2 2 2 1 ...
$ daddadspecies : Factor w/ 2 levels "hook","kaal": NA 2 NA NA NA 1 1 1 1 2 ...
$ momcross : Factor w/ 4 levels "H","HK","KH",..: 2 4 3 3 3 3 3 3 3 4 ...
$ dadcross : Factor w/ 4 levels "H","HK","KH",..: 1 2 1 4 1 3 3 3 3 2 ...
$ crosstype : Factor w/ 10 levels "HKxH","HKxHK",..: 1 9 6 7 6 8 8 8 8 9 ...
$ species : Factor w/ 2 levels "hook","kaal": 1 2 1 1 1 1 1 1 1 2 ...
$ dadsp : Factor w/ 2 levels "hook","kaal": 1 1 1 2 1 1 1 1 1 1 ...
$ smompop : Factor w/ 5 levels "hook879WKG","hookNA",..: 2 4 2 2 2 2 2 2 2 5 ...
$ mompid : Factor w/ 11 levels "3587WP.14","879WKG.10-1",..: 11 8 11 11 11 11 11 11 11 9 ...
$ dadpid : Factor w/ 11 levels "3587WP.A","879WKG.10-1",..: 11 10 5 6 4 10 10 10 10 10 ...
Take out F2s before analysis
table(sdspol$crosstype)
HKxH HKxHK HKxK HxHK HxKH KHxH KHxK KHxKH KxHK KxKH
3 9 2 5 3 3 4 8 3 3
sdspol <- sdspol %>% filter(!(crosstype %in% c("HKxHK","KHxKH")))
(sdspol.plot <- ggplot(sdspol, aes(x=vpf, y=viable_seeds, color=crosstype, linetype=pollen_donor, shape=pollen_donor)) +
geom_point(size=3) +
#geom_smooth(method = "lm", se=F, aes(group=crosstype)) +
# geom_quantile(aes(group=1))+
xlab("Viable pollen grains per flower") +
ylab("Mean number of viable seeds") +
scale_color_brewer("Backcross Type", type="qual", palette="Paired") +
#scale_shape_manual("Maternal Plant", values=c(15, 18, 16, 17))+
scale_shape_manual("Pollen donor", values=c(15, 19))+
theme_minimal()+
theme(legend.key.size = unit(1.7,"lines")))
#ggsave("viableseeds_viablepollen.png", sdspol.plot, type="cairo", width=6, height=6)
Variables: viable pollen grains per flower, mean number of viable seeds.
sdspol.mean <- sdspol %>% group_by(crosstype) %>% summarize_at(c("viable_seeds", "vpf"), mean)
sdspol.mean.mod <- lm(viable_seeds ~ vpf, data=sdspol.mean)
summary(aov(sdspol.mean.mod))
Df Sum Sq Mean Sq F value Pr(>F)
vpf 1 0.148 0.1476 0.057 0.819
Residuals 6 15.459 2.5765
sdspol.overall.mod <- lm(viable_seeds ~ vpf, data=sdspol)
summary(aov(sdspol.overall.mod))
Df Sum Sq Mean Sq F value Pr(>F)
vpf 1 45.67 45.67 4.578 0.0428 *
Residuals 24 239.43 9.98
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sdspol.mod <- lm(viable_seeds ~ vpf * crosstype, data=sdspol)
summary(aov(sdspol.mod))
Df Sum Sq Mean Sq F value Pr(>F)
vpf 1 45.67 45.67 12.514 0.00538 **
crosstype 7 101.11 14.44 3.957 0.02485 *
vpf:crosstype 7 101.83 14.55 3.986 0.02429 *
Residuals 10 36.50 3.65
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(sdspol.mod)
Call:
lm(formula = viable_seeds ~ vpf * crosstype, data = sdspol)
Residuals:
Min 1Q Median 3Q Max
-3.01790 -0.88314 0.08254 0.89785 2.02652
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.458e+01 4.083e+00 -3.572 0.005076 **
vpf 2.029e-03 4.162e-04 4.875 0.000647 ***
crosstypeHKxK 1.506e+01 4.638e+00 3.247 0.008760 **
crosstypeHxHK 1.588e+01 5.373e+00 2.956 0.014375 *
crosstypeHxKH 8.514e+00 4.753e+00 1.791 0.103539
crosstypeKHxH 1.243e+01 5.706e+00 2.178 0.054385 .
crosstypeKHxK 1.487e+01 4.976e+00 2.988 0.013622 *
crosstypeKxHK 1.996e+01 6.037e+00 3.307 0.007921 **
crosstypeKxKH 2.228e+01 5.434e+00 4.100 0.002146 **
vpf:crosstypeHKxK -9.577e-04 8.926e-04 -1.073 0.308498
vpf:crosstypeHxHK -2.142e-03 6.142e-04 -3.488 0.005841 **
vpf:crosstypeHxKH -3.249e-04 5.578e-04 -0.583 0.573112
vpf:crosstypeKHxH -7.853e-04 8.528e-04 -0.921 0.378769
vpf:crosstypeKHxK -1.876e-03 6.830e-04 -2.746 0.020613 *
vpf:crosstypeKxHK -2.477e-03 9.109e-04 -2.719 0.021592 *
vpf:crosstypeKxKH -2.888e-03 8.356e-04 -3.457 0.006159 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.91 on 10 degrees of freedom
Multiple R-squared: 0.872, Adjusted R-squared: 0.68
F-statistic: 4.541 on 15 and 10 DF, p-value: 0.01014