R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

Loading libraries and data

library(readxl)
library(ggplot2)
library(ggpubr)
library(emmeans)
library(multcomp)
library(multcompView)
setwd("C:/Users/backhaua/Desktop/")
dat <- read_excel("qpcr_v2_04_03.xlsx", sheet = "all_data", range = cell_cols("A:O") )

New names:
* `` -> ...1

head(dat)

dat$genoT <- as.factor(paste0(dat$geno1,"_",dat$geno2))
level_order <- c('A', 'C', 'B' , 'S')

Initial plot of all data (2^-DeltaDelta values in all steps)

ggplot(dat, aes(x = factor(POS, level = level_order), y = DD_2,color=WAD)) + geom_point() + facet_grid(dat$genoT, scales = "free") +theme_bw()

I decided to first proceed with NILs as they are important to me and data looks somewhat cleaner

Step1: Looking at data: density plots

dat_NILs <- dat[which(dat$geno2 == "NIL"),]
NILs <- ggplot(dat_NILs, aes(x = factor(POS, level = level_order), y = DD_2,color=genoT)) + geom_point() + facet_grid(rows = vars(WAD), scales = "free")

ggdensity(dat_NILs$DD_2)

qqnorm(dat_NILs$DD_2, pch = 1, frame = FALSE)
qqline(dat_NILs$DD_2, col = "steelblue", lwd = 2)

shapiro.test(dat_NILs$DD_2)

    Shapiro-Wilk normality test

data:  dat_NILs$DD_2
W = 0.45034, p-value = 5.82e-15

Conclusion from density and qq plots (confirmed by shapiro test): Data is skewed. I decide to log2 transform it

plotting density again of Log2 tranformed data

#log2 transform data
dat_NILs$DD_2 <-log2(dat_NILs$DD_2)

ggdensity(dat_NILs$DD_2)

qqnorm(dat_NILs$DD_2, pch = 1, frame = FALSE)
qqline(dat_NILs$DD_2, col = "steelblue", lwd = 2)

shapiro.test(dat_NILs$DD_2)

    Shapiro-Wilk normality test

data:  dat_NILs$DD_2
W = 0.95752, p-value = 0.01603

Looks much better, but still sig dif from normal distribution (shapiro.test). I think its okay?

#step2: ANOVA III to test for interactions

Anova(lm(DD_2 ~  POS * WAD * genoT , data = dat_NILs), type='III', singular.ok = TRUE)

Note: model has aliased coefficients
      sums of squares computed by model comparison

Anova Table (Type III tests)

Response: DD_2
              Sum Sq Df F values    Pr(>F)    
POS           17.704  2  18.1715 1.171e-06 ***
WAD           92.407  2  94.8493 < 2.2e-16 ***
genoT          1.929  1   3.9591    0.0521 .  
POS:WAD       25.907  5  10.6368 5.377e-07 ***
POS:genoT      0.760  2   0.7804    0.4637    
WAD:genoT     24.841  2  25.4975 2.327e-08 ***
POS:WAD:genoT 22.193  5   9.1119 3.289e-06 ***
Residuals     24.356 50                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# I also tried to add REP as block effect (but not sure I wrote it correctly to achieve what I want)
Anova(lm(DD_2 ~  POS * WAD * genoT + REP:WAD, data = dat_NILs), type='III', singular.ok = TRUE)

Note: model has aliased coefficients
      sums of squares computed by model comparison

Anova Table (Type III tests)

Response: DD_2
              Sum Sq Df F values    Pr(>F)    
POS           17.704  2  22.7092 2.300e-07 ***
WAD           48.913  2  62.7436 3.340e-13 ***
genoT          1.641  1   4.2112   0.04658 *  
POS:WAD       25.251  5  12.9561 1.393e-07 ***
POS:genoT      0.760  2   0.9753   0.38565    
WAD:genoT     22.881  2  29.3501 1.227e-08 ***
WAD:REP        8.375  9   2.3873   0.02816 *  
POS:WAD:genoT 21.624  5  11.0954 8.425e-07 ***
Residuals     15.981 41                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Either way, the 3-way interaction is significant so I decide to start with simpler analysis, deciphering the interactions

Looking only at Paragon

are there sig dif between the positions and/or across WAddingtons?

dat_PAR <- dat[which(dat$genoT == "PAR_NIL"),]
PAR <- ggplot(dat_PAR, aes(x = factor(POS, level = level_order), y = DD_2,color=REP)) + geom_point() + facet_grid(rows = vars(WAD), scales = "free") + theme_bw()
PAR

dat_PAR$DD_2 <- log2(dat_PAR$DD_2)

Anova(lm(DD_2 ~  POS * WAD + REP:WAD, data = dat_PAR), type='III', singular.ok = TRUE)

Note: model has aliased coefficients
      sums of squares computed by model comparison

Anova Table (Type III tests)

Response: DD_2
           Sum Sq Df F values    Pr(>F)    
POS       17.7035  2  27.8577 4.315e-06 ***
WAD       11.7840  1  37.0860 1.202e-05 ***
POS:WAD   25.2481  5  15.8919 6.913e-06 ***
WAD:REP    5.5063  7   2.4756   0.06011 .  
Residuals  5.4017 17                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q: At this point, should I plot the transformed data or the raw data? In the end I would probably plot emmeans values anyway? The interaction between position and waddington is again significant, thus I divided the data further by focusing on only one WAD stage at a time. Rep is not sig so I leave it out for post hoc tests

#Only DR

dat_PAR_DR <- dat_PAR[which(dat_PAR$WAD == "DR"),]
Anova(lm(DD_2 ~  POS + REP, data = dat_PAR_DR), type='III', singular.ok = TRUE)

Anova Table (Type III tests)

Response: DD_2
             Sum Sq Df F value   Pr(>F)   
(Intercept) 10.1785  1 66.5978 0.001227 **
POS         17.7035  2 57.9171 0.001114 **
REP          0.9229  2  3.0194 0.158765   
Residuals    0.6113  4                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Position is significant so I ran Tukey/sidak as PH test

flo <- lm(formula = DD_2 ~  POS    , data = dat_PAR_DR,)
marginal = emmeans(flo,  ~ POS, adjust="tukey")
cld(marginal, alpha=0.05, Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

 POS emmean    SE df lower.CL upper.CL .group
 A     2.76 0.292  6     1.80     3.71  a    
 B     5.22 0.292  6     4.26     6.17   b   
 S     6.07 0.292  6     5.11     7.02   b   

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 3 estimates 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 

Conclusion: Expression is sig lower in the apical section, while the basal and stem section are the same (at DR stage) Next, I am testing GP and TS stage

#ONly GP

cld(marginal, alpha=0.05, by = "REP" ,Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

REP = R1:
 POS emmean    SE df lower.CL upper.CL .group
 A   -1.457 0.347  6   -3.005   0.0906  a    
 C   -1.208 0.347  6   -2.756   0.3400  a    
 B   -0.384 0.347  6   -1.932   1.1633  a    
 S    3.004 0.347  6    1.457   4.5522   b   

REP = R3:
 POS emmean    SE df lower.CL upper.CL .group
 A   -2.036 0.347  6   -3.584  -0.4882  a    
 C   -1.787 0.347  6   -3.334  -0.2387  a    
 B   -0.963 0.347  6   -2.511   0.5846  a    
 S    2.426 0.347  6    0.878   3.9735   b   

REP = R4:
 POS emmean    SE df lower.CL upper.CL .group
 A   -0.612 0.347  6   -2.160   0.9360  a    
 C   -0.362 0.347  6   -1.910   1.1854  a    
 B    0.461 0.347  6   -1.087   2.0088  a    
 S    3.850 0.347  6    2.302   5.3976   b   

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 12 estimates 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 

REP is significant but the trend is the same in all reps (only S sig) so I decided to take it out again for post hoc analysis (interaction would not be sig, but I can’t run this)

#post-hoc
flo <- lm(formula = DD_2 ~  POS   , data = dat_PAR_GP,)
marginal = emmeans(flo,  ~ POS, adjust="tukey")
cld(marginal, alpha=0.05 ,Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

 POS emmean    SE df lower.CL upper.CL .group
 A   -1.368 0.481  8    -2.90    0.167  a    
 C   -1.119 0.481  8    -2.65    0.416  a    
 B   -0.296 0.481  8    -1.83    1.240  a    
 S    3.093 0.481  8     1.56    4.628   b   

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 4 estimates 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 

At GP, only the stem is sig different. Basal seems to have a very larger convidence interval so I wonder if I should maybe process the last rep of this condition. When comparing reps, I found that rep4 is sig different from 1 and 3 is not sig dif from either

cld(marginal, alpha=0.05, by = "POS" ,Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

POS = A:
 REP emmean    SE df lower.CL upper.CL .group
 R3  -2.036 0.347  6   -3.584  -0.4882  a    
 R1  -1.457 0.347  6   -3.005   0.0906  ab   
 R4  -0.612 0.347  6   -2.160   0.9360   b   

POS = B:
 REP emmean    SE df lower.CL upper.CL .group
 R3  -0.963 0.347  6   -2.511   0.5846  a    
 R1  -0.384 0.347  6   -1.932   1.1633  ab   
 R4   0.461 0.347  6   -1.087   2.0088   b   

POS = C:
 REP emmean    SE df lower.CL upper.CL .group
 R3  -1.787 0.347  6   -3.334  -0.2387  a    
 R1  -1.208 0.347  6   -2.756   0.3400  ab   
 R4  -0.362 0.347  6   -1.910   1.1854   b   

POS = S:
 REP emmean    SE df lower.CL upper.CL .group
 R3   2.426 0.347  6    0.878   3.9735  a    
 R1   3.004 0.347  6    1.457   4.5522  ab   
 R4   3.850 0.347  6    2.302   5.3976   b   

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 12 estimates 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 

#ONly TS

dat_PAR_TS <- dat_PAR[which(dat_PAR$WAD == "TS"),]
Anova(lm(DD_2 ~  POS + REP, data = dat_PAR_TS), type='III', singular.ok = TRUE)

Anova Table (Type III tests)

Response: DD_2
             Sum Sq Df F value    Pr(>F)    
(Intercept)  47.067  1 98.4168 2.255e-05 ***
POS         125.137  3 87.2205 6.580e-06 ***
REP           0.479  3  0.3342    0.8014    
Residuals     3.348  7                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#rep not sig so left out in PH

#post-hoc
flo <- lm(formula = DD_2 ~  POS   , data = dat_PAR_TS,)
marginal = emmeans(flo,  ~ POS, adjust="tukey")
cld(marginal, alpha=0.05,Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

 POS emmean    SE df lower.CL upper.CL .group
 A   -4.584 0.309 10   -5.520   -3.647  a    
 C   -1.926 0.357 10   -3.007   -0.844   b   
 B    0.794 0.357 10   -0.287    1.876    c  
 S    3.038 0.309 10    2.102    3.974     d 

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 4 estimates 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 

At TS stage, the 4 positions are all sig different. This is what we would have expected from RNA-seq where we saw a step-wise increase from apical to basal section of VRT-2A. New information is that the stem contains high expression. It might be that differences where sig in TS but not GP, because sample collection was much easier for bigger spikes in TS and less plants were pooled per sample. Also I had more reps. In DR, I didnt cut central area which made it easier.

In summary: VRT-2A is sig different expressed between apical and base in DR and TS. It is sig between apcial and stem in all conditions. At TS, central is also sig dif from apical and basal.

#Waddington Next, I want to look at the difference between the three waddington stages. Does expression change over time and is it the same change in all positions?

Anova(lm(DD_2 ~  POS * WAD, data = dat_PAR), type='III', singular.ok = TRUE)

Note: model has aliased coefficients
      sums of squares computed by model comparison

Anova Table (Type III tests)

Response: DD_2
          Sum Sq Df F values    Pr(>F)    
POS       17.704  2   19.476 9.429e-06 ***
WAD       92.407  2  101.657 1.919e-12 ***
POS:WAD   25.907  5   11.400 1.028e-05 ***
Residuals 10.908 24                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

flo <- lm(formula = DD_2 ~  WAD, data = dat_PAR,)
marginal = emmeans(flo,  ~ WAD, adjust="tukey")
cld(marginal, alpha=0.05,Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

 WAD  emmean    SE df lower.CL upper.CL .group
 TS  -0.6840 0.659 32    -2.34    0.975  a    
 GP   0.0776 0.711 32    -1.71    1.870  a    
 DR   4.6798 0.821 32     2.61    6.749   b   

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 3 estimates 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 

So I ran ANOVA by waddington (ignoring pos for now) and found that DR has sig higher expression overall. But the interaction between position and waddington was sig so I should analyse the 2-way ANOVA (And its also interesting)

flo <- lm(formula = DD_2 ~  POS * WAD, data = dat_PAR,)
marginal = emmeans(flo,  ~ POS * WAD, adjust="tukey")
cld(marginal, alpha=0.05,Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

 POS WAD emmean    SE df lower.CL upper.CL .group 
 A   TS  -4.584 0.337 24   -5.648   -3.520  a     
 C   TS  -1.926 0.389 24   -3.154   -0.697   b    
 A   GP  -1.368 0.389 24   -2.597   -0.140   b    
 C   GP  -1.119 0.389 24   -2.348    0.110   bc   
 B   GP  -0.296 0.389 24   -1.524    0.933   bc   
 B   TS   0.794 0.389 24   -0.434    2.023    cd  
 A   DR   2.758 0.389 24    1.529    3.986     de 
 S   TS   3.038 0.337 24    1.974    4.102      e 
 S   GP   3.093 0.389 24    1.865    4.322      e 
 B   DR   5.217 0.389 24    3.988    6.445       f
 S   DR   6.065 0.389 24    4.836    7.294       f
 C   DR  nonEst    NA NA       NA       NA        

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 12 estimates 
P value adjustment: tukey method for comparing a family of 12 estimates 
significance level used: alpha = 0.05 

In Essence: In the apical region, VRT-2A expression decreases significanlty over all time points. In the central region, dif is not sig between GP and TS (DR has none). In the basal region and the stem, expression is sig higher at DR compared to GP and TS.

Comparison of Paragon to POL_NIL

Because I know now that the expression varies between waddington stages I start by analying the expression differences between PAR and POL one waddington stage at a time

#DR stage

dat_PAR_POL_DR <- dat_PAR_POL[which(dat_PAR_POL$WAD == "DR"),]
PARPOL_DR <- ggplot(dat_PAR_POL_DR, aes(x = factor(POS), y = DD_2,color=genoT)) + geom_point()  + theme_bw()
PARPOL_DR

dat_PAR_POL_DR$DD_2 <- log2(dat_PAR_POL_DR$DD_2)
Anova(lm(DD_2 ~  genoT * POS, data = dat_PAR_POL_DR), type='III', singular.ok = TRUE)

Anova Table (Type III tests)

Response: DD_2
             Sum Sq Df F value    Pr(>F)    
(Intercept) 22.8160  1 62.4224 4.271e-06 ***
genoT        1.9286  1  5.2763   0.04041 *  
POS         17.7035  2 24.2176 6.128e-05 ***
genoT:POS    0.7603  2  1.0401   0.38320    
Residuals    4.3861 12                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

flo <- lm(formula = DD_2 ~  genoT * POS, data = dat_PAR_POL_DR,)
marginal = emmeans(flo,  ~ genoT * POS, adjust="tukey")
cld(marginal, alpha=0.05,Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

 genoT   POS emmean    SE df lower.CL upper.CL .group
 PAR_NIL A     2.76 0.349 12     1.66     3.85  a    
 POL_NIL A     3.89 0.349 12     2.80     4.99  ab   
 PAR_NIL B     5.22 0.349 12     4.12     6.31   bc  
 PAR_NIL S     6.07 0.349 12     4.97     7.16    c  
 POL_NIL B     6.27 0.349 12     5.18     7.37    c  
 POL_NIL S     8.03 0.349 12     6.93     9.13     d 

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 6 estimates 
P value adjustment: tukey method for comparing a family of 6 estimates 
significance level used: alpha = 0.05 

Interaction is not sig. POL and PAR are only sig different in stem. But, in POL_NIL the apical, basal and stem are all sig different.

#GP stage

dat_PAR_POL_GP <- dat_PAR_POL[which(dat_PAR_POL$WAD == "GP"),]
PARPOL_GP <- ggplot(dat_PAR_POL_GP, aes(x = factor(POS, levels = level_order), y = DD_2,color=genoT)) + geom_point()  + theme_bw()
PARPOL_GP

dat_PAR_POL_GP$DD_2 <- log2(dat_PAR_POL_GP$DD_2)
Anova(lm(DD_2 ~  genoT * POS, data = dat_PAR_POL_GP), type='III', singular.ok = TRUE)

Anova Table (Type III tests)

Response: DD_2
            Sum Sq Df F value    Pr(>F)    
(Intercept)  5.617  1 10.3688  0.004295 ** 
genoT       34.663  1 63.9864 1.168e-07 ***
POS         38.269  3 23.5475 9.048e-07 ***
genoT:POS    6.994  3  4.3037  0.016977 *  
Residuals   10.834 20                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

flo <- lm(formula = DD_2 ~  genoT * POS, data = dat_PAR_POL_GP,)
marginal = emmeans(flo,  ~ genoT * POS, adjust="tukey")
cld(marginal, alpha=0.05,Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

 genoT   POS emmean    SE df lower.CL upper.CL .group
 PAR_NIL A   -1.368 0.425 20    -2.66  -0.0743  a    
 PAR_NIL C   -1.119 0.425 20    -2.41   0.1751  a    
 PAR_NIL B   -0.296 0.425 20    -1.59   0.9984  a    
 PAR_NIL S    3.093 0.425 20     1.80   4.3873   bc  
 POL_NIL A    3.128 0.368 20     2.01   4.2490   b   
 POL_NIL C    4.029 0.368 20     2.91   5.1498   bc  
 POL_NIL B    4.907 0.368 20     3.79   6.0277    cd 
 POL_NIL S    5.800 0.368 20     4.68   6.9206     d 

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 8 estimates 
P value adjustment: tukey method for comparing a family of 8 estimates 
significance level used: alpha = 0.05 

Interaction, genotype and position are all sig. I decided to run the post-hoc on all data together. Expression is sig higher in POL_NIL in all postions Within POL, apical is sig different from basal and stem. stem and basal are not sig different, neiter are basal vs central.

#TS stage

dat_PAR_POL_TS <- dat_PAR_POL[which(dat_PAR_POL$WAD == "TS"),]
PARPOL_TS <- ggplot(dat_PAR_POL_TS, aes(x = factor(POS, levels = level_order), y = DD_2,color=genoT)) + geom_point()  + theme_bw()
PARPOL_TS

dat_PAR_POL_TS$DD_2 <- log2(dat_PAR_POL_TS$DD_2)
Anova(lm(DD_2 ~  genoT * POS, data = dat_PAR_POL_TS), type='III', singular.ok = TRUE)

Anova Table (Type III tests)

Response: DD_2
             Sum Sq Df F value    Pr(>F)    
(Intercept)  84.043  1 165.591 1.626e-10 ***
genoT        76.632  1 150.990 3.443e-10 ***
POS         127.425  3  83.689 9.168e-11 ***
genoT:POS    28.130  3  18.475 9.947e-06 ***
Residuals     9.136 18                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

flo <- lm(formula = DD_2 ~  genoT * POS, data = dat_PAR_POL_TS,)
marginal = emmeans(flo,  ~ genoT * POS, adjust="tukey")
cld(marginal, alpha=0.05,Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

 genoT   POS emmean    SE df lower.CL upper.CL .group
 PAR_NIL A   -4.584 0.356 18   -5.683   -3.485  a    
 PAR_NIL C   -1.926 0.411 18   -3.194   -0.657   b   
 PAR_NIL B    0.794 0.411 18   -0.474    2.063    c  
 POL_NIL A    2.102 0.411 18    0.833    3.371    cd 
 POL_NIL C    2.559 0.411 18    1.291    3.828    cd 
 POL_NIL B    2.622 0.411 18    1.353    3.891    cd 
 PAR_NIL S    3.038 0.356 18    1.939    4.137     de
 POL_NIL S    4.766 0.411 18    3.497    6.035      e

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 8 estimates 
P value adjustment: tukey method for comparing a family of 8 estimates 
significance level used: alpha = 0.05 

In TS, the difference between POL and PAR is sig dif in apical and central secion but not in basal or stem (although basal PAR is also not sig different from apical or central POL section but it is sig dif to the stem of POL). In POL, only the stem is sig different from all other conditions.

#Only POL

I also wanted to compare just POL_NILs

Anova(lm(DD_2 ~  POS * WAD, data = dat_POL), type='III', singular.ok = TRUE)

Note: model has aliased coefficients
      sums of squares computed by model comparison

Anova Table (Type III tests)

Response: DD_2
           Sum Sq Df F values    Pr(>F)    
POS       25.8847  2  25.0222 8.725e-07 ***
WAD        4.8446  2   4.6832   0.01832 *  
POS:WAD    4.3022  5   1.6635   0.17880    
Residuals 13.4481 26                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

flo <- lm(formula = DD_2 ~  WAD, data = dat_POL,)
marginal = emmeans(flo,  ~ WAD, adjust="tukey")
cld(marginal, alpha=0.05,Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

 WAD emmean    SE df lower.CL upper.CL .group
 TS    3.01 0.408 34     1.99     4.04  a    
 GP    4.47 0.353 34     3.58     5.35   b   
 DR    6.07 0.471 34     4.88     7.25    c  

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 3 estimates 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 

In POl, all waddington stages are sig dif. Expression is decreasing over time, which is consistent with what we knew before

Overall, some positions were also sig dif

flo <- lm(formula = DD_2 ~  POS, data = dat_POL,)
marginal = emmeans(flo,  ~ POS, adjust="tukey")
cld(marginal, alpha=0.05,Letters=letters, adjust="tukey")

Note: adjust = "tukey" was changed to "sidak"
because "tukey" is only appropriate for one set of pairwise comparisons

 POS emmean    SE df lower.CL upper.CL .group
 A     3.05 0.424 33     1.93     4.17  a    
 C     3.40 0.507 33     2.06     4.73  a    
 B     4.63 0.424 33     3.51     5.75  ab   
 S     6.16 0.424 33     5.04     7.28   b   

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 4 estimates 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 

Stem is sig different. But I find it hard to interpret basal. It is not sig dif from apcial and central, but neither from stem.

SO I ran the Dunnetts test to compare everything to apical (because I really want to know if A is dif to C, B or S). Afterall, thats what I found in the RNA-seq

dat_POL$POS <- as.factor(dat_POL$POS)
dat_POL$WAD <- as.factor(dat_POL$WAD)
library(DescTools)
DunnettTest(DD_2 ~ POS, data = dat_POL)

  Dunnett's test for comparing several treatments with a control :  
    95% family-wise confidence level

$A
         diff     lwr.ci   upr.ci    pval    
B-A 1.5821667  0.1020387 3.062295  0.0341 *  
C-A 0.3497857 -1.2812359 1.980807  0.9152    
S-A 3.1093333  1.6292054 4.589461 4.6e-05 ***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Then A to B comarision is sig, I guess because there are less comparisons? Maybe we could use Dunnett to always compare S,B and C to A?