Supplementary_for_GSE144277
Anastasia D. Belskaya, Zoia S. Fesenko, Anna B. Volnova, Raul G. Gainetdinov, and Anastasia N. Vaganova
2026-04-23
library(corrplot)
library(edgeR)
library(BioNERO)
library(SummarizedExperiment)
library(org.Mm.eg.db)
library(clusterProfiler)
# Lists for genes of interest with Ensembl codes and gene names.
DA_list <- mapIds(org.Mm.eg.db, c("Th", "Ddc", "Drd1", "Drd2", "Drd3", "Drd4","Drd5", "Slc6a3", "Slc18a2", "Comt", "Maoa", "Maob"), "ENSEMBL", "SYMBOL")
DA_list_s <- c("Th", "Ddc", "Drd1", "Drd2", "Drd3", "Drd4","Drd5", "Slc6a3", "Slc18a2", "Comt", "Maoa", "Maob")Firstly, we downloaded the table of raw counts from GEO and identified DEGs between the two groups. Since the original dataset included multiple brain structures, only the cerebellum data was extracted for further analysis. DEG analysis was performed as described in the Materials and Methods section.
Subsequently, the data were CPM-normalized for basic correlation analysis and WGCNA.
GSE144277_all.cnts.txt <- read.delim("~/Documents/RSF/Cerebel_2026/GSE144277_all.cnts.txt.gz")
GSE144277_groups <- c("Model", "Model", "Model", "WT", "WT", "WT", "WT", "WT", "Model", "Model",
"Model", "Model", "Model", "WT", "WT", "Model", "Model", "WT")
GSE144277_cer <- GSE144277_all.cnts.txt[, c(4, 5, 11, 12, 17, 18, 19, 24, 25, 26, 31, 32, 38, 39, 45, 46, 51, 52)]
#DEG analysis with edgeR
y <- DGEList(counts=GSE144277_cer, group=GSE144277_groups)
design <- model.matrix(~0+GSE144277_groups)
keep <- filterByExpr(y)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
y <- estimateDisp(y, design, robust=TRUE)
plotBCV(y)
con <- makeContrasts(GSE144277_groupsModel - GSE144277_groupsWT, levels=design)
fit <- glmQLFit(y, design)
qlf <- glmQLFTest(fit, contrast = con)
# Here we could demonstratem that expression levels of dopaminergic genes is not significantly cahanged by the ASD condition modeling
topTags(qlf, n = Inf)$table[rownames(topTags(qlf, n = Inf)$table) %in% DA_list, ]## logFC logCPM F PValue FDR
## ENSMUSG00000039358 -0.72191580 -0.3154121 2.128338084 0.1612350 0.9999851
## ENSMUSG00000040147 -0.14363158 4.3857926 1.438033873 0.2453805 0.9999851
## ENSMUSG00000000214 0.26519504 2.6788479 1.278805960 0.2723681 0.9999851
## ENSMUSG00000000326 0.06884957 6.2666057 0.317689093 0.5796733 0.9999851
## ENSMUSG00000025037 0.02239864 5.9588141 0.103425814 0.7513110 0.9999851
## ENSMUSG00000021478 -0.07366686 -1.1067223 0.037627981 0.8483206 0.9999851
## ENSMUSG00000020182 -0.03095362 2.6647168 0.010424696 0.9197594 0.9999851
## ENSMUSG00000032259 0.07288081 0.7988088 0.006370803 0.9372284 0.9999851plotMDS(y)
# CPM-normalization and filtration
GSE144277_cer <- cpm(GSE144277_cer)
is_expressed <- rowMeans(GSE144277_cer > 0.1) > 0.5
GSE144277_cer <- GSE144277_cer[is_expressed, ]Correlation analysis was performed using Spearman’s rank correlation test for the WT and model groups separately.
M <- cor(t(GSE144277_cer[rownames(GSE144277_cer) %in% DA_list, GSE144277_groups == "WT"]), method = "spearman")
# We filtered our results for correlation level
M[abs(M) < 0.6] <- 0
valid_genes <- intersect(DA_list, rownames(M))
M <- M[valid_genes, valid_genes]
rownames(M) <- mapIds(org.Mm.eg.db, rownames(M), "SYMBOL", "ENSEMBL")
colnames(M) <- mapIds(org.Mm.eg.db, colnames(M), "SYMBOL", "ENSEMBL")
res_p <- cor.mtest(t(GSE144277_cer[rownames(GSE144277_cer) %in% DA_list, GSE144277_groups == "WT"]), conf.level = 0.95)
res_p$p[valid_genes, valid_genes]## ENSMUSG00000000214 ENSMUSG00000020182 ENSMUSG00000021478
## ENSMUSG00000000214 0.00000000 0.3330551669 0.96298254
## ENSMUSG00000020182 0.33305517 0.0000000000 0.10639799
## ENSMUSG00000021478 0.96298254 0.1063979893 0.00000000
## ENSMUSG00000032259 0.34674082 0.0007389596 0.04090084
## ENSMUSG00000025496 0.10104275 0.0518337071 0.55269468
## ENSMUSG00000039358 0.24962771 0.0002931176 0.07804085
## ENSMUSG00000025094 0.18469052 0.0037394178 0.07724370
## ENSMUSG00000000326 0.20519177 0.9061271244 0.55586242
## ENSMUSG00000025037 0.27971729 0.6387768810 0.52238536
## ENSMUSG00000040147 0.03584324 0.5872431432 0.70182944
## ENSMUSG00000032259 ENSMUSG00000025496 ENSMUSG00000039358
## ENSMUSG00000000214 3.467408e-01 0.10104275 2.496277e-01
## ENSMUSG00000020182 7.389596e-04 0.05183371 2.931176e-04
## ENSMUSG00000021478 4.090084e-02 0.55269468 7.804085e-02
## ENSMUSG00000032259 0.000000e+00 0.07886215 1.159278e-06
## ENSMUSG00000025496 7.886215e-02 0.00000000 5.720124e-02
## ENSMUSG00000039358 1.159278e-06 0.05720124 0.000000e+00
## ENSMUSG00000025094 1.972839e-04 0.04901729 1.831969e-04
## ENSMUSG00000000326 4.466022e-01 0.26708803 4.740369e-01
## ENSMUSG00000025037 8.726564e-01 0.96647408 8.252746e-01
## ENSMUSG00000040147 4.644548e-01 0.17093584 4.176708e-01
## ENSMUSG00000025094 ENSMUSG00000000326 ENSMUSG00000025037
## ENSMUSG00000000214 0.1846905243 0.2051918 0.2797173
## ENSMUSG00000020182 0.0037394178 0.9061271 0.6387769
## ENSMUSG00000021478 0.0772437023 0.5558624 0.5223854
## ENSMUSG00000032259 0.0001972839 0.4466022 0.8726564
## ENSMUSG00000025496 0.0490172872 0.2670880 0.9664741
## ENSMUSG00000039358 0.0001831969 0.4740369 0.8252746
## ENSMUSG00000025094 0.0000000000 0.2627930 0.7727422
## ENSMUSG00000000326 0.2627930093 0.0000000 0.1369759
## ENSMUSG00000025037 0.7727421587 0.1369759 0.0000000
## ENSMUSG00000040147 0.1668913416 0.0337900 0.0620423
## ENSMUSG00000040147
## ENSMUSG00000000214 0.03584324
## ENSMUSG00000020182 0.58724314
## ENSMUSG00000021478 0.70182944
## ENSMUSG00000032259 0.46445475
## ENSMUSG00000025496 0.17093584
## ENSMUSG00000039358 0.41767078
## ENSMUSG00000025094 0.16689134
## ENSMUSG00000000326 0.03379000
## ENSMUSG00000025037 0.06204230
## ENSMUSG00000040147 0.00000000res_p$p[abs(M) < 0.6] <- 1
rownames(res_p$p) <- mapIds(org.Mm.eg.db, rownames(res_p$p), "SYMBOL", "ENSEMBL")
colnames(res_p$p) <- mapIds(org.Mm.eg.db, colnames(res_p$p), "SYMBOL", "ENSEMBL")
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 mice",
mar = c(0, 0, 3, 0))
M <- cor(t(GSE144277_cer[rownames(GSE144277_cer) %in% DA_list, GSE144277_groups == "Model"]), method = "spearman")
M[abs(M) < 0.6] <- 0
valid_genes <- intersect(DA_list, rownames(M))
M <- M[valid_genes, valid_genes]
rownames(M) <- mapIds(org.Mm.eg.db, rownames(M), "SYMBOL", "ENSEMBL")
colnames(M) <- mapIds(org.Mm.eg.db, colnames(M), "SYMBOL", "ENSEMBL")
res_p <- cor.mtest(t(GSE144277_cer[rownames(GSE144277_cer) %in% DA_list, GSE144277_groups == "Model"]), conf.level = 0.95)
res_p$p[valid_genes, valid_genes]## ENSMUSG00000000214 ENSMUSG00000020182 ENSMUSG00000021478
## ENSMUSG00000000214 0.00000000 3.002910e-01 0.0768001137
## ENSMUSG00000020182 0.30029101 0.000000e+00 0.0005939673
## ENSMUSG00000021478 0.07680011 5.939673e-04 0.0000000000
## ENSMUSG00000032259 0.12304206 9.109913e-06 0.0002674862
## ENSMUSG00000025496 0.22339504 1.826296e-01 0.1268529078
## ENSMUSG00000039358 0.08009416 4.625905e-05 0.0003519500
## ENSMUSG00000025094 0.28646794 6.373282e-08 0.0005209141
## ENSMUSG00000000326 0.35916992 9.313314e-01 0.4569627499
## ENSMUSG00000025037 0.60796293 7.669106e-01 0.9976664713
## ENSMUSG00000040147 0.22636971 6.442128e-03 0.0364386711
## ENSMUSG00000032259 ENSMUSG00000025496 ENSMUSG00000039358
## ENSMUSG00000000214 1.230421e-01 0.22339504 8.009416e-02
## ENSMUSG00000020182 9.109913e-06 0.18262960 4.625905e-05
## ENSMUSG00000021478 2.674862e-04 0.12685291 3.519500e-04
## ENSMUSG00000032259 0.000000e+00 0.11921004 2.042868e-08
## ENSMUSG00000025496 1.192100e-01 0.00000000 1.063119e-01
## ENSMUSG00000039358 2.042868e-08 0.10631192 0.000000e+00
## ENSMUSG00000025094 1.651606e-06 0.16824250 2.958162e-05
## ENSMUSG00000000326 5.461331e-01 0.06255606 5.859975e-01
## ENSMUSG00000025037 5.239022e-01 0.42554029 5.136382e-01
## ENSMUSG00000040147 7.884064e-03 0.09664545 1.396015e-02
## ENSMUSG00000025094 ENSMUSG00000000326 ENSMUSG00000025037
## ENSMUSG00000000214 2.864679e-01 0.35916992 0.6079629
## ENSMUSG00000020182 6.373282e-08 0.93133143 0.7669106
## ENSMUSG00000021478 5.209141e-04 0.45696275 0.9976665
## ENSMUSG00000032259 1.651606e-06 0.54613309 0.5239022
## ENSMUSG00000025496 1.682425e-01 0.06255606 0.4255403
## ENSMUSG00000039358 2.958162e-05 0.58599752 0.5136382
## ENSMUSG00000025094 0.000000e+00 0.73029490 0.7134279
## ENSMUSG00000000326 7.302949e-01 0.00000000 0.3657311
## ENSMUSG00000025037 7.134279e-01 0.36573109 0.0000000
## ENSMUSG00000040147 5.878397e-03 0.59706127 0.1446239
## ENSMUSG00000040147
## ENSMUSG00000000214 0.226369713
## ENSMUSG00000020182 0.006442128
## ENSMUSG00000021478 0.036438671
## ENSMUSG00000032259 0.007884064
## ENSMUSG00000025496 0.096645451
## ENSMUSG00000039358 0.013960152
## ENSMUSG00000025094 0.005878397
## ENSMUSG00000000326 0.597061272
## ENSMUSG00000025037 0.144623896
## ENSMUSG00000040147 0.000000000res_p$p[abs(M) < 0.6] <- 1
rownames(res_p$p) <- mapIds(org.Mm.eg.db, rownames(res_p$p), "SYMBOL", "ENSEMBL")
colnames(res_p$p) <- mapIds(org.Mm.eg.db, colnames(res_p$p), "SYMBOL", "ENSEMBL")
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 model mice",
mar = c(0, 0, 3, 0))
WGCNA was performed to identify dopaminergic genes involved in co-expression modules.
GSE144277_se <- SummarizedExperiment(assays=list(counts=cpm(GSE144277_cer)),colData=data.frame(Group = GSE144277_groups, row.names = colnames(GSE144277_cer)))
GSE144277_exp_filt <- replace_na(GSE144277_se)
GSE144277_exp_filt <- remove_nonexp(GSE144277_exp_filt, method = "median", min_exp = 0.5)
GSE144277_exp_filt <- filter_by_variance(GSE144277_exp_filt, n = 5000)
GSE144277_exp_filt@NAMES[GSE144277_exp_filt@NAMES %in% DA_list]## [1] "ENSMUSG00000000326"GSE144277_exp_filt <- ZKfiltering(GSE144277_exp_filt, cor_method = "pearson")
GSE144277_sft <- SFT_fit(GSE144277_exp_filt, net_type = "signed hybrid", cor_method = "pearson")## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 3 0.9010 1.0300 0.924 710.0 710.0 1130
## 2 4 0.7650 0.6260 0.886 570.0 544.0 991
## 3 5 0.4100 0.3210 0.647 472.0 454.0 885
## 4 6 0.0725 0.1090 0.489 400.0 383.0 801
## 5 7 0.0156 -0.0507 0.441 345.0 328.0 732
## 6 8 0.1440 -0.1730 0.517 302.0 284.0 674
## 7 9 0.2960 -0.2790 0.586 266.0 248.0 625
## 8 10 0.4110 -0.3680 0.648 237.0 216.0 581
## 9 11 0.4880 -0.4440 0.681 213.0 190.0 543
## 10 12 0.5510 -0.5090 0.710 192.0 169.0 509
## 11 13 0.5970 -0.5630 0.735 175.0 151.0 479
## 12 14 0.6370 -0.6020 0.765 159.0 135.0 451
## 13 15 0.6650 -0.6400 0.787 146.0 121.0 426
## 14 16 0.6930 -0.6730 0.806 134.0 109.0 403
## 15 17 0.7160 -0.7020 0.820 124.0 98.7 383
## 16 18 0.7240 -0.7350 0.818 115.0 89.9 363
## 17 19 0.7410 -0.7580 0.829 106.0 82.0 346
## 18 20 0.7540 -0.7780 0.842 99.1 74.4 329GSE144277_power <- GSE144277_sft$power
GSE144277_net <- exp2gcn(
GSE144277_exp_filt, net_type = "signed hybrid", SFTpower = GSE144277_power,
cor_method = "pearson"
)## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.GSE144277_MEtrait <- module_trait_cor(exp = GSE144277_exp_filt, MEs = GSE144277_net$MEs)
GSE144277_hubs <- get_hubs_gcn(GSE144277_exp_filt, GSE144277_net)
GSE144277_hubs[GSE144277_hubs$Gene %in% DA_list,]## [1] Gene Module kWithin
## <0 rows> (or 0-length row.names)GSE144277_net$genes_and_modules[GSE144277_net$genes_and_modules$Genes %in% DA_list, ]## Genes Modules
## 4324 ENSMUSG00000000326 blackplot_module_trait_cor(GSE144277_MEtrait)
We demonstrated that Comt and Maoa were* was involved in a co-expression module that was not associated with either genotype. Subsequently, we analyzed the two genotypes separately to more precisely identify genotype-specific features.
GSE144277_exp_ctrl <- replace_na(GSE144277_se[,GSE144277_se$Group == "WT"])
new_symbols <- mapIds(org.Mm.eg.db, rownames(GSE144277_exp_ctrl), "SYMBOL", "ENSEMBL")
rownames(GSE144277_exp_ctrl) <- new_symbols
#We renamed genes and applied symbols for more understandable results.
GSE144277_exp_ctrl <- GSE144277_exp_ctrl[!is.na(rownames(GSE144277_exp_ctrl)), ]
GSE144277_exp_ctrl <- GSE144277_exp_ctrl[!duplicated(rownames(GSE144277_exp_ctrl)), ]
GSE144277_exp_ctrl <- remove_nonexp(GSE144277_exp_ctrl, method = "median", min_exp = 0.5)
# Here we applied low threshold for expression variability.
GSE144277_exp_ctrl <- filter_by_variance(GSE144277_exp_ctrl, n = 5000)
GSE144277_exp_ctrl <- ZKfiltering(GSE144277_exp_ctrl, cor_method = "pearson")
GSE144277_sft <- SFT_fit(GSE144277_exp_ctrl, net_type = "signed hybrid", cor_method = "pearson")## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 3 0.89100 1.5300 0.861 794 792 1170
## 2 4 0.82800 1.1000 0.792 653 645 1020
## 3 5 0.74700 0.8250 0.721 554 556 915
## 4 6 0.64100 0.6270 0.654 479 484 828
## 5 7 0.49700 0.4690 0.559 421 425 757
## 6 8 0.36200 0.3470 0.488 375 377 697
## 7 9 0.24500 0.2500 0.489 336 338 646
## 8 10 0.12800 0.1640 0.506 305 304 602
## 9 11 0.03570 0.0840 0.489 278 276 564
## 10 12 0.00199 0.0195 0.506 254 252 530
## 11 13 0.00713 -0.0360 0.550 234 230 500
## 12 14 0.04820 -0.0961 0.569 217 211 473
## 13 15 0.11300 -0.1610 0.596 201 194 451
## 14 16 0.16700 -0.2210 0.614 188 179 431
## 15 17 0.22400 -0.2660 0.664 175 165 412
## 16 18 0.27000 -0.3180 0.681 164 154 395
## 17 19 0.31300 -0.3580 0.709 155 143 380
## 18 20 0.34900 -0.3980 0.726 146 134 365power <- GSE144277_sft$power
GSE144277_net_ctrl <- exp2gcn(
GSE144277_exp_ctrl, net_type = "signed hybrid", SFTpower = power,
cor_method = "pearson"
)## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.hubs <- get_hubs_gcn(GSE144277_exp_ctrl, GSE144277_net_ctrl)
# For checking the module in which Comt is involved. It is "darkorange"
GSE144277_net_ctrl$genes_and_modules[GSE144277_net_ctrl$genes_and_modules$Genes %in% DA_list_s, ]## Genes Modules
## 4339 Comt darkorange#Basic network characteristics.
module_genes <- GSE144277_net_ctrl$genes_and_modules$Genes[GSE144277_net_ctrl$genes_and_modules$Modules == "darkorange"]
module_adj <- GSE144277_net_ctrl$adjacency[module_genes, module_genes]
mod_density <- sum(module_adj) / (nrow(module_adj) * (ncol(module_adj) - 1))
mod_cc <- mean(WGCNA::clusterCoef(module_adj))
mod_connectivity <- mean(rowSums(module_adj))
cat("Genes count:", length(module_genes), "\n",
"Density:", mod_density, "\n",
"Clustering Coeff:", mod_cc, "\n",
"Mean Connectivity:", mod_connectivity)## Genes count: 980
## Density: 0.5479507
## Clustering Coeff: 0.6409519
## Mean Connectivity: 536.4438# We also modified hubs object to add genes of interest in the list for the visualization
manual_hubs <- rbind(hubs, data.frame(Gene = "Comt", Module = c("darkorange"), kWithin = 900))
GSE144277_net_ctrl$genes_and_modules[GSE144277_net_ctrl$genes_and_modules$Genes %in% DA_list_s, ]## Genes Modules
## 4339 Comt darkorangeGSE144277_ctrl_edges_filtered_to <- get_edge_list(
GSE144277_net_ctrl, module = c("darkorange"),
filter = TRUE, method = "pvalue",
nSamples = ncol(GSE144277_exp_ctrl)
)
custom_annotation_df <- data.frame(
Genes = unique(c(GSE144277_ctrl_edges_filtered_to$Gene1, GSE144277_ctrl_edges_filtered_to$Gene2)),
Status = "forestgreen" #this color will be for "WT" here and in the another dataset.
)
plot_gcn(
edgelist_gcn =GSE144277_ctrl_edges_filtered_to[GSE144277_ctrl_edges_filtered_to$Weight > 0.95, ],
net = GSE144277_net_ctrl,
color_by = custom_annotation_df,
hubs = hubs
)
plot_gcn(
edgelist_gcn =GSE144277_ctrl_edges_filtered_to[GSE144277_ctrl_edges_filtered_to$Weight > 0.5 & GSE144277_ctrl_edges_filtered_to$Gene1 %in% c(DA_list_s) | GSE144277_ctrl_edges_filtered_to$Gene2 %in% c(DA_list_s) , ],
net = GSE144277_net_ctrl,
color_by = custom_annotation_df,
hubs = manual_hubs[c(1:100, 494),], top_n_hubs = 22)
GSE144277_exp_Model <- replace_na(GSE144277_se[,GSE144277_se$Group == "Model"])
#We renamed genes and applied symbols for more understandable results.
new_symbols <- mapIds(org.Mm.eg.db, rownames(GSE144277_exp_Model), "SYMBOL", "ENSEMBL")
rownames(GSE144277_exp_Model) <- new_symbols
GSE144277_exp_Model <- GSE144277_exp_Model[!is.na(rownames(GSE144277_exp_Model)), ]
GSE144277_exp_Model <- GSE144277_exp_Model[!duplicated(rownames(GSE144277_exp_Model)), ]
GSE144277_exp_Model <- remove_nonexp(GSE144277_exp_Model, method = "median", min_exp = 0.5)
# Here we applied low threshold for expression variability.
GSE144277_exp_Model <- filter_by_variance(GSE144277_exp_Model, n = 5000)
GSE144277_exp_Model@NAMES[GSE144277_exp_Model@NAMES %in% DA_list]## character(0)GSE144277_exp_Model <- ZKfiltering(GSE144277_exp_Model, cor_method = "pearson")
GSE144277_sft <- SFT_fit(GSE144277_exp_Model, net_type = "signed hybrid", cor_method = "pearson")## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 3 0.5110 0.5040 0.436 861 983 1390
## 2 4 0.2040 0.2230 0.260 715 808 1250
## 3 5 0.0077 0.0365 0.238 611 674 1150
## 4 6 0.0465 -0.0887 0.296 532 573 1060
## 5 7 0.1580 -0.1740 0.362 470 494 985
## 6 8 0.2620 -0.2390 0.392 420 428 922
## 7 9 0.3500 -0.2930 0.415 378 374 867
## 8 10 0.4130 -0.3390 0.435 344 329 818
## 9 11 0.4620 -0.3750 0.439 314 290 774
## 10 12 0.5170 -0.4080 0.467 289 260 735
## 11 13 0.5540 -0.4420 0.489 266 232 700
## 12 14 0.5740 -0.4720 0.489 247 210 668
## 13 15 0.6100 -0.4960 0.525 230 188 638
## 14 16 0.6410 -0.5140 0.553 214 171 611
## 15 17 0.6610 -0.5340 0.574 201 155 586
## 16 18 0.6750 -0.5520 0.588 188 142 563
## 17 19 0.7030 -0.5650 0.621 177 130 541
## 18 20 0.7120 -0.5880 0.630 167 118 522power <- GSE144277_sft$power
GSE144277_net_Model <- exp2gcn(GSE144277_exp_Model,
net_type = "signed",
SFTpower = 6,
cor_method = "pearson")## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.# For checking the module in which dopaminergic genes are involved. It is "cyan".
GSE144277_net_Model$genes_and_modules[GSE144277_net_Model$genes_and_modules$Genes %in% DA_list_s, ]## Genes Modules
## 4335 Comt cyan#Basic network characteristics.
module_genes <- GSE144277_net_Model$genes_and_modules$Genes[GSE144277_net_Model$genes_and_modules$Modules == "cyan"]
module_adj <- GSE144277_net_Model$adjacency[module_genes, module_genes]
mod_density <- sum(module_adj) / (nrow(module_adj) * (ncol(module_adj) - 1))
mod_cc <- mean(WGCNA::clusterCoef(module_adj))
mod_connectivity <- mean(rowSums(module_adj))
cat("Genes count:", length(module_genes), "\n",
"Density:", mod_density, "\n",
"Clustering Coeff:", mod_cc, "\n",
"Mean Connectivity:", mod_connectivity)## Genes count: 2100
## Density: 0.4825281
## Clustering Coeff: 0.5917742
## Mean Connectivity: 1012.826GSE144277_Model_edges_filtered_to <- get_edge_list(
GSE144277_net_Model, module = "cyan",
filter = TRUE, method = "pvalue",
nSamples = ncol(GSE144277_exp_Model)
)
hubs <- get_hubs_gcn(GSE144277_exp_Model, GSE144277_net_Model)
# We also modified hubs object to add genes of interest in the list for the visualization
manual_hubs <- rbind(hubs, data.frame(Gene = "Comt", Module = c("cyan"), kWithin = 900))
custom_annotation_df <- data.frame(
Genes = unique(c(GSE144277_Model_edges_filtered_to$Gene1, GSE144277_Model_edges_filtered_to$Gene2)),
Status = "purple"
)
plot_gcn(
edgelist_gcn =GSE144277_Model_edges_filtered_to[GSE144277_Model_edges_filtered_to$Weight > 0.99, ],
net = GSE144277_net_Model,
color_by = custom_annotation_df,
hubs = hubs
)
#General network topology in WT mice
plot_gcn(
edgelist_gcn =GSE144277_Model_edges_filtered_to[GSE144277_Model_edges_filtered_to$Weight > 0.5 & GSE144277_Model_edges_filtered_to$Gene1 %in% c(DA_list_s) | GSE144277_Model_edges_filtered_to$Gene2 %in% c(DA_list_s) , ],
net = GSE144277_net_Model,
color_by = custom_annotation_df,
hubs = manual_hubs[c(1:100, 500),], top_n_hubs = 22)
#Topology for network components interacting with Comt and Maoa genes in the modelFinally, we extracted genes that interact with Comt or Maoa* in the identified networks and performed a functional analysis.
Het_coex <- unique(c(GSE144277_Model_edges_filtered_to[GSE144277_Model_edges_filtered_to$Weight > 0.5 & GSE144277_Model_edges_filtered_to$Gene1 %in% c(DA_list_s) | GSE144277_Model_edges_filtered_to$Gene2 %in% c(DA_list_s) , 1], (GSE144277_Model_edges_filtered_to[GSE144277_Model_edges_filtered_to$Weight > 0.5 & GSE144277_Model_edges_filtered_to$Gene1 %in% c(DA_list_s) | GSE144277_Model_edges_filtered_to$Gene2 %in% c(DA_list_s) , 2])))
WT_coex <- unique(c(GSE144277_ctrl_edges_filtered_to[GSE144277_ctrl_edges_filtered_to$Weight > 0.5 & GSE144277_ctrl_edges_filtered_to$Gene1 %in% c(DA_list_s) | GSE144277_ctrl_edges_filtered_to$Gene2 %in% c(DA_list_s) , 1], (GSE144277_ctrl_edges_filtered_to[GSE144277_ctrl_edges_filtered_to$Weight > 0.5 & GSE144277_ctrl_edges_filtered_to$Gene1 %in% c(DA_list_s) | GSE144277_ctrl_edges_filtered_to$Gene2 %in% c(DA_list_s) , 2])))
cc <- clusterProfiler::compareCluster(list(WT = WT_coex, Model = Het_coex,
Common = intersect(Het_coex, WT_coex)),
fun=enrichGO, ont="BP", OrgDb="org.Mm.eg.db", keyType = "SYMBOL")
dotplot(cc, showCategory = 8)