#Здесь мы запустим все пакеты, которые понадобятся позже
library(ggplot2)
library(ggpubr)
library(rstatix)
library(dplyr)
library(org.Hs.eg.db)
library(SummarizedExperiment)
library(BioNERO)
library(clusterProfiler)
library(GOSemSim)
library(edgeR)
library(psych)
library(BRETIGEA)

Сейчас мы загрузим набор данных и убедимся, что существенных различий между группами в нём нет. Это хороший пример, чтобы искать неочевидные закономерности с помощью WGCNA

urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?format=file&type=rnaseq_counts"
path <- paste(urld, "acc=GSE87194", "file=GSE87194_raw_counts_GRCh38.p13_NCBI.tsv.gz", sep="&")
GSE87194 <- as.matrix(data.table::fread(path, header=T, colClasses="integer"), rownames=1)
groups <- c(rep("Ctrl", 19), rep("Shz", 19))
y <- DGEList(counts=GSE87194,group=groups)
keep <- filterByExpr(y)
y <- y[keep,,keep.lib.sizes=FALSE]
y <- normLibSizes(y)
design <- model.matrix(~0+groups)
y <- estimateDisp(y, design, robust=TRUE)
y$common.dispersion
## [1] 0.1039908
# здесь можно посмотреть результат снижения размерностей. Разделить образцы на группы не удалось.
group_factor <- as.factor(y$samples$group) 
colors_to_use <- c("magenta3", "forestgreen")[group_factor]
plotMDS(y, 
        col = colors_to_use, 
        pch = 16,            
        cex = 2,             
        gene.selection = "common") 

#Убедимся, что действительно нет никакой разницы
con <- makeContrasts(groupsShz - groupsCtrl, levels=design)
fit <- glmQLFit(y, design)
qlf<-glmQLFTest(fit,contrast=con)
res <- topTags(qlf, n = Inf)$table
res$diffexpressed <- "NS"
res$diffexpressed[res$logFC > 1 & res$FDR < 0.05] <- "UP"
res$diffexpressed[res$logFC < -1 & res$FDR < 0.05] <- "DOWN"
ggplot(data = res, aes(x = logFC, y = -log10(FDR), col = diffexpressed)) +
geom_point(alpha = 0.4, size = 1.5) +
theme_minimal() +
scale_color_manual(values = c("DOWN" = "cyan", "NS" = "grey", "UP" = "magenta3")) +
geom_vline(xintercept = c(-1, 1), col = "black", linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), col = "black", linetype = "dashed") +
labs(title = "Volcano Plot (edgeR)", x = "log2 Fold Change", y = "-log10 FDR")

Сделаем небольшую, но важную вещь. Поскольку мы хотим видеть названия генов, нужно переименовать названия строк в таблице. Сейчас они в формате ENTREZ ID. Кроме того, таблицу нужно загрузить повторно, потому что нам понадобятся нормализованные данные. В этот раз таблица будет загружена в варианте, содержащем данные в TPM.

#Вот как называются строки и выглядят значения в таблице сейчас
head(GSE87194)
##           GSM2324383 GSM2324384 GSM2324385 GSM2324386 GSM2324387 GSM2324388
## 100287102          4          5          4          6          1          6
## 653635            85         38         84        104         39         91
## 102466751          3          1          3          3          1          4
## 107985730          3          1          4          1          3          1
## 100302278          0          0          0          0          0          0
## 645520             0          0          0          0          0          1
##           GSM2324389 GSM2324390 GSM2324391 GSM2324392 GSM2324393 GSM2324394
## 100287102          4          3          2          8          5          6
## 653635            80         62         40        241        168        122
## 102466751          3          1          1         11          3          5
## 107985730          0          1          0          3          5          3
## 100302278          0          0          0          0          0          0
## 645520             1          0          0          1          1          1
##           GSM2324395 GSM2324396 GSM2324397 GSM2324398 GSM2324399 GSM2324400
## 100287102          9          7         10          6          8          7
## 653635           156        133        128        155         99        148
## 102466751          8          7          4          5          2          2
## 107985730          3          3          0          3          4          2
## 100302278          0          0          0          0          0          0
## 645520             1          0          1          1          0          0
##           GSM2324401 GSM2324402 GSM2324403 GSM2324404 GSM2324405 GSM2324406
## 100287102          7          9          3          6          4         10
## 653635           121        102         31         65         91         91
## 102466751          5          3          1          1          2          3
## 107985730          2          0          1          0          4          4
## 100302278          0          0          0          0          0          0
## 645520             1          0          0          0          0          1
##           GSM2324407 GSM2324408 GSM2324409 GSM2324410 GSM2324411 GSM2324412
## 100287102          1          5          7          5          4          6
## 653635            25         70         63         39         63        106
## 102466751          1          2          1          0          2          5
## 107985730          0          1          4          3          0          1
## 100302278          0          0          0          0          0          1
## 645520             0          0          0          0          0          1
##           GSM2324413 GSM2324414 GSM2324415 GSM2324416 GSM2324417 GSM2324418
## 100287102          7          5         11         14          3          6
## 653635           111        131        129        167         62        126
## 102466751          2          3          6          6          3          4
## 107985730          2          1          6          3          3          5
## 100302278          0          0          0          0          0          0
## 645520             1          1          1          4          0          0
##           GSM2324419 GSM2324420
## 100287102          6          3
## 653635           114        144
## 102466751          2          4
## 107985730          4          2
## 100302278          0          0
## 645520             1          0
path <- paste(urld, "acc=GSE87194", "file=GSE87194_norm_counts_TPM_GRCh38.p13_NCBI.tsv.gz", sep="&")
GSE87194 <- as.matrix(data.table::fread(path, header=T, colClasses="integer"), rownames=1)
GSE87194 <- as.data.frame(GSE87194)
rows <- rownames(GSE87194)
rows <- mapIds(org.Hs.eg.db, rows, "SYMBOL", "ENTREZID")
## 'select()' returned 1:1 mapping between keys and columns
GSE87194 <- cbind(GSE87194, rows)
GSE87194 <- na.omit(GSE87194)
GSE87194 <- aggregate(. ~ rows, data = GSE87194, FUN = sum)
rownames(GSE87194) <- GSE87194[,1]
GSE87194 <- GSE87194[,-1]
#Вот как называются строки и выглядят показатели экспрессии в новом варианте
head(GSE87194)
##          GSM2324383 GSM2324384 GSM2324385 GSM2324386 GSM2324387 GSM2324388
## A1BG        0.61320    0.36460    0.31010     0.5211     0.4069    0.67420
## A1BG-AS1    2.17000    0.90050    0.69880     1.6840     1.1710    2.52300
## A1CF        0.08084    0.07189    0.05325     0.1189     0.0000    0.02564
## A2M        25.47000   19.12000   12.98000    25.6500    21.0100   28.70000
## A2M-AS1     3.27200    2.42200    1.38300     3.0930     2.7670    2.86800
## A2ML1       1.00400    0.65730    0.93320     0.6606     0.7479    0.79350
##          GSM2324389 GSM2324390 GSM2324391 GSM2324392 GSM2324393 GSM2324394
## A1BG        0.54370    0.34620    0.25370     0.5401    0.37280    0.66390
## A1BG-AS1    1.65500    0.89530    0.73010     2.3450    2.17100    1.80700
## A1CF        0.06432    0.09479    0.02251     0.1138    0.04812    0.06312
## A2M        34.40000   12.84000   13.13000    35.0700   32.64000   25.61000
## A2M-AS1     3.04200    1.25100    2.21500     1.1290    1.84400    1.69100
## A2ML1       0.77340    0.69220    0.82810     0.5905    0.71030    1.05200
##          GSM2324395 GSM2324396 GSM2324397 GSM2324398 GSM2324399 GSM2324400
## A1BG        0.71970    0.61950    1.02800    0.53280     0.4782     0.6369
## A1BG-AS1    1.54900    2.65200    2.60900    2.60600     1.7460     1.9760
## A1CF        0.07058    0.01912    0.07132    0.05673     0.1697     0.1356
## A2M        38.79000   17.93000   34.18000   20.22000    27.2000    28.2600
## A2M-AS1     1.86600    2.02000    1.35600    0.95590     2.6960     2.0890
## A2ML1       1.15200    1.29700    1.18100    0.55660     1.7220     1.1580
##          GSM2324401 GSM2324402 GSM2324403 GSM2324404 GSM2324405 GSM2324406
## A1BG         0.5708     0.6672    0.33690     0.1686     0.3002    0.41350
## A1BG-AS1     1.6410     1.6480    1.42700     0.2759     0.8122    1.48300
## A1CF         0.1223     0.0000    0.06644     0.1088     0.1184    0.05436
## A2M         30.9000    17.0600   23.68000     9.6330    20.6400   28.10000
## A2M-AS1      1.9250     2.1370    1.43100     2.5850     2.5020    3.28800
## A2ML1        1.2970     0.5714    0.51250     0.3524     0.5751    0.86200
##          GSM2324407 GSM2324408 GSM2324409 GSM2324410 GSM2324411 GSM2324412
## A1BG         0.6641    0.28250    0.21810    0.51820    0.48910    0.32040
## A1BG-AS1     1.1720    1.21100    0.75550    0.86390    1.17900    2.32100
## A1CF         0.1048    0.09525    0.05629    0.01839    0.05787    0.06752
## A2M         12.0000   17.19000   10.84000   13.51000   14.06000   44.16000
## A2M-AS1      2.5500    2.24300    1.44300    0.72380    2.18400    3.28300
## A2ML1        0.3368    0.93820    0.54280    0.67000    0.33990    1.09100
##          GSM2324413 GSM2324414 GSM2324415 GSM2324416 GSM2324417 GSM2324418
## A1BG        0.59540     0.4613     1.1370     0.5589    0.52850    0.69990
## A1BG-AS1    2.23800     2.2360     1.7270     2.3580    1.22100    2.29500
## A1CF        0.06038     0.1336     0.0000     0.1133    0.01137    0.06515
## A2M        32.92000    21.8400    32.8400    30.2700   11.06000   14.39000
## A2M-AS1     2.66200     1.3210     2.1750     2.7710    0.77710    1.28200
## A2ML1       1.15200     0.7078     0.9916     1.1090    0.51970    1.10500
##          GSM2324419 GSM2324420
## A1BG        0.73940    0.75510
## A1BG-AS1    1.97600    1.97800
## A1CF        0.06998    0.01005
## A2M        37.17000   20.98000
## A2M-AS1     2.84500    1.45700
## A2ML1       1.10000    1.25900

Теперь можно приступить собственно к взвешенному анализу сетей ко-экспрессии генов, WGCNA

colData <- data.frame(Group=groups, row.names = colnames(GSE87194))
GSE87194_se <- SummarizedExperiment(assays=list(counts=cpm(GSE87194)),colData=colData)
exp_filt <- replace_na(GSE87194_se)
exp_filt <- remove_nonexp(exp_filt, method = "median", min_exp = 5)
exp_filt <- filter_by_variance(exp_filt, n = 1250)
plot_heatmap(GSE87194_se, type = "samplecor", show_rownames = FALSE)

plot_heatmap(exp_filt, type = "samplecor", show_rownames = FALSE)

Мы отфильтровали для анализа 1250 генов, экспрессирующихся на уровне не ниже 5 TPM и наиболее вариабельных. Две тепловых карты показывают, что до фильтрации (1) сходство между образцами гораздо менее заметно, чем после (2).

exp_filt <- ZKfiltering(exp_filt, cor_method = "pearson")
## Number of removed samples: 2
exp_filt <- PC_correction(exp_filt)
sft <- SFT_fit(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.650 -0.640          0.561   89.80    63.700  248.0
## 2      4    0.733 -0.788          0.668   61.10    38.500  200.0
## 3      5    0.808 -0.868          0.758   44.20    24.700  166.0
## 4      6    0.810 -0.947          0.758   33.40    17.100  141.0
## 5      7    0.828 -0.993          0.780   26.00    12.200  122.0
## 6      8    0.839 -1.040          0.797   20.70     8.610  106.0
## 7      9    0.824 -1.090          0.790   16.80     6.400   92.7
## 8     10    0.843 -1.100          0.817   13.80     4.750   81.8
## 9     11    0.867 -1.120          0.854   11.50     3.570   72.6
## 10    12    0.871 -1.140          0.877    9.73     2.860   64.7
## 11    13    0.888 -1.150          0.906    8.28     2.250   57.9
## 12    14    0.898 -1.140          0.926    7.10     1.750   52.0
## 13    15    0.903 -1.160          0.933    6.14     1.420   46.8
## 14    16    0.908 -1.180          0.946    5.34     1.140   42.3
## 15    17    0.903 -1.190          0.955    4.67     0.983   38.3
## 16    18    0.904 -1.200          0.959    4.11     0.829   34.8
## 17    19    0.909 -1.210          0.966    3.64     0.701   31.6
## 18    20    0.909 -1.210          0.967    3.23     0.610   28.8
power <- sft$power
net <- exp2gcn(
exp_filt, net_type = "signed hybrid", SFTpower = power,
cor_method = "pearson"
)
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
MEtrait <- module_trait_cor(exp = exp_filt, MEs = net$MEs)
plot_module_trait_cor(MEtrait)

Три модуля, отмеченных звёздочками ассоциированы с тем или иным состоянием. Модуль включает гены и корреляции между ними. Давайте посмотрим, что они из собой представляют.

plot_expression_profile(
exp = exp_filt,
net = net,
plot_module = TRUE,
modulename = "brown"
)

plot_expression_profile(
exp = exp_filt,
net = net,
plot_module = TRUE,
modulename = "red"
)

plot_expression_profile(
exp = exp_filt,
net = net,
plot_module = TRUE,
modulename = "magenta"
)

Это профили экспрессии всех генов в модуле. В том состоянии, с которым есть ассоциация, колебания их уровня более согласованы. Ассоциации слабые, хотя и статистически значимые, так что и согласованность выражена не сильно. Жирная линия - общий тренд, усреднённые значения уровня экспрессии генов в образце. Напомню, что коэкспрессия генов кластера “brown” характерна для нормы, а остальных двух, “red” и “magenta” - для пациентов с шизофренией.

То же самое можно представить как тепловые карты (heatmaps)

module_genes <- net$genes_and_modules[net$genes_and_modules$Modules == "brown", ]$Genes
exp_module <- exp_filt[rownames(exp_filt) %in% module_genes, ]
plot_heatmap(exp_module, type = "expr", log_trans = TRUE, show_rownames = FALSE, cluster_cols = FALSE)

module_genes <- net$genes_and_modules[net$genes_and_modules$Modules == "red", ]$Genes
exp_module <- exp_filt[rownames(exp_filt) %in% module_genes, ]
plot_heatmap(exp_module, type = "expr", log_trans = TRUE, show_rownames = FALSE, cluster_cols = FALSE)

module_genes <- net$genes_and_modules[net$genes_and_modules$Modules == "magenta", ]$Genes
exp_module <- exp_filt[rownames(exp_filt) %in% module_genes, ]
plot_heatmap(exp_module, type = "expr", log_trans = TRUE, show_rownames = FALSE, cluster_cols = FALSE)

Обратите внимание, что в той группе, с которой ассоциирован тот или иной модуль, он представлен на тепловой карте несколько более однородно.

Теперь давайте визуализируем наши модули.

hubs <- get_hubs_gcn(exp_filt, net) #так мы находим гены-хабы, те гены, у которых больше всего связей.
edges_filtered_to <- get_edge_list(
net, module = c("magenta", "red", "brown"),
filter = TRUE, method = "pvalue",
nSamples = ncol(exp_filt)
)
plot_gcn(
edgelist_gcn =edges_filtered_to,
net = net,
color_by = "module",
hubs = hubs
)

plot_gcn(
edgelist_gcn =edges_filtered_to[edges_filtered_to$Weight > 0.8, ],
net = net,
color_by = "module",
hubs = hubs
)

Похоже, что у пациентов с шизофренией активируется какая-то сеть, управляемая малой некодирующей РНК SNORD115-32 и меняется активность какого-то процесса, в который вовлечены треониновые транспортные РНК. Последнее не так уж и странно, сейчас стало известно, что тРНК учатсвуют не только в синтезе белка и их активность действительльно может меняться при шизофрении.

Теперь перейдём ко второму важному блоку задач на сегодня - анализу профиля коэкспрессии.

#Подготовим два отдельных фрейма
GSE87194_shz <- data.frame(t(GSE87194[, c(20:38)]))
GSE87194_ctrl <- data.frame(t(GSE87194[, c(1:19)]))
#Мы будем изучать связи гена GRIN2A, его мутации связаны с повышенным риском шизофрении. Для этого выделим коэкспрессирующиеся гены с помощью корреляции Кендалла.
corrs <- corr.test(x = GSE87194_shz, y = GSE87194_shz$GRIN2A, use = "pairwise",method="kendall",alpha=.05,ci=TRUE,minlength=5,normal=F)
GSE87194_shz_cluster  <- cbind(corrs[["r"]], corrs[["p"]], corrs[["p.adj"]])
colnames(GSE87194_shz_cluster) <- c("corr", "p.value", "p.val_BH.adj")
GSE87194_shz_cluster <- as.data.frame(na.omit(GSE87194_shz_cluster))
GSE87194_shz_cluster <- GSE87194_shz_cluster[rev(order(GSE87194_shz_cluster$corr)),]
GSE87194_shz_cluster <- GSE87194_shz_cluster[!(row.names(GSE87194_shz_cluster)=="GRIN2A"),]
corrs <- corr.test(x = GSE87194_ctrl, y = GSE87194_ctrl$GRIN2A, use = "pairwise",method="kendall",alpha=.05,ci=TRUE,minlength=5,normal=F)
GSE87194_ctrl_cluster  <- cbind(corrs[["r"]], corrs[["p"]], corrs[["p.adj"]])
colnames(GSE87194_ctrl_cluster) <- c("corr", "p.value", "p.val_BH.adj")
GSE87194_ctrl_cluster <- as.data.frame(na.omit(GSE87194_ctrl_cluster))
GSE87194_ctrl_cluster <- GSE87194_ctrl_cluster[rev(order(GSE87194_ctrl_cluster$corr)),]
GSE87194_ctrl_cluster <- GSE87194_ctrl_cluster[!(row.names(GSE87194_ctrl_cluster)=="GRIN2A"),]
cc <- compareCluster(list(Ctrl = rownames(GSE87194_ctrl_cluster[GSE87194_ctrl_cluster$corr >0.7,]), Shz = rownames(GSE87194_shz_cluster[GSE87194_shz_cluster$corr >0.7,])),
fun=enrichGO, ont="BP", OrgDb="org.Hs.eg.db", keyType = "SYMBOL")
dotplot(cc)

Ещё одна методика оценки функциональности гена – сравнение семантического сходства генов, коэкспрессирующихся с ним. Размер кластеров получился очень разным. Для примера мы возмём 100 генов с наиболее выраженной коэкспрессией с GRIN2A. Также нам понадобится случайный набор генов.

hsGO2 <- godata('org.Hs.eg.db', keytype = "SYMBOL", ont="BP", computeIC=FALSE) 
## preparing gene to GO mapping data...
Shz_sim <- (mgeneSim(rownames(GSE87194_shz_cluster[c(1:100),]), semData=hsGO2, measure="Wang", combine="BMA", verbose=FALSE))
Ctrl_sim <- (mgeneSim(rownames(GSE87194_ctrl_cluster[c(1:100),]), semData=hsGO2, measure="Wang", combine="BMA", verbose=FALSE))
Random_sim <- (mgeneSim(sample(rownames(GSE87194_ctrl_cluster), 100), semData=hsGO2, measure="Wang", combine="BMA", verbose=FALSE))
#Это матрицы с попарным сравнением всех генов. Нам понадобится удалить повторы и сделать из матриц вектора.
Shz_sim[lower.tri(Shz_sim, diag = T)] <- NA
Ctrl_sim[lower.tri(Ctrl_sim, diag = T)] <- NA
Random_sim[lower.tri(Random_sim, diag = T)] <- NA
Shz_sim <- na.omit(c(Shz_sim))
Ctrl_sim <- na.omit(c(Ctrl_sim))
Random_sim <- na.omit(c(Random_sim))
SemSim <- data.frame(Wang = c(Ctrl_sim, Shz_sim, Random_sim), Group = c(rep("Ctrl", length(Ctrl_sim)), rep("Shz", length(Shz_sim)), rep("Random", length(Random_sim))))
stat.test <- SemSim %>%
rstatix::wilcox_test(Wang ~ Group)
stat.test$y.position <- c(1.1, 1.2, 1.3)
stat.test
## # A tibble: 3 × 10
##   .y.   group1 group2    n1    n2 statistic     p p.adj p.adj.signif y.position
##   <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <chr>             <dbl>
## 1 Wang  Ctrl   Random  3160   630  1033487  0.129 0.387 ns                  1.1
## 2 Wang  Ctrl   Shz     3160  3160  5078468  0.238 0.476 ns                  1.2
## 3 Wang  Random Shz      630  3160   974646. 0.408 0.476 ns                  1.3
ggplot(data = SemSim, aes(x = factor(Group, level = c("Ctrl", "Shz", "Random")), y = Wang))+
geom_boxplot(aes(fill=Group))+
scale_fill_manual(values = c("purple", "orchid1", "coral1"))+
theme(axis.title = element_text(size = 16), axis.text = element_text(size = 14))+
theme(legend.title = element_text(size = 16),
legend.text = element_text(size = 14))+
stat_pvalue_manual(stat.test, label = "p.adj.signif", tip.length = 0.01)+
xlab("")

На графике только в группе пациентов с шизофренией связь между генами выше, чем в случайном наборе. По-видимому, в данной группе GRIN2A действительно имеет более сложные связи, чем это характерно для нормы.

Последний тест на сегодня - деконволюция. Мы сделаем его с помощью готового пакета BRETIGEA, но можно использовать и данные других исследований единиченых клеток. Деконволюция позволяет оценить изменения клеточного состава в исследуемых образцах.

cell_proportions <- as.data.frame(brainCells(cpm(GSE87194),
species = "human",
nMarker = 50, scale = TRUE))
cell_proportions <- as.data.frame(cbind(cell_proportions, groups))
myy <- data.frame(Cells = c(rep("ast", 38), rep("end", 38),rep("mic", 38),rep("neu", 38),rep("oli", 38),rep("opc", 38)), Level = c(cell_proportions$ast, cell_proportions$end, cell_proportions$mic, cell_proportions$neu, cell_proportions$oli, cell_proportions$opc), Group = rep(cell_proportions$groups, 6))
ggplot(myy, aes(fill=Group, y=Level, x=Cells)) +
geom_boxplot()+
theme(axis.title = element_text(size = 16), axis.text = element_text(size = 14))+
scale_fill_manual(values = c("purple", "orchid1"))

Шкала Y отражает отклонение пропорции данных клеток от среднего в образцах - в данном случае, существенных изменений нет.