The statistical approaches for assessing flavonoid biosynthesis pathway genes under two-spotted spider mite treatment

library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggfortify)
library(arm)

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## Loading required package: Matrix

## Loading required package: lme4

## 
## arm (Version 1.15-3, built: 2026-4-15)

## Working directory is D:/Ph.d/course Fall 1/BiostatisticR/Notebook/Session3

library(cowplot)

library(readr)
Test4 <- read_csv("D:/Ph.d/course Fall 1/BiostatisticR/Notebook/Session3/Gene.csv")

## Rows: 40 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Plant.ID, Treatment
## dbl (2): Replicate, CHS.express
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(Test4)

## Rows: 40
## Columns: 4
## $ Plant.ID    <chr> "P1", "P2", "P3", "P4", "P5", "P6", "P7", "P8", "P9", "P10…
## $ Treatment   <chr> "Control", "Control", "Control", "Control", "Control", "Co…
## $ Replicate   <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, …
## $ CHS.express <dbl> 1.898, 2.078, 2.004, 1.661, 1.217, 1.680, 2.053, 1.204, 1.…

ggplot(Test4, aes(x = Treatment, y = CHS.express)) +
  geom_boxplot() +
  theme_bw()

ggplot(Test4, aes(x = Treatment, y = CHS.express)) +
  geom_boxplot() +
  theme_bw() +
  coord_flip()

model_gene <- lm(CHS.express ~ Treatment, data = Test4)
autoplot(model_gene, smooth.colour = NA)

## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Removed 40 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 40 rows containing missing values or values outside the scale range
## (`geom_line()`).

anova(model_gene)

## Analysis of Variance Table
## 
## Response: CHS.express
##           Df Sum Sq  Mean Sq F value  Pr(>F)  
## Treatment  3 0.8825 0.294155  2.9862 0.04384 *
## Residuals 36 3.5462 0.098506                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model_gene)

## 
## Call:
## lm(formula = CHS.express ~ Treatment, data = Test4)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.57820 -0.20155  0.06475  0.22945  0.60140 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           1.72940    0.09925  17.425   <2e-16 ***
## TreatmentMite_Ten    -0.05620    0.14036  -0.400   0.6912    
## TreatmentMite_Thirty  0.23520    0.14036   1.676   0.1025    
## TreatmentMite_Twenty  0.29190    0.14036   2.080   0.0447 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3139 on 36 degrees of freedom
## Multiple R-squared:  0.1993, Adjusted R-squared:  0.1325 
## F-statistic: 2.986 on 3 and 36 DF,  p-value: 0.04384

mean_gene <- Test4 %>%
  group_by(Treatment) %>%
  summarise(mean.CHS = mean(CHS.express))
control_mean <- mean(Test4$CHS.express[Test4$Treatment == "Control"])

sumDat <- Test4 %>%
  group_by(Treatment) %>%
  summarise(meanCHS = mean(CHS.express))
control_mean <- sumDat$meanCHS[sumDat$Treatment == "Control"]

ggplot(Test4, aes(x = CHS.express, y = Treatment,
colour = Treatment)) +
geom_vline(xintercept = control_mean,
           linetype = "dashed", colour = "grey40") +
geom_segment(data = sumDat,
             aes(x = meanCHS, xend = control_mean,
                 yend = Treatment),
             colour = "black") +
geom_point(alpha = 0.4, size = 2) +
geom_point(data = sumDat, aes(x = meanCHS),
           shape = 18, size = 5) +
scale_y_discrete(limits = c("Control", "Mite_Ten",
                           "Mite_Twenty", "Mite_Thirty")) +
scale_colour_manual(values = c(
  Control = "red",
  Mite_Ten = "green4",
  Mite_Twenty = "turquoise3",
  Mite_Thirty = "purple"
)) +
labs(x = "CHS expression", y = "Treatment",
     colour = "Treatment") +
theme_bw()

sumDat <- Test4 %>%
  group_by(Treatment) %>%
  summarise(meanCHS = mean(CHS.express))

ggplot(Test4, aes(x = CHS.express, y = Treatment,
colour = Treatment)) +
geom_point(alpha = 0.4, size = 2) +
geom_point(data = sumDat, aes(x = meanCHS),
           shape = 18, size = 5) +
scale_y_discrete(limits = c("Control", "Mite_Ten",
                           "Mite_Twenty", "Mite_Thirty")) +
scale_colour_manual(values = c(
  Control = "red",
  Mite_Ten = "green4",
  Mite_Twenty = "turquoise3",
  Mite_Thirty = "purple"
)) +
labs(x = "CHS expression", y = "Treatment",
     colour = "Treatment") +
theme_bw()