library(ggplot2)
library(corrplot)
library(dplyr)
library(tidyr)
library(kableExtra)
library(corrplot)

DA_list <- c("Comt", "Ddc", "Drd1", "Drd2", "Drd3", "Drd4", "Drd5", "Slc18a2", "Comt", "Maoa", "Maob")

PCR results

With the applied primers, we received results for Comt, Drd1, Drd2, Maob, Drd5, Th, Ddc, and Maoa, raw 2−ΔΔCt data are represented below in the Table S2-1. The Shapiro–Wilk test results are summarized in the Table S2-2.

results_pcr <- read.csv("results_pcr.csv", sep="")
knitr::kable(results_pcr, 
             caption = "Table S4-1. RT-PCR results",
             digits = 2,
             align = 'c')
Table S4-1. RT-PCR results
Comt Drd1 Drd2 Maob Drd5 Th Ddc Maoa Group
1.53 0.99 0.98 1.30 1.05 0.65 1.13 1.01 WT
0.97 1.34 1.01 1.01 0.78 0.62 0.98 0.80 WT
0.86 0.78 1.10 0.68 0.80 1.31 1.15 0.81 WT
1.02 1.45 1.22 1.09 1.11 2.72 0.92 0.80 WT
0.82 0.91 0.76 1.17 0.51 0.99 1.17 1.36 WT
0.69 1.74 1.57 0.80 2.69 1.21 0.74 1.12 WT
1.36 0.43 0.63 1.10 2.69 0.57 0.99 1.26 WT
0.76 2.60 1.62 0.87 2.23 1.64 0.57 1.65 Het
1.44 1.98 1.56 0.80 0.99 1.00 0.80 1.26 Het
1.21 2.10 1.36 0.66 0.98 0.66 1.42 1.06 Het
2.18 0.46 0.98 0.92 0.41 1.18 1.22 1.05 Het
3.42 1.14 1.06 1.17 0.74 1.47 1.34 1.03 Het
1.29 0.83 1.19 0.94 0.64 1.11 1.04 1.18 Het
0.53 5.44 2.55 0.72 1.33 0.98 0.96 0.95 Het
1.20 1.09 1.15 0.96 0.96 0.58 1.05 1.15 Het
1.40 6.51 1.75 0.97 0.94 1.26 0.40 1.00 KO
1.10 0.82 0.99 1.12 0.75 0.91 1.20 1.65 KO
0.59 1.17 1.08 0.95 1.48 1.74 0.48 1.24 KO
0.79 1.55 1.31 1.00 1.35 0.87 0.77 1.44 KO
0.48 1.67 1.12 1.03 1.03 0.89 0.66 0.94 KO
0.44 2.42 1.70 1.12 1.08 0.84 0.82 1.08 KO
2.20 0.61 0.53 0.97 0.74 0.60 1.80 1.69 KO 
shapiro_tab <- results_pcr %>%
select(where(is.numeric), "Group") %>% 
  pivot_longer(cols = -one_of("Group"), names_to = "Gene", values_to = "Val") %>%
  group_by(`Group`, Gene) %>%
  summarise(
    n_samples = n(),
    p_val = if(n() >= 3) shapiro.test(Val)$p.value else NA_real_,
    .groups = "drop"
  ) %>%
    pivot_wider(names_from = Gene, 
              values_from = p_val)

table_output <- shapiro_tab %>%
  mutate(across(where(is.numeric), ~ round(., 4))) %>%
  kable(booktabs = TRUE, align = "c", escape = FALSE,
        caption = "Table S4-2. P-values Shapiro-Wilk (P values below 0.05 for non-normal distribution are in bold)") %>%
  kable_styling(full_width = F)

for (i in 2:ncol(shapiro_tab)) {
  bold_logic <- shapiro_tab[[i]] < 0.05
  table_output <- table_output %>%
    column_spec(i, bold = bold_logic)
}
table_output
Table S4-2. P-values Shapiro-Wilk (P values below 0.05 for non-normal distribution are in bold)
Group n_samples Comt Ddc Drd1 Drd2 Drd5 Maoa Maob Th
Het 8 0.1227 0.8935 0.0514 0.0428 0.1292 0.0492 0.8590 0.8038
KO 7 0.1867 0.2481 0.0066 0.6759 0.5200 0.3834 0.0965 0.1158
WT 7 0.4009 0.3394 0.9667 0.8790 0.0226 0.2058 0.6692 0.0264 
run_anova_summary <- function(data) {
    group_col <- names(data)[ncol(data)]
    genes <- names(data)[-ncol(data)]
    
    for (gene in genes) {
        formula <- as.formula(paste(gene, "~", group_col))
        res_anova <- aov(formula, data = data)
        p_val <- summary(res_anova)[[1]][["Pr(>F)"]][1]
        
        cat(paste0("Gene: ", gene, "\n"))
        
        if (p_val > 0.05) {
            cat("Result: No significant differences (p =", round(p_val, 4), ")\n")
        } else {
            cat("Result: SIGNIFICANT (p =", round(p_val, 4), ")\n")
            # Print the ANOVA table for detail
            print(summary(res_anova))
        }
        cat("-----------------------------------\n")
    }
}

run_anova_summary(results_pcr[, c(1, 4, 7, 9)])
## Gene: Comt
## Result: No significant differences (p = 0.2982 )
## -----------------------------------
## Gene: Maob
## Result: No significant differences (p = 0.1586 )
## -----------------------------------
## Gene: Ddc
## Result: No significant differences (p = 0.583 )
## -----------------------------------
run_KW_summary <- function(data) {
    group_col <- names(data)[ncol(data)]
    genes <- names(data)[-ncol(data)]
    
    for (gene in genes) {
        formula <- as.formula(paste(gene, "~", group_col))
        res_KW <- kruskal.test(formula, data = data)
        p_val <- res_KW$p.value
    
    cat(paste0("Gene: ", gene, "\n"))
    
    if (p_val > 0.05) {
      cat("Result: No significant differences (p =", round(p_val, 4), ")\n")
    } else {
      cat("Result: SIGNIFICANT (Kruskal-Wallis, p =", round(p_val, 4), ")\n")
      print(res_KW)
  
      cat("Recommended: Run Dunn's test for this gene.\n")
    }
    cat("-----------------------------------\n")
  }
}

run_KW_summary(results_pcr[, c(2, 3, 5, 6, 8, 9)])
## Gene: Drd1
## Result: No significant differences (p = 0.3883 )
## -----------------------------------
## Gene: Drd2
## Result: No significant differences (p = 0.232 )
## -----------------------------------
## Gene: Drd5
## Result: No significant differences (p = 0.6542 )
## -----------------------------------
## Gene: Th
## Result: No significant differences (p = 0.855 )
## -----------------------------------
## Gene: Maoa
## Result: No significant differences (p = 0.2951 )
## -----------------------------------
Fig1_tab <- results_pcr %>%
select(where(is.numeric), "Group") %>% 
  pivot_longer(cols = -one_of("Group"), names_to = "Gene", values_to = "Val")
ggplot(Fig1_tab, aes(x = Gene, y = Val, fill = factor(Group, levels = c("WT", "Het", "KO")))) +
 geom_boxplot(alpha=0.6, outlier.shape = NA, width = 0.9)+
  geom_jitter(pch = 21, size=2, position = position_dodge(width = 0.8))+
scale_fill_manual(values = c("forestgreen", "purple", "magenta2"))+
theme(axis.text=element_text(size=16), axis.title.y = element_text(size=18), axis.title.x = element_text(size=18)
)+
ylab(expression("Normalized expression"))+
theme(legend.title = element_blank())+
theme(legend.text = element_text(size=16))+
theme(axis.text.x = element_text(face = "italic"))

Given that the majority of our data values exhibited a non-normal distribution, we determined that parametric methods were unsuitable. Consequently, we applied the Spearman rank correlation test for all subsequent analyses to accurately assess the relationships between variables without assuming linearity.

M <- cor(results_pcr[results_pcr$Group == "WT", c(1:8)], method = "spearman")
M[abs(M) < 0.5] <- 0
valid_genes <- intersect(DA_list, rownames(M))
M <- M[valid_genes, valid_genes]

res_p <- cor.mtest(results_pcr[results_pcr$Group == "WT", c(1:8)], conf.level = 0.95)
res_p$p <- res_p$p[valid_genes, valid_genes]
res_p$p[abs(M) < 0.5] <- 1


corrplot(M,
method = "color",
p.mat = res_p$p,
sig.level = 0.05,
insig = "blank",
addCoef.col = "black",
number.digits = 1,
number.cex = 0.7,
tl.col = "black",
cl.lim = c(-1, 1),
col = colorRampPalette(c("magenta4", "white", "green4"))(200),
title = "Coexpression pattern in WT rats cerebellum",
mar = c(0, 0, 3, 0))

M <- cor(results_pcr[results_pcr$Group == "Het", c(1:8)], method = "spearman")
M[abs(M) < 0.5] <- 0
valid_genes <- intersect(DA_list, rownames(M))
M <- M[valid_genes, valid_genes]

res_p <- cor.mtest(results_pcr[results_pcr$Group == "Het", c(1:8)], conf.level = 0.95)
res_p$p <- res_p$p[valid_genes, valid_genes]
res_p$p[abs(M) < 0.5] <- 1

corrplot(M,
method = "color",
p.mat = res_p$p,
sig.level = 0.05,
insig = "blank",
addCoef.col = "black",
number.digits = 1,
number.cex = 0.7,
tl.col = "black",
cl.lim = c(-1, 1),
col = colorRampPalette(c("magenta4", "white", "green4"))(200),
title = "Coexpression pattern in DAT heterozygous knock-out rats cerebellum",
mar = c(0, 0, 3, 0))

M <- cor(results_pcr[results_pcr$Group == "KO", c(1:8)], method = "spearman")
M[abs(M) < 0.5] <- 0
valid_genes <- intersect(DA_list, rownames(M))
M <- M[valid_genes, valid_genes]

res_p <- cor.mtest(results_pcr[results_pcr$Group == "KO" , c(1:8)], conf.level = 0.95)
res_p$p <- res_p$p[valid_genes, valid_genes]
res_p$p[abs(M) < 0.5] <- 1

corrplot(M,
method = "color",
p.mat = res_p$p,
sig.level = 0.05,
insig = "blank",
addCoef.col = "black",
number.digits = 1,
number.cex = 0.7,
tl.col = "black",
cl.lim = c(-1, 1),
col = colorRampPalette(c("magenta4", "white", "green4"))(200),
title = "Coexpression pattern in DAT knock-out rats cerebellum",
mar = c(0, 0, 3, 0))