#Bc. Jaroslav Rohel, příloha diplomové práce: Endopolyploidie v přírodních populacích diploidních rostlin a jejich polyploidních příbuzných. Masarykova univerzita, Přírodovědecká fakulta. Brno 2026. Dostupné z: https://is.muni.cz/auth/th/vboor/
#Elektronická příloha 5: Analýzy dat 


library(lme4)
library(emmeans)
library(lmerTest)
library(dplyr)

dat <- read.csv("dat.csv", header=T, sep=",")

#Faktory ovlivňující míru endopolyploidie--------
mm<-lm(sqrt(ei)~Xmonoploid+Genus*part*ploidy*position, data=dat)
anova(mm)
mm2 <- update(mm, ~.-Genus:part:ploidy:position)
anova(mm2) #odstraním nevýznamnou čtyřnásobnou interakci
mm3 <- update(mm2, ~.-part:ploidy:position)
anova(mm3) #odstraním nevýznamnou trojnásobnou interakci
#nic dalšího už úplně odstranit nemůžu - jsou nějaké nevýznamné dvojnásobné, ale ty jsou zahrnuty i v trojnásobných - není to s tím tak jendoduché
summary(mm3)

plot(mm3)






mmm<-lm(sqrt(ei)~Genus*part*ploidy*position, data=dat)
anova(mmm)
summary(mmm)

plot(mm)


#Rozdíl endopolyploidie v čepeli a řapíku----------------

mm <- lmer(sqrt(ei) ~ part * position + (1 + part | taxon) + (1 | taxon:code) + (1 | taxon:code:position), data = dat) 
anova(mm)
summary(mm)


#heteroskedsticita rozptylu
plot(fitted(mm), resid(mm))
abline(h = 0, col="red")

#normalita reziduí
qqnorm(resid(mm))
qqline(resid(mm), col="red")



emm <- emmeans(mm, ~ part | position) #using Satterthwaite's method stejně jako anova
pairs(emm, adjust = "tukey")


#Vliv pozice (stáří) listu---------------

#pro řapík
dat_rapik_hel_contr <- subset(dat, part=="rapik")
dat_rapik_hel_contr <- droplevels(dat_rapik_hel_contr)
dat_rapik_hel_contr$position<- factor(dat_rapik_hel_contr$position, levels = c("dolni", "prostredni", "horni"), ordered = TRUE)
#dat_rapik_hel_contr$position<- factor(dat_rapik_hel_contr$position, levels = c("horni", "prostredni", "dolni"), ordered = TRUE)
contrasts(dat_rapik_hel_contr$position) <- contr.helmert(3)
contrasts(dat_rapik_hel_contr$position)
ml <- lmer(sqrt(ei) ~ position + (1 | taxon/code), data=dat_rapik_hel_contr) 
anova(ml)
summary(ml)


#heteroskedsticita rozptylu
plot(fitted(ml), resid(ml))
abline(h = 0, col="red")

#normalita reziduí
qqnorm(resid(ml))
qqline(resid(ml), col="red")


#pro čepel
dat_cepel <- read.csv("dat_cepel.csv", header=T, sep=",")

dat_cepel_hel_contr <- dat_cepel
dat_cepel_hel_contr <- droplevels(dat_cepel_hel_contr)
dat_cepel_hel_contr$position<- factor(dat_cepel_hel_contr$position, levels = c("dolni", "prostredni", "horni"), ordered = TRUE)
#dat_cepel_hel_contr$position<- factor(dat_cepel_hel_contr$position, levels = c("horni", "prostredni", "dolni"), ordered = TRUE)
contrasts(dat_cepel_hel_contr$position) <- contr.helmert(3)
contrasts(dat_cepel_hel_contr$position)
ml <- lmer(sqrt(ei) ~ position + (1 | taxon/code), data=dat_cepel_hel_contr) 
anova(ml)
summary(ml)


#heteroskedsticita rozptylu
plot(fitted(ml), resid(ml))
abline(h = 0, col="red")

#normalita reziduí
qqnorm(resid(ml))
qqline(resid(ml), col="red")

#Vliv velikosti genomu---------

mg2 <- lmer(sqrt(ei) ~ X2C_log * part * position +  (1 | code), data = dat)
anova(mg2)
mg3 <- update(mg2, . ~ . - X2C_log:part:position)
anova(mg3)
mg4 <- update(mg3, . ~ . - X2C_log:position)
anova(mg4)
mg5 <- update(mg4, . ~ . - X2C_log:part)
anova(mg5)
mg6 <- update(mg5, . ~ . - part:position)
anova(mg6)
summary(mg6)


#heteroskedsticita rozptylu
plot(fitted(mg6), resid(mg6))
abline(h = 0, col="red")

#normalita reziduí
qqnorm(resid(mg6))
qqline(resid(mg6), col="red")




#Vliv přírodních podmínek----------
dat_veg <- read.csv("dat_veg.csv", header=T, sep=",")
#byl použit následující model pro jednotlivé kombinace všech části a pozice listů  
#a jednotlivě pro všechny vypočtené indikační hodnoty: "L.CZ_prum", "T.CZ_prum", "M.CZ_prum", "N.CZ_prum", "R.CZ_prum", "S.CZ_prum", "LIGHT_EU_prum", "TEMPERATURE_EU_prum", "MOISTURE_EU_prum", "NUTRIENTS_EU_prum", "REACTION_EU_prum"

#příklad modelu pro Nutrients CZ, prostredni rapik
dat_veg_pr <- subset(dat_veg, position=="prostredni" & part=="rapik" )

m<-lmer(sqrt(ei)~N.CZ_prum + (N.CZ_prum|taxon), data=dat_veg_pr)
anova(m)
summary(m) 


#heteroskedsticita rozptylu
plot(fitted(m), resid(m))
abline(h = 0, col="red")

#normalita reziduí
qqnorm(resid(m))
qqline(resid(m), col="red")


#vegetace - korekce na mnohonásobná porovnání Li & Ji
prediktory_CZ <- data.frame(var=c("L.CZ_prum", "T.CZ_prum", "M.CZ_prum", "N.CZ_prum", "R.CZ_prum", "S.CZ_prum", "LIGHT_EU_prum"))

prediktory_EU <- data.frame(var=c("LIGHT_EU_prum", "TEMPERATURE_EU_prum", "MOISTURE_EU_prum", "NUTRIENTS_EU_prum", "REACTION_EU_prum"))

prediktory_all <- data.frame(var=c("L.CZ_prum", "T.CZ_prum", "M.CZ_prum", "N.CZ_prum", "R.CZ_prum", "S.CZ_prum", "LIGHT_EU_prum", "TEMPERATURE_EU_prum", "MOISTURE_EU_prum", "NUTRIENTS_EU_prum", "REACTION_EU_prum"))


#výpočet correlation matrix
correlation_matrix <- cor(dat_veg[,prediktory_CZ$var], use = "complete.obs")  #tady je potřeba napsat které prediktory chci

# get eigenvalues
eigenvalues <- eigen(correlation_matrix)$values

# Apply Li and Ji's method to calculate the effective number of tests (M_eff)
effective_tests <- sum(eigenvalues > 1) + sum(eigenvalues[eigenvalues <= 1])

# Step 4: Print the results
cat("Eigenvalues:\n", eigenvalues, "\n")
cat("Effective number of tests (M_eff):", effective_tests, "\n")



#Vliv velikosti listu---------
df <- read.csv("df_velikost_listu.csv", header=T, sep=",")

m2 <- lmer(sqrt(ei) ~ velikost_listu_cm_log*part + (1| taxon/position), data=df)
anova(m2)
summary(m2)


#heteroskedsticita rozptylu
plot(fitted(m2), resid(m2))
abline(h = 0, col="red")

#normalita reziduí
qqnorm(resid(m2))
qqline(resid(m2), col="red")

emtrends(m2, specs = ~ part,  var = "velikost_listu_cm_log")
tr <- emtrends(m2, ~ part, var = "velikost_listu_cm_log")
summary(tr, infer = c(TRUE, TRUE))




#Vliv životní formy----------
endoploidy_zivotFormy <- read.csv("endoploidy_zivotFormy.csv", header=T, sep=",")

#pro všechna pletiva a všechny pozice dohromady
endoploidy_zivotFormy_f <-endoploidy_zivotFormy
m <- lmer(sqrt(ei) ~ Životní_forma_._2023_08_18 -1 + (1 | taxon), data = endoploidy_zivotFormy_f)
anova(m)
summary(m)

unique(endoploidy_zivotFormy_f$Životní_forma_._2023_08_18)

#pro prostřední řapík
endoploidy_zivotFormy_f <- subset(endoploidy_zivotFormy, endoploidy_zivotFormy$part=="rapik" & endoploidy_zivotFormy$position=="prostredni")
m <- lmer(sqrt(ei) ~ Životní_forma_._2023_08_18 - 1 + (1 | taxon), data = endoploidy_zivotFormy_f)
anova(m)
summary(m)

#pro prostřední čepel
endoploidy_zivotFormy_f <- subset(endoploidy_zivotFormy, endoploidy_zivotFormy$part=="cepel" & endoploidy_zivotFormy$position=="prostredni")
m <- lmer(sqrt(ei) ~ Životní_forma_._2023_08_18 - 1 + (1 | taxon), data = endoploidy_zivotFormy_f)
anova(m)
summary(m)


#heteroskedsticita rozptylu
plot(fitted(m), resid(m))
abline(h = 0, col="red")

#normalita reziduí
qqnorm(resid(m))
qqline(resid(m), col="red")



#postupné zjednodušování modelu,
#pro všechna pletiva dohromady
endoploidy_zivotFormy_f <-endoploidy_zivotFormy

#pro prostřední řapík
#endoploidy_zivotFormy_f <- subset(endoploidy_zivotFormy, endoploidy_zivotFormy$part=="rapik" & endoploidy_zivotFormy$position=="prostredni")
#pro prostřední čepel
#endoploidy_zivotFormy_f <- subset(endoploidy_zivotFormy, endoploidy_zivotFormy$part=="rapik" & endoploidy_zivotFormy$position=="prostredni")

endoploidy_zivotFormy_f$Životní_forma_._2023_08_18 <- as.factor(endoploidy_zivotFormy_f$Životní_forma_._2023_08_18)

m <- lmer(sqrt(ei) ~ Životní_forma_._2023_08_18 - 1 + (1 | taxon), data = endoploidy_zivotFormy_f, REML=FALSE)
anova(m)
summary(m)
endoploidy_zivotFormy_f$Životní_forma2 <- endoploidy_zivotFormy_f$Životní_forma_._2023_08_18
levels(endoploidy_zivotFormy_f$Životní_forma2)
levels(endoploidy_zivotFormy_f$Životní_forma2)[4:5] <- "makronanofanerofyt"

m2<- lmer(sqrt(ei) ~ Životní_forma2 -1 + (1 | taxon), data = endoploidy_zivotFormy_f, REML=FALSE)
anova(m, m2) #modely se neliší, mohu je zjednodušit
summary(m2)
endoploidy_zivotFormy_f$Životní_forma3 <- endoploidy_zivotFormy_f$Životní_forma2
levels(endoploidy_zivotFormy_f$Životní_forma3)
levels(endoploidy_zivotFormy_f$Životní_forma3)[c(3, 4)] <- "makronanofanerofyt_chamaefyt"

m3<- lmer(sqrt(ei) ~ Životní_forma3 -1 + (1 | taxon), data = endoploidy_zivotFormy_f, REML=FALSE)
anova(m2, m3) #modely se neliší, mohu je zjednodušit
summary(m3)

endoploidy_zivotFormy_f$Životní_forma4 <- endoploidy_zivotFormy_f$Životní_forma3
levels(endoploidy_zivotFormy_f$Životní_forma4)
levels(endoploidy_zivotFormy_f$Životní_forma4)[c(1, 2)] <- "geofyt_hemikryptofyt"

m4<-lmer(sqrt(ei) ~ Životní_forma4 -1 + (1 | taxon), data = endoploidy_zivotFormy_f, REML=FALSE)
anova(m3, m4) #modely se neliší, mohu je zjednodušit
summary(m4)

endoploidy_zivotFormy_f$Životní_forma5 <- endoploidy_zivotFormy_f$Životní_forma4
levels(endoploidy_zivotFormy_f$Životní_forma5)
levels(endoploidy_zivotFormy_f$Životní_forma5)[c(1, 3)] <- "geo_tero_hemikryptofyt"

m5<-lm(sqrt(ei)~Životní_forma5 -1, data=endoploidy_zivotFormy_f, REML=FALSE)
anova(m4, m5) #modely už se liší, nemohu zjedodušit!
summary(m5)


#jiná varianta sloučení
endoploidy_zivotFormy_f$Životní_forma6 <- endoploidy_zivotFormy_f$Životní_forma3
levels(endoploidy_zivotFormy_f$Životní_forma6)
levels(endoploidy_zivotFormy_f$Životní_forma6)[c(1, 4)] <- "geo_tero_fyt"

m6<-lm(sqrt(ei)~Životní_forma6 -1, data=endoploidy_zivotFormy_f, REML=FALSE)
anova(m4, m6) #modely se (ne)liší, nemohu zjedodušit
summary(m6)



endoploidy_zivotFormy_f$Životní_forma7 <- endoploidy_zivotFormy_f$Životní_forma6
levels(endoploidy_zivotFormy_f$Životní_forma7)
levels(endoploidy_zivotFormy_f$Životní_forma7)[c(1, 3)] <- "geo_tero_hemikryptofyt"
m7<-lm(sqrt(ei)~Životní_forma7 -1, data=endoploidy_zivotFormy_f, REML=FALSE)
anova(m6, m7) #modely už se liší, nemohu zjedodušit!

#vyřadit dřeviny
endoploidy_zivotFormy_f <- endoploidy_zivotFormy
dreviny <- c("Ginkgo","Betula","Salix")
endoploidy_zivotFormy_f_nedreviny <- subset(endoploidy_zivotFormy_f, !rod %in% dreviny)
endoploidy_zivotFormy_f_nedreviny$Životní_forma_._2023_08_18 <- droplevels(endoploidy_zivotFormy_f_nedreviny$Životní_forma_._2023_08_18)

m <- lmer(sqrt(ei) ~ Životní_forma_._2023_08_18 + (1 | taxon), data = endoploidy_zivotFormy_f_nedreviny)
anova(m)
summary(m)

#velikost genomu
endoploidy_zivotFormy_f <- endoploidy_zivotFormy
endoploidy_zivotFormy_f$X2C_log <- log(endoploidy_zivotFormy_f$X2C)
m <- lmer(sqrt(ei) ~  X2C_log + Životní_forma_._2023_08_18 + (1 | taxon), data = endoploidy_zivotFormy_f)
anova(m)
summary(m)


#Tukey - testování rozdílnosti mezi jednotlivých životními formami + pro vykreslení k boxplotům
endoploidy_zivotFormy_f <- subset(endoploidy_zivotFormy, endoploidy_zivotFormy$part=="rapik" & endoploidy_zivotFormy$position=="prostredni")
model_lmm <- lmer(sqrt(ei) ~ Životní_forma_._2023_08_18  -1 + (1 | taxon), data = endoploidy_zivotFormy_f)

# 2) ANOVA tabulka fixních efektů
anova(model_lmm)

# 3) Post-hoc test (Tukey)
emm <- emmeans(model_lmm, "Životní_forma_._2023_08_18")
pairs_res <- pairs(emm, adjust = "tukey")



#Variabilita v EI mezi částmi, pozicemi a velikostmi listu---------


#mezi částmi a pozicemi listů - postupně pro všechny kombinace řapíku a čepele, data už jsou vyfiltorvaná že tan nejsou dvojice bez řapíku
z <- lm(sqrt(ei)~taxon, data=dat[which(dat$age_pletivo=="horni cepel"), ])
anova(z)
summary(z)

#heteroskedsticita rozptylu
plot(fitted(z), resid(z))
abline(h = 0, col="red")

#normalita reziduí
qqnorm(resid(z))
qqline(resid(z), col="red")

#mezí částmi a velikostmi listu
#pro řapíky
data2v <- read.csv("data2v.csv", header=T, sep=",")
data2v_n <- subset(data2v, rod!="Gagea")
z <- lm(sqrt(ei)~taxon, data=data2v_n[which(data2v_n$velikost_pletivo=="velky rapik"), ]) #a postupně nahrazovat za různé veliksoti řapíku
anova(z)
summary(z)

#pro čepele stejných rosltin jako pro které mám řapíky
taxony <- unique(data2v_n$taxon[data2v_n$velikost_pletivo == "velky rapik" |    data2v_n$velikost_pletivo == "stredni rapik" |    data2v_n$velikost_pletivo == "maly rapik"])
data2v_n_filt <- data2v_n[data2v_n$taxon %in% taxony, ]
levels(data2v_n_filt$taxon)
data2v_n_filt <- droplevels(data2v_n_filt)
levels(data2v_n_filt$taxon)

z <- lm(sqrt(ei)~taxon, data=data2v_n_filt[which(data2v_n_filt$velikost_pletivo=="velky cepel"), ]) #a postupně nahrazovat za různé veliksoti řapíku a čepele
anova(z)
summary(z)




#variabilita v poziciích, ale pro ty druhy pro které mám velikost
#řapíky
data2v_n$age_pletivo <- paste(data2v_n$position, data2v_n$part)
z <- lm(sqrt(ei)~taxon, data=data2v_n[which(data2v_n$age_pletivo=="horni rapik"), ]) #a postupně nahrazovat za různé veliksoti řapíku a čepele
anova(z)
summary(z)

#čepele
data2v_n_filt$age_pletivo <- paste(data2v_n_filt$position, data2v_n_filt$part)
z <- lm(sqrt(ei)~taxon, data=data2v_n_filt[which(data2v_n_filt$age_pletivo=="horni cepel"), ]) #a postupně nahrazovat za různé veliksoti řapíku a čepele
anova(z)
summary(z)


#Korelace mezi částmi listu a pozicemi


dat_wide <- read.csv("dat_wide.csv", header=T, sep=",")


qqnorm(dat_wide$rapik_dolni)
qqline(dat_wide$rapik_dolni)
shapiro.test(dat_wide$rapik_dolni) # není normální rozdělení, p<0.05,  -aby to bylo normální, nechci zamítnout nulovou hypotézu

#data nejsou normální, ale jsou blízká normálním. Zároveň mám dost měření - přes 30, můžu tak použít centrální limitní větu a porovnat to parametrickým spearmenem (zákaldní), a nedělat to neparametrickým Pearsonem

#korelace, use="pairwise.complete.obs" zajišťuje, že se při výpočtu korelace ignorují chybějící hodnoty.
cor(dat_wide %>% select(starts_with("rapik")), use="pairwise.complete.obs")
cor(dat_wide %>% select(starts_with("cepel")), use="pairwise.complete.obs")

cor(dat_wide$rapik_dolni, dat_wide$cepel_dolni, use="pairwise.complete.obs")
cor(dat_wide$rapik_prostredni, dat_wide$cepel_prostredni, use="pairwise.complete.obs")
cor(dat_wide$rapik_horni, dat_wide$cepel_horni, use="pairwise.complete.obs")

cor(dat_wide %>% select(starts_with("rapik") | starts_with("cepel")), use="pairwise.complete.obs")



#Liší se jádra endopolyploidních diploidů a polyploidů ve finální velikosti genomu?------

dat_gen_wide <- read.csv("dat_gen_wide.csv", header=T, sep=",")

#chci otestovat, že se diploid a polyploid liší v reálném genomu:
# Vypočítat rozdíl -ODMOCNINOVĚ TRANSFORMUJU VELIKOSTI GENOMU, ze kterých počítám rozdíl
dat_gen_wide <- dat_gen_wide %>%  mutate(diff = log(P) - log(D))
# Histogram
hist(dat_gen_wide$diff, main = "Rozdíl genomu (P-D)", xlab = "Rozdíl [Mbp]")
# Q-Q plot - měly by být nromální rozdíly mezi skupinami
qqnorm(dat_gen_wide$diff)
qqline(dat_gen_wide$diff, col = "red")
# Shapiro-Wilk test
shapiro.test(dat_gen_wide$diff) #zamítáme H0, že data mají nromální rozložení
#oetstuju to - nenormální => wilcox test
wilcox.test(log10(dat_gen_wide$D), log10(dat_gen_wide$P), paired = TRUE) #p<0.0001 - zamítáme nulovou hypotézu, že diploid a polyploid se neliší v realizované velikosti genomu
#nemohu použíét t.test - data nejsou normální
t.test(dat_gen_wide$D, dat_gen_wide$P, paired = TRUE)

median(dat_gen_wide$diff)
median(dat_gen_wide$D)
median(dat_gen_wide$P)

#v datehc byla jedna ulítlá hodnota - viz qqline, odstraním ji
dat_gen_wide$diff #je to řádek 19
dat_gen_wide[19,] # pro rumex, tam byl ten 20 ploidní polyploid - ten nemá úplně smysl porovnávat, odstraním
dat_gen_wide2 <- dat_gen_wide[-19,]
dat_gen_wide2$diff

# Histogram - teď už to vypadá lépe
hist(dat_gen_wide2$diff, main = "Rozdíl genomu (P-D)", xlab = "Rozdíl [Mbp]")
# Q-Q plot - měly by být nromální rozdíly mezi skupinami
qqnorm(dat_gen_wide2$diff)
qqline(dat_gen_wide2$diff, col = "red") #to už vypadá normálně
# Shapiro-Wilk test
shapiro.test(dat_gen_wide2$diff) #neliší se od normálního rozdělení, mohu použít t test
t.test(dat_gen_wide2$D, dat_gen_wide2$P, paired = TRUE)
t.test(log10(dat_gen_wide2$D), log10(dat_gen_wide2$P), paired = TRUE)


#na netransformovaných dfatech - wilcox test protže nebylo normální
wilcox.test(dat_gen_wide$D, dat_gen_wide$P, paired = TRUE) 


mean(dat_gen_wide2$D)
mean(dat_gen_wide2$P)
dat_gen_wide2 <- (dat_gen_wide2 %>%  mutate(rozdil_skut = P - D))
dat_gen_wide2 <- (dat_gen_wide2 %>%  mutate(rozdil_skut_log = log10(P) - log10(D)))
mean(dat_gen_wide2$rozdil_skut)
print(dat_gen_wide2, n=30) #jediná dvojice, kde diploid má většií skutečnou velikost genomu než polyploid je draba


#porovnání rozdílů mezi rozdíly X2C a skutečným obsahem DNA v řapíku

dat_gen_wide2 <- dat_gen_wide2 %>%  mutate(rozdil_2C = X2C_P - X2C_D)
dat_gen_wide2 <- dat_gen_wide2 %>%  mutate(rozdil_2C_od_skut = rozdil_2C - rozdil_skut)


#bez log transforamce dat
# Histogram - teď už to vypadá lépe
hist(dat_gen_wide2$rozdil_2C_od_skut, main = "Rozdíl genomu (P-D)", xlab = "Rozdíl [Mbp]")
# Q-Q plot - měly by být nromální rozdíly mezi skupinami
qqnorm(dat_gen_wide2$rozdil_2C_od_skut)
qqline(dat_gen_wide2$rozdil_2C_od_skut, col = "red") #to už vypadá normálně
# Shapiro-Wilk test
shapiro.test(dat_gen_wide2$rozdil_2C_od_skut) #neliší se od normálního rozdělení, mohu použít t test
t.test(dat_gen_wide2$rozdil_2C, dat_gen_wide2$rozdil_skut, paired = TRUE) #neprokázal jsem, že by se rozdíly v X2C a v průměrném obsahu DNA lišily od nuly



