Backcross Seeds GLM
John Powers
Dec 13, 2019
Modified from hybridViability.Rmd (Oct 2017)
Purpose
Test for cytonuclear incompatabilities in two sympatric moth-pollinated plant species, Schiedea kaalae and S. hookeri.
Data Import
#DATA CHANGED
#changed from within to hybrid in hybridseeds.csv and f1&hybridseeds.csv: 615 597 hook 794 2 892 hybrid K HK 16
#need to fix dadpop = closed in hybridseeds.csv
setwd("~/MyDocs/MEGA/UCI/Schiedea/Analysis/HybridSeeds")
sds <- read.table("f1&hybridseeds2.csv", header=T, sep="\t")
sds <- within(sds, yesno <- as.numeric(seednum>0))
sdsm2 <- sds[grepl("x",sds$momfullcross),]
sdsm2.split <- cbind(sdsm2$crossid, as.data.frame(matrix(unlist(matrix(lapply(matrix(unlist(strsplit(as.character(sdsm2$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(sdsm2.split) <- c("crossid","mommompop","mommomid","momdadpop","momdadid")
sdsd2 <- sds[grepl("x",sds$dadfullcross),]
sdsd2.split <- cbind(sdsd2$crossid, as.data.frame(matrix(unlist(matrix(lapply(matrix(unlist(strsplit(as.character(sdsd2$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(sdsd2.split) <- c("crossid","dadmompop","dadmomid","daddadpop","daddadid")
sdsm1 <- sds[!grepl("x",sds$momfullcross),]
sdsm1.split <- cbind(sdsm1$crossid, as.data.frame(matrix(unlist(matrix(lapply(as.character(sdsm1$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(sdsm1.split) <- c("crossid","mompop","momid")
sdsd1 <- sds[!grepl("x",sds$dadfullcross),]
sdsd1.split <- cbind(sdsd1$crossid, as.data.frame(matrix(unlist(matrix(lapply(as.character(sdsd1$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(sdsd1.split) <- c("crossid","dadpop","dadid")
sds <- merge(sds, sdsm2.split, by="crossid", all.x=T)
sds <- merge(sds, sdsd2.split, by="crossid", all.x=T)
sds <- merge(sds, sdsm1.split, by="crossid", all.x=T)
sds <- merge(sds, sdsd1.split, by="crossid", all.x=T)
pop2sp <- function(x) { factor(ifelse(x %in% c("904","3587", "892"), "kaal", ifelse(is.na(x),NA,"hook")))}
sds$mommomspecies <- pop2sp(sds$mommompop)
sds$momdadspecies <- pop2sp(sds$momdadpop)
sds$dadmomspecies <- pop2sp(sds$dadmompop)
sds$daddadspecies <- pop2sp(sds$daddadpop)
sds$momcross <- factor(toupper(paste0(substr(sds$mommomspecies,0,1), substr(sds$momdadspecies,0,1))))
sds$dadcross <- factor(toupper(paste0(substr(sds$dadmomspecies,0,1), substr(sds$daddadspecies,0,1))))
sds$momcross <- factor(ifelse(sds$momcross=="NANA", toupper(substr(as.character(pop2sp(sds$mompop)),0,1)), as.character(sds$momcross)))
sds$dadcross <- factor(ifelse(sds$dadcross=="NANA", toupper(substr(as.character(pop2sp(sds$dadpop)),0,1)), as.character(sds$dadcross)))
is.na(sds$dadcross) <- is.na(sds$momcross) <- sds$crosstype=="control"
#table(sds$crosstype)
#with(sds, plot(seednum~crosstype, las=2))
sds$dadcross <- factor(sds$dadcross, c("H", "HK", "KH", "K"))
sds$momcross <- factor(sds$momcross, c("H", "HK", "KH", "K"))
#table(sds$momcross, sds$dadcross)
#table(sds$momcross, sds$dadcross, sds$crosstype) #errors in crosstype! don't use this variable
sds <- sds[sds$crosstype!="control",]
sds <- sds[!(sds$momcross=="KH" & sds$dadcross=="HK"),]#setwd("~/MyDocs/MEGA/UCI/Schiedea/Analysis/HybridSeeds")
f2sds <- read.csv("f2seedsEH.csv", sep="\t")
f2sdsm2 <- f2sds[grepl("x",f2sds$momfullcross),]
f2sdsm2.split <- cbind(f2sdsm2$pollinationid, as.data.frame(matrix(unlist(matrix(lapply(matrix(unlist(strsplit(as.character(f2sdsm2$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(f2sdsm2.split) <- c("pollinationid","mommompop","mommomid","momdadpop","momdadid")
f2sdsd2 <- f2sds[grepl("x",f2sds$dadfullcross),]
f2sdsd2.split <- cbind(f2sdsd2$pollinationid, as.data.frame(matrix(unlist(matrix(lapply(matrix(unlist(strsplit(as.character(f2sdsd2$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(f2sdsd2.split) <- c("pollinationid","dadmompop","dadmomid","daddadpop","daddadid")
f2sdsm1 <- f2sds[!grepl("x",f2sds$momfullcross),]
f2sdsm1.split <- cbind(f2sdsm1$pollinationid, as.data.frame(matrix(unlist(matrix(lapply(as.character(f2sdsm1$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(f2sdsm1.split) <- c("pollinationid","mompop","momid")
f2sdsd1 <- f2sds[!grepl("x",f2sds$dadfullcross),]
f2sdsd1.split <- cbind(f2sdsd1$pollinationid, as.data.frame(matrix(unlist(matrix(lapply(as.character(f2sdsd1$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(f2sdsd1.split) <- c("pollinationid","dadpop","dadid")
f2sds <- merge(f2sds, f2sdsm2.split, by="pollinationid", all.x=T)
f2sds <- merge(f2sds, f2sdsd2.split, by="pollinationid", all.x=T)
f2sds <- merge(f2sds, f2sdsm1.split, by="pollinationid", all.x=T)
f2sds <- merge(f2sds, f2sdsd1.split, by="pollinationid", all.x=T)
popcols <- c("mompop","dadpop","mommompop","momdadpop","dadmompop","daddadpop")
f2sds[popcols] <- lapply(f2sds[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")))}
f2sds$mommomspecies <- pop2sp(f2sds$mommompop)
f2sds$momdadspecies <- pop2sp(f2sds$momdadpop)
f2sds$dadmomspecies <- pop2sp(f2sds$dadmompop)
f2sds$daddadspecies <- pop2sp(f2sds$daddadpop)
f2sds$momcross <- factor(toupper(paste0(substr(f2sds$mommomspecies,0,1), substr(f2sds$momdadspecies,0,1))))
f2sds$dadcross <- factor(toupper(paste0(substr(f2sds$dadmomspecies,0,1), substr(f2sds$daddadspecies,0,1))))
f2sds$momcross <- factor(ifelse(f2sds$momcross=="NANA", toupper(substr(as.character(pop2sp(f2sds$mompop)),0,1)), as.character(f2sds$momcross)))
f2sds$dadcross <- factor(ifelse(f2sds$dadcross=="NANA", toupper(substr(as.character(pop2sp(f2sds$dadpop)),0,1)), as.character(f2sds$dadcross)))
# is.na(f2sds$dadcross) <- is.na(f2sds$momcross) <- f2sds$crosstype=="control"
f2sds$dadcross <- factor(f2sds$dadcross, c("H", "HK", "KH", "K"))
f2sds$momcross <- factor(f2sds$momcross, c("H", "HK", "KH", "K"))
f2sds$crosstype <- with(f2sds, factor(paste(momcross, dadcross, sep="x")))
f2sds$species <- factor(ifelse(f2sds$mompop %in% kaalpops, "kaal", "hook"))
f2sds$dadsp <- factor(ifelse(f2sds$dadpop %in% kaalpops, "kaal", "hook"))
f2sds <- within(f2sds, smompop <- as.factor(paste(species,mompop,sep="")))
f2sds <- within(f2sds, mompid <- as.factor(paste(mompop,momid,sep=".")))
f2sds <- within(f2sds, dadpid <- as.factor(paste(dadpop,dadid,sep=".")))
f2sds$mompop <- as.character(f2sds$mompop)
f2sds$mompop[is.na(f2sds$mompop)] <- as.character(f2sds$momcross[is.na(f2sds$mompop)])
f2sds$mompop <- as.factor(f2sds$mompop)
f2sds$dadpop <- as.character(f2sds$dadpop)
f2sds$dadpop[is.na(f2sds$dadpop)] <- as.character(f2sds$dadcross[is.na(f2sds$dadpop)])
f2sds$dadpop <- as.factor(f2sds$dadpop)
# check final structure
#str(f2sds)
#the momcross for F2s is momcross x dadcross and the dadcross is pollen_donor. So put "mom" in front of everything
f2sds.fix <- f2sds %>%
rename_at(vars(matches("mom|dad|cross|species")), function(x) paste0("mom",x)) %>%
dplyr::rename(momcross = momcrosstype, dadcross= pollen_donor, seednum=viable_seeds) %>%
mutate(yesno = as.numeric(seednum > 0),
crosstype=paste0(momcross, "x", dadcross),
momfullcross =paste(mommomfullcross, momdadfullcross),
dadfullcross = dadcross,
crossid=paste0("P",pollinationid))
f2sds <- f2sds.fix
#Inventory of attempted crosses
#with(f2sds, table(momcross, dadcross))#intersect(colnames(sds), colnames(f2sds))
sds$crossid <- as.character(sds$crossid)
sds.all <- bind_rows(sds, f2sds) %>% mutate_if(is.character, as.factor)
str(sds.all) 'data.frame': 1872 obs. of 51 variables:
$ crossid : Factor w/ 1872 levels "1","100","1000",..: 1 536 1133 1205 1286 316 372 436 500 526 ...
$ momfullcross : Factor w/ 111 levels "3587-1","3587-10-1",..: 65 65 65 65 65 66 66 66 66 66 ...
$ dadfullcross : Factor w/ 102 levels "3587","3587-10",..: 20 20 20 20 20 90 32 32 89 89 ...
$ crosstype : Factor w/ 23 levels "FF","FH","FK",..: 2 2 2 2 2 3 2 2 3 3 ...
$ seednum : int 3 1 2 0 9 1 3 0 0 0 ...
$ yesno : num 1 1 1 0 1 1 1 0 0 0 ...
$ mommompop : Factor w/ 10 levels "3587","3587WP",..: 3 3 3 3 3 3 3 3 3 3 ...
$ mommomid : Factor w/ 17 levels "1","10","10-1",..: 17 17 17 17 17 17 17 17 17 17 ...
$ momdadpop : Factor w/ 10 levels "3587WP","879",..: 6 6 6 6 6 6 6 6 6 6 ...
$ momdadid : Factor w/ 33 levels "10-1","10-1ID1",..: 23 23 23 23 23 24 24 24 24 24 ...
$ dadmompop : Factor w/ 4 levels "3587","879","892",..: NA NA NA NA NA NA NA NA NA NA ...
$ dadmomid : Factor w/ 12 levels "1","10","10-1",..: NA NA NA NA NA NA NA NA NA NA ...
$ daddadpop : Factor w/ 3 levels "879","892","904": NA NA NA NA NA NA NA NA NA NA ...
$ daddadid : Factor w/ 39 levels ".202ID5","10-1ID1",..: NA NA NA NA NA NA NA NA NA NA ...
$ mompop : Factor w/ 7 levels "3587","794","879",..: NA NA NA NA NA NA NA NA NA NA ...
$ momid : Factor w/ 32 levels "1","10","10-1",..: NA NA NA NA NA NA NA NA NA NA ...
$ dadpop : Factor w/ 11 levels "3587","794","866",..: 3 3 3 3 3 7 4 4 7 7 ...
$ dadid : Factor w/ 28 levels "1","10","10-1",..: 14 14 14 14 14 18 5 5 17 17 ...
$ mommomspecies : Factor w/ 2 levels "hook","kaal": 1 1 1 1 1 1 1 1 1 1 ...
$ momdadspecies : Factor w/ 2 levels "hook","kaal": 2 2 2 2 2 2 2 2 2 2 ...
$ dadmomspecies : Factor w/ 2 levels "hook","kaal": NA NA NA NA NA NA NA NA NA NA ...
$ daddadspecies : Factor w/ 2 levels "hook","kaal": NA NA NA NA NA NA NA NA NA NA ...
$ momcross : Factor w/ 14 levels "H","HK","HKxH",..: 2 2 2 2 2 2 2 2 2 2 ...
$ dadcross : Factor w/ 4 levels "H","HK","K","KH": 1 1 1 1 1 3 1 1 3 3 ...
$ pollinationid : int NA NA NA NA NA NA NA NA NA NA ...
$ mommomfullcross : Factor w/ 29 levels "3587-14","3587-7 x 879 10-1 ID 1",..: NA NA NA NA NA NA NA NA NA NA ...
$ momdadfullcross : Factor w/ 31 levels "3587-7 x 879 G-2 ID 4",..: NA NA NA NA NA NA NA NA NA NA ...
$ momcrossid : Factor w/ 45 levels "101-1","111-1",..: NA NA NA NA NA NA NA NA NA NA ...
$ date_emasculation: Factor w/ 23 levels "","2018-10-31",..: NA NA NA NA NA NA NA NA NA NA ...
$ date_pollination : Factor w/ 21 levels "2018-11-06","2018-11-07",..: NA NA NA NA NA NA NA NA NA NA ...
$ inviable_seeds : int NA NA NA NA NA NA NA NA NA NA ...
$ comments : Factor w/ 23 levels "","Branch dying at bagging; bagged",..: NA NA NA NA NA NA NA NA NA NA ...
$ mommommompop : Factor w/ 4 levels "3587WP","879WKG",..: NA NA NA NA NA NA NA NA NA NA ...
$ mommommomid : Factor w/ 10 levels "1","10-1","2",..: NA NA NA NA NA NA NA NA NA NA ...
$ mommomdadpop : Factor w/ 3 levels "879WKG","892WKG",..: NA NA NA NA NA NA NA NA NA NA ...
$ mommomdadid : Factor w/ 11 levels "1-ID-1","10-1-ID",..: NA NA NA NA NA NA NA NA NA NA ...
$ momdadmompop : Factor w/ 4 levels "3587WP","879WKG",..: NA NA NA NA NA NA NA NA NA NA ...
$ momdadmomid : Factor w/ 10 levels "1","10","10-1",..: NA NA NA NA NA NA NA NA NA NA ...
$ momdaddadpop : Factor w/ 3 levels "879WKG","892WKG",..: NA NA NA NA NA NA NA NA NA NA ...
$ momdaddadid : Factor w/ 16 levels "10-1-ID","10-ID-7",..: NA NA NA NA NA NA NA NA NA NA ...
$ mommommomspecies : Factor w/ 2 levels "hook","kaal": NA NA NA NA NA NA NA NA NA NA ...
$ mommomdadspecies : Factor w/ 2 levels "hook","kaal": NA NA NA NA NA NA NA NA NA NA ...
$ momdadmomspecies : Factor w/ 2 levels "hook","kaal": NA NA NA NA NA NA NA NA NA NA ...
$ momdaddadspecies : Factor w/ 2 levels "hook","kaal": NA NA NA NA NA NA NA NA NA NA ...
$ mommomcross : Factor w/ 4 levels "H","HK","KH",..: NA NA NA NA NA NA NA NA NA NA ...
$ momdadcross : Factor w/ 4 levels "H","HK","KH",..: NA NA NA NA NA NA NA NA NA NA ...
$ momspecies : Factor w/ 2 levels "hook","kaal": NA NA NA NA NA NA NA NA NA NA ...
$ momdadsp : Factor w/ 2 levels "hook","kaal": NA NA NA NA NA NA NA NA NA NA ...
$ momsmompop : Factor w/ 5 levels "hook879WKG","hookNA",..: NA NA NA NA NA NA NA NA NA NA ...
$ mommompid : Factor w/ 12 levels "3587WP.14","879WKG.10-1",..: NA NA NA NA NA NA NA NA NA NA ...
$ momdadpid : Factor w/ 12 levels "3587WP.A","879WKG.10-1",..: NA NA NA NA NA NA NA NA NA NA ...sds <- sds.all
sds$dadcross <- factor(sds$dadcross, c("H", "HK", "KH", "K"))
sds$momcross <- factor(sds$momcross, c("H", "HK", "HxKH", "HxHK", "HKxH", "KHxH","HKxHK",
"KHxKH", "HKxK", "KHxK", "KxKH", "KxHK", "KH", "K" ))Inventory
Capsules
#number of capsules (multiple capsules per plant)
with(sds, table(momcross, dadcross)) dadcross
momcross H HK KH K
H 112 61 105 156
HK 117 97 0 67
HxKH 11 0 0 0
HxHK 26 0 0 0
HKxH 28 0 0 0
KHxH 25 0 0 0
HKxHK 18 0 0 20
KHxKH 12 0 0 23
HKxK 0 0 0 7
KHxK 0 0 0 22
KxKH 0 0 0 18
KxHK 0 0 0 12
KH 191 0 195 177
K 177 37 69 89sds_n <- sds %>% group_by(momcross, dadcross) %>% tally()
sds_capsule <- sds %>% group_by(momcross, dadcross) %>% summarize_at(vars(yesno), mean)Maternal plants
sds_mom_n <- sds %>% group_by(momcross, dadcross, momfullcross) %>% tally
with(sds_mom_n, table(momcross, dadcross)) dadcross
momcross H HK KH K
H 10 10 9 10
HK 13 11 0 8
HxKH 3 0 0 0
HxHK 5 0 0 0
HKxH 4 0 0 0
KHxH 4 0 0 0
HKxHK 4 0 0 3
KHxKH 5 0 0 4
HKxK 0 0 0 2
KHxK 0 0 0 5
KxKH 0 0 0 3
KxHK 0 0 0 4
KH 19 0 19 19
K 14 6 11 14Maternal plants x Paternal plants
sds_momdad_n <- sds %>% group_by(momcross, dadcross, momfullcross, dadfullcross) %>% tally
with(sds_momdad_n, table(momcross, dadcross)) dadcross
momcross H HK KH K
H 28 14 26 30
HK 20 12 0 15
HxKH 3 0 0 0
HxHK 5 0 0 0
HKxH 4 0 0 0
KHxH 4 0 0 0
HKxHK 4 0 0 3
KHxKH 5 0 0 4
HKxK 0 0 0 2
KHxK 0 0 0 5
KxKH 0 0 0 3
KxHK 0 0 0 4
KH 26 0 22 24
K 41 10 16 36#lump by plant
f2sds_plant_n <- f2sds %>% group_by(mommomfullcross, momdadfullcross, momcrossid, dadcross) %>% tally()
#hist(f2sds_plant_n$n)
#leave out crosstype (errors), crossid, pollinationid, momcrossid, dependent variables, dates, comments
sds_plant_n <- sds %>% group_by(momfullcross, dadfullcross, mommompop, mommomid, momdadpop, momdadid, dadmompop, dadmomid, daddadpop, daddadid, mompop, momid, dadpop, dadid, mommomspecies, momdadspecies, dadmomspecies, daddadspecies, momcross, dadcross, mommomfullcross, momdadfullcross, mommommompop, mommommomid, mommomdadpop, momdadmomid, momdaddadpop, momdaddadid, mommommomspecies, mommomdadspecies, momdaddadspecies, mommomcross, momdadcross, momspecies, momdadsp, momsmompop, mommompid, momdadpid) %>% tally()
hist(sds_plant_n$n)
nrow(sds_plant_n)
#these two variables (momfullcross, dadfullcross) are sufficient to describe the rest
sds %>% group_by(momfullcross, dadfullcross) %>% tally %>% nrow
sds %>% group_by(momcross, dadcross, momfullcross, dadfullcross) %>% tally
#there are some plant x plant with multiple momcrossid's - lump those
sds %>% group_by(momfullcross, dadfullcross, momcrossid) %>% tally %>% nrow
sds %>% group_by(momfullcross, momcrossid) %>% tally %>% filter(momfullcross %in% momfullcross[duplicated(momfullcross)])
#there are some plant x plant with multiple crosstypes - these are the crosstype errors but we don't use this variable so OK
sds %>% group_by(momfullcross, dadfullcross, crosstype) %>% tally %>% filter(momfullcross %in% momfullcross[duplicated(momfullcross)])
#the plant x plant crosses with the highest tallies are often missing the dad plant ID (just has dad population from Thu)Plot of raw data
(bothplot <- ggplot(sds, aes(y=seednum, x=dadcross, fill=dadcross)) +
facet_wrap(vars(momcross), nrow=1, strip.position="bottom")+
#geom_boxplot(outlier.size=0.5) +
geom_violin(scale="width")+
stat_summary(fun.y = median, fun.ymin = median, fun.ymax = median, geom = "crossbar", width = 0.9) +
geom_text(data=sds_n, aes(y=-2, label=n))+
geom_text(data=sds_capsule, aes(y=-1, label=scales::percent(yesno, accuracy=1)))+
labs(x="\nCapsule formation (%)\nNumber of pollinations\nFemale plant\nPollen donor", y="Viable seeds") +
scale_y_continuous(expand = expand_scale(add=c(0.5,0.5)), breaks = scales::pretty_breaks(n = 6)) +
scale_fill_manual("Paternal plant", values=c("grey90","gray70","gray50","gray30"))+
theme_classic()+ theme(axis.text = element_text(colour="black", size=rel(1)),
text=element_text(size=14),
axis.ticks.x = element_blank(),
panel.grid.major.y = element_line(color="grey80"))+
guides(fill=F))
#ggsave("seedplot_all.pdf", bothplot, width=24, height=10)Models
Overall model
Number of seeds per pollination
f1.md.qpoi <- glm(seednum~momcross*dadcross, data=sds, family="quasipoisson")
Anova(f1.md.qpoi) Analysis of Deviance Table (Type II tests)
Response: seednum
LR Chisq Df Pr(>Chisq)
momcross 391.03 13 < 2.2e-16 ***
dadcross 42.63 3 2.952e-09 ***
momcross:dadcross 118.10 9 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cld.mod <- f1.md.qpoi %>%
ref_grid %>% emmeans(~ momcross*dadcross) %>%
summary %>%
mutate(response = exp(emmean), uSE = exp(emmean+SE), lSE = exp(emmean-SE))
cld.mod[,-c(3:7)] momcross dadcross response uSE lSE
1 H H 10.8571429 11.560398 10.1966687
2 HK H 2.0170940 2.325937 1.7492599
3 HxKH H 3.1818182 4.606138 2.1979298
4 HxHK H 0.3846154 0.768419 0.1925109
5 HKxH H 3.6428571 4.524335 2.9331178
6 KHxH H 3.8400000 4.801097 3.0712982
7 HKxHK H 1.9444444 2.814848 1.3431861
8 KHxKH H 3.8333333 5.293190 2.7761038
9 HKxK H NA NA NA
10 KHxK H NA NA NA
11 KxKH H NA NA NA
12 KxHK H NA NA NA
13 KH H 1.9109948 2.142943 1.7041525
14 K H 4.4802260 4.842313 4.1452139
15 H HK 5.0655738 5.737204 4.4725686
16 HK HK 2.1443299 2.495723 1.8424122
17 HxKH HK NA NA NA
18 HxHK HK NA NA NA
19 HKxH HK NA NA NA
20 KHxH HK NA NA NA
21 HKxHK HK NA NA NA
22 KHxKH HK NA NA NA
23 HKxK HK NA NA NA
24 KHxK HK NA NA NA
25 KxKH HK NA NA NA
26 KxHK HK NA NA NA
27 KH HK NA NA NA
28 K HK 1.5675676 2.089450 1.1760357
29 H KH 6.4571429 7.023342 5.9365885
30 HK KH NA NA NA
31 HxKH KH NA NA NA
32 HxHK KH NA NA NA
33 HKxH KH NA NA NA
34 KHxH KH NA NA NA
35 HKxHK KH NA NA NA
36 KHxKH KH NA NA NA
37 HKxK KH NA NA NA
38 KHxK KH NA NA NA
39 KxKH KH NA NA NA
40 KxHK KH NA NA NA
41 KH KH 2.8307692 3.107135 2.5789850
42 K KH 1.1304348 1.448327 0.8823167
43 H K 4.6987179 5.094325 4.3338326
44 HK K 1.5970149 1.973316 1.2924728
45 HxKH K NA NA NA
46 HxHK K NA NA NA
47 HKxH K NA NA NA
48 KHxH K NA NA NA
49 HKxHK K 2.3500000 3.233787 1.7077499
50 KHxKH K 1.8260871 2.559540 1.3028101
51 HKxK K 2.5714286 4.307326 1.5351161
52 KHxK K 1.1363636 1.760418 0.7335318
53 KxKH K 3.8888889 5.051625 2.9937805
54 KxHK K 3.5000000 4.906053 2.4969156
55 KH K 2.0000000 2.246716 1.7803766
56 K K 5.6853933 6.266356 5.1582925Hurdle model
Fraction of pollinations that produced capsules
All data
f1.md.hr <- glm(yesno~momcross*dadcross, data=sds, family="binomial")
Anova(f1.md.hr) Analysis of Deviance Table (Type II tests)
Response: yesno
LR Chisq Df Pr(>Chisq)
momcross 103.132 13 4.091e-16 ***
dadcross 1.716 3 0.6334
momcross:dadcross 75.887 9 1.056e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cld.mod <- f1.md.hr %>%
ref_grid %>% emmeans(~ momcross*dadcross) %>%
summary%>%
mutate(response=plogis(emmean), uSE = plogis(emmean+SE), lSE = plogis(emmean-SE))
cld.mod[1:26,-c(3:7)] momcross dadcross response uSE lSE
1 H H 0.7946429 0.8302011 0.7538461
2 HK H 0.5641026 0.6092755 0.5178400
3 HxKH H 0.2727273 0.4246224 0.1600532
4 HxHK H 0.1923077 0.2814114 0.1264523
5 HKxH H 0.5000000 0.5933821 0.4066179
6 KHxH H 0.7600000 0.8349292 0.6647170
7 HKxHK H 0.2777778 0.3942991 0.1851638
8 KHxKH H 0.3333333 0.4798173 0.2132375
9 HKxK H NA NA NA
10 KHxK H NA NA NA
11 KxKH H NA NA NA
12 KxHK H NA NA NA
13 KH H 0.4869110 0.5230827 0.4508759
14 K H 0.7288136 0.7609115 0.6941386
15 H HK 0.7049180 0.7597961 0.6433869
16 HK HK 0.6391753 0.6863718 0.5891279
17 HxKH HK NA NA NA
18 HxHK HK NA NA NA
19 HKxH HK NA NA NA
20 KHxH HK NA NA NA
21 HKxHK HK NA NA NA
22 KHxKH HK NA NA NA
23 HKxK HK NA NA NA
24 KHxK HK NA NA NA
25 KxKH HK NA NA NA
26 KxHK HK NA NA NAggplot(as.data.frame(cld.mod), aes(y=response, x=momcross, fill=dadcross)) +
geom_col(position=position_dodge2()) +
geom_linerange(aes(ymin=lSE, ymax=uSE), position=position_dodge2(0.9)) +
labs(x="Maternal plant", y="Capsule formation") +
#geom_text(aes(label=trimws(.group), x=momcross, hjust=length(trimws(.group))/30), position=position_dodge2(0.9), vjust=-4) +
geom_text(aes(label=dadcross, y=0.05), position=position_dodge2(0.9)) +
scale_y_continuous(expand = expand_scale(mult=c(0,0)), breaks = scales::pretty_breaks(n = 5), limits=c(0,1)) +
scale_fill_manual(values=c("grey80","gray60","gray20","gray40"))+
theme_classic() + theme(legend.text=element_text(size=rel(1)), legend.position="bottom", axis.text = element_text(colour="black", size=rel(1)), text=element_text(size=14), axis.ticks.x = element_blank()) + guides(fill=F)
Contrasts for cytonuclear incompatibility
This analysis uses the maternal plant ID (momfullcross) as a random effect.
f1.m.hr <- glm(yesno~momcross, data=sds, family="binomial")
f1.m.mixed.hr <- glmmTMB(yesno~momcross + (1|momfullcross), data=sds, family="binomial")
K = matrix(rep(0,14*4), 4)
colnames(K) <- substr(names(coef(f1.m.hr)),9,13)
K[1,c("HKxK","KxHK")] <- K[2,c("KHxH","HxKH")] <- K[3,c("HKxH","HxHK")] <- K[4,c("KHxK","KxKH")] <- c(1, -1)
rownames(K) <- c("HK x K - K x HK","KH x H - H x KH","HK x H - H x HK","KH x K - K x KH")
library(multcomp)
summary(glht(f1.m.mixed.hr, linfct=K),test=adjusted("none"))
Simultaneous Tests for General Linear Hypotheses
Fit: glmmTMB(formula = yesno ~ momcross + (1 | momfullcross), data = sds,
family = "binomial", ziformula = ~0, dispformula = ~1)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
HK x K - K x HK == 0 -1.5923 1.6520 -0.964 0.33512
KH x H - H x KH == 0 2.4233 1.1561 2.096 0.03607 *
HK x H - H x HK == 0 1.7457 0.9252 1.887 0.05918 .
KH x K - K x KH == 0 -4.0245 1.3611 -2.957 0.00311 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- none method)cld.mod <- f1.m.mixed.hr %>%
ref_grid %>% emmeans(~ momcross) %>%
summary%>%
mutate(response=plogis(emmean), uSE = plogis(emmean+SE), lSE = plogis(emmean-SE))
cld.mod[,-c(2:6)] momcross response uSE lSE
1 H 0.7044810 0.7545873 0.64890434
2 HK 0.5664087 0.6384150 0.49148474
3 HxKH 0.2490996 0.4556974 0.11617475
4 HxHK 0.1643308 0.2794974 0.09064807
5 HKxH 0.5297901 0.6784332 0.37566787
6 KHxH 0.7891619 0.8818788 0.65235759
7 HKxHK 0.4927057 0.6412917 0.34539726
8 KHxKH 0.1599799 0.2700769 0.08927466
9 HKxK 0.7584987 0.9057453 0.50654331
10 KHxK 0.2670157 0.4127048 0.15884621
11 KxKH 0.9532337 0.9852959 0.86111424
12 KxHK 0.9391591 0.9812133 0.82021544
13 KH 0.5649653 0.6196988 0.50860089
14 K 0.6539257 0.7107177 0.59237878grn_pur <- brewer.pal(10,"Paired")[c(8,2)]
cld.mod$pair <- factor(c(1,2,3,4,4,3,5,5,6,7,7,6,2,1))[as.integer(rownames(cld.mod))] #the reciprocal pairs
contr_order <- c("HxHK","HKxH","HxKH","KHxH","KxKH","KHxK","KxHK","HKxK")
otherlevels <- levels(cld.mod$momcross)[!(levels(cld.mod$momcross) %in% contr_order)]
cld.mod$momcross2 <- factor(cld.mod$momcross, levels = c(contr_order, otherlevels))
levels(cld.mod$momcross2) <- gsub("."," x ",levels(cld.mod$momcross2), fixed=T)
(em_cont_r <- ggplot(as.data.frame(cld.mod[cld.mod$momcross2 %in% contr_order,]),
aes(y=response, x=momcross2, fill=pair)) +
geom_col() +
geom_linerange(aes(ymin=lSE, ymax=uSE)) +
labs(x="Cross", y="Capsule formation") +
scale_y_continuous(expand = expand_scale(mult=c(0,0)), labels=scales::percent, breaks = scales::pretty_breaks(n = 5), limits=c(0,1)) +
scale_fill_manual(values=setNames(rep(grn_pur,each=2), c(4, 7, 3, 6)))+ #7,4,6,1
theme_classic() + theme(legend.text=element_text(size=rel(1)), legend.position="bottom", axis.text = element_text(colour="black", size=rel(1)), text=element_text(size=14), axis.ticks.x = element_blank()) + guides(fill=F))
ggsave("hybrid_capsule_F2.png", em_cont_r, width=7, height = 0.75*7, units="in")Number of seeds if capsule formed
All data
f1.md.qpoi.tr <- glm(seednum~momcross*dadcross, data=subset(sds,seednum>0), family="quasipoisson")
Anova(f1.md.qpoi.tr) Analysis of Deviance Table (Type II tests)
Response: seednum
LR Chisq Df Pr(>Chisq)
momcross 409.99 13 < 2.2e-16 ***
dadcross 47.79 3 2.364e-10 ***
momcross:dadcross 68.47 9 3.035e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cld.mod <- f1.md.qpoi.tr %>%
ref_grid %>% emmeans(~ momcross*dadcross) %>%
summary %>%
mutate(response = exp(emmean), uSE = exp(emmean+SE), lSE = exp(emmean-SE))
cld.mod[,-c(3:7)] momcross dadcross response uSE lSE
1 H H 13.662921 14.328336 13.028409
2 HK H 3.575758 3.983336 3.209883
3 HxKH H 11.666667 15.441061 8.814881
4 HxHK H 2.000000 3.378832 1.183841
5 HKxH H 7.285714 8.585766 6.182516
6 KHxH H 5.052632 5.984382 4.265952
7 HKxHK H 7.000000 9.264636 5.288929
8 KHxKH H 11.500000 14.685224 9.005651
9 HKxK H NA NA NA
10 KHxK H NA NA NA
11 KxKH H NA NA NA
12 KxHK H NA NA NA
13 KH H 3.924731 4.280605 3.598443
14 K H 6.147287 6.520146 5.795750
15 H HK 7.186047 7.896940 6.539148
16 HK HK 3.354839 3.763624 2.990454
17 HxKH HK NA NA NA
18 HxHK HK NA NA NA
19 HKxH HK NA NA NA
20 KHxH HK NA NA NA
21 HKxHK HK NA NA NA
22 KHxKH HK NA NA NA
23 HKxK HK NA NA NA
24 KHxK HK NA NA NA
25 KxKH HK NA NA NA
26 KxHK HK NA NA NA
27 KH HK NA NA NA
28 K HK 2.761905 3.433770 2.221500
29 H KH 9.287671 9.898392 8.714631
30 HK KH NA NA NA
31 HxKH KH NA NA NA
32 HxHK KH NA NA NA
33 HKxH KH NA NA NA
34 KHxH KH NA NA NA
35 HKxHK KH NA NA NA
36 KHxKH KH NA NA NA
37 HKxK KH NA NA NA
38 KHxK KH NA NA NA
39 KxKH KH NA NA NA
40 KxHK KH NA NA NA
41 KH KH 4.246154 4.556675 3.956794
42 K KH 3.714286 4.481445 3.078453
43 H K 7.330000 7.792987 6.894520
44 HK K 3.242424 3.806194 2.762159
45 HxKH K NA NA NA
46 HxHK K NA NA NA
47 HKxH K NA NA NA
48 KHxH K NA NA NA
49 HKxHK K 3.916667 4.988428 3.075173
50 KHxKH K 10.500000 13.561677 8.129526
51 HKxK K 3.600000 5.321663 2.435329
52 KHxK K 3.571429 4.975941 2.563355
53 KxKH K 4.117647 5.020251 3.377324
54 KxHK K 3.818182 4.931512 2.956195
55 KH K 3.847826 4.202346 3.523214
56 K K 7.333333 7.894356 6.812180ggplot(as.data.frame(cld.mod), aes(y=response, x=momcross, fill=dadcross)) +
geom_col(position=position_dodge2()) +
geom_linerange(aes(ymin=lSE, ymax=uSE), position=position_dodge2(0.9)) +
labs(x="Maternal plant", y="Seeds per capsule") +
#geom_text(aes(label=trimws(.group), x=momcross, hjust=length(trimws(.group))/30), position=position_dodge2(0.9), vjust=-2) +
geom_text(aes(label=dadcross, y=0.5), position=position_dodge2(0.9)) +
scale_y_continuous(expand = expand_scale(mult=c(0,.1)), breaks = scales::pretty_breaks(n = 20))+
scale_fill_manual(values=c("grey80","gray60","gray20","gray40"))+
theme_classic() + theme(legend.text=element_text(size=rel(1)), legend.position="bottom", axis.text = element_text(colour="black", size=rel(1)), text=element_text(size=14), axis.ticks.x = element_blank()) + guides(fill=F)
Contrasts for cytonuclear incompatibility
This analysis uses the maternal plant ID (momfullcross) as a random effect.
f1.m.qpoi.tr <- glm(seednum~momcross, data=subset(sds,seednum>0), family="quasipoisson")
f1.m.mixed.tr <- glmmTMB(seednum~momcross + (1|momfullcross), data=subset(sds,seednum>0), family="nbinom1")
K = matrix(rep(0,14*4), 4)
colnames(K) <- substr(names(coef(f1.m.qpoi.tr)),9,13)
K[1,c("HKxK","KxHK")] <- K[2,c("KHxH","HxKH")] <- K[3,c("HKxH","HxHK")] <- K[4,c("KHxK","KxKH")] <- c(1, -1)
rownames(K) <- c("HK x K - K x HK","KH x H - H x KH","HK x H - H x HK","KH x K - K x KH")
summary(glht(f1.m.mixed.tr, linfct=K),test=adjusted("none"))
Simultaneous Tests for General Linear Hypotheses
Fit: glmmTMB(formula = seednum ~ momcross + (1 | momfullcross), data = subset(sds,
seednum > 0), family = "nbinom1", ziformula = ~0, dispformula = ~1)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
HK x K - K x HK == 0 -0.1346 0.5302 -0.254 0.7996
KH x H - H x KH == 0 -0.8975 0.5311 -1.690 0.0911 .
HK x H - H x HK == 0 0.7976 0.5373 1.485 0.1377
KH x K - K x KH == 0 -0.1927 0.4687 -0.411 0.6810
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- none method)cld.mod <- f1.m.mixed.tr %>%
ref_grid %>% emmeans(~ momcross) %>%
summary%>%
mutate(response = exp(emmean), uSE = exp(emmean+SE), lSE = exp(emmean-SE))
cld.mod[,-c(2:6)] momcross response uSE lSE
1 H 7.236977 8.015444 6.534115
2 HK 3.315782 3.806946 2.887988
3 HxKH 11.790326 18.791430 7.397616
4 HxHK 2.177252 3.461508 1.369468
5 HKxH 4.834133 6.361299 3.673596
6 KHxH 4.805795 6.199394 3.725472
7 HKxHK 4.005334 5.263089 3.048154
8 KHxKH 7.662486 11.027103 5.324489
9 HKxK 3.646162 5.575178 2.384587
10 KHxK 3.476069 5.015923 2.408940
11 KxKH 4.214594 5.647415 3.145298
12 KxHK 4.171463 5.725304 3.039332
13 KH 3.711861 4.112329 3.350391
14 K 5.073075 5.669737 4.539204grn_pur <- brewer.pal(10,"Paired")[c(8,2)]
cld.mod$pair <- factor(c(1,2,3,4,4,3,5,5,6,7,7,6,2,1))[as.integer(rownames(cld.mod))] #the reciprocal pairs
contr_order <- c("HxHK","HKxH","HxKH","KHxH","KxKH","KHxK","KxHK","HKxK")
otherlevels <- levels(cld.mod$momcross)[!(levels(cld.mod$momcross) %in% contr_order)]
cld.mod$momcross2 <- factor(cld.mod$momcross, levels = c(contr_order, otherlevels))
levels(cld.mod$momcross2) <- gsub("."," x ",levels(cld.mod$momcross2), fixed=T)
(em_cont_r <- ggplot(as.data.frame(cld.mod[cld.mod$momcross2 %in% contr_order,]),
aes(y=response, x=momcross2, fill=pair)) +
geom_col() +
geom_linerange(aes(ymin=lSE, ymax=uSE)) +
labs(x="Cross", y="Seeds per capsule") +
scale_y_continuous(expand = expand_scale(mult=c(0,0)), breaks = scales::pretty_breaks(n = 15)) +
scale_fill_manual(values=setNames(rep(grn_pur,each=2), c(4, 7, 3, 6)))+ #7,4,6,1
theme_classic() + theme(legend.text=element_text(size=rel(1)), legend.position="bottom", axis.text = element_text(colour="black", size=rel(1)), text=element_text(size=14), axis.ticks.x = element_blank()) + guides(fill=F))
ggsave("hybrid_seedspercapsule_F2.png", em_cont_r, width=7, height = 0.75*7, units="in")