-
Notifications
You must be signed in to change notification settings - Fork 0
/
3. Single-cell hdWGCNA
359 lines (275 loc) · 11 KB
/
3. Single-cell hdWGCNA
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
########## NOTE: No novel code/scripts are generated from our analysis.
########## All codes used are adopted from the corresponding packages' public tutorials/vignettes used.
library(Seurat)
library(tidyverse)
library(cowplot)
library(patchwork)
library(WGCNA)
library(hdWGCNA)
library(igraph)
theme_set(theme_cowplot())
set.seed(12345)
##### Subsetting only the neutrophils from HC, KD, and MIS-C
Idents(Neutrophils) <- "disease"
HCKDMISC <- subset(Neutrophils, idents = c("HC", "KD", "MIS-C"))
##### SET UP SEURAT OBJECT FOR WGCNA
Idents(HCKDMISC) <- "RNA_snn_res.0.1"
seurat_obj <- SetupForWGCNA(
HCKDMISC,
gene_select = "fraction", # the gene selection approach
fraction = 0.05, # fraction of cells that a gene needs to be expressed in order to be included
wgcna_name = "Neutrophil hdWGCNA" # the name of the hdWGCNA experiment
)
head(Neutrophils)
DimPlot(Neutrophils, label=TRUE) + umap_theme() + ggtitle('Neutrophils') + NoLegend()
DimPlot(Neutrophils, label=TRUE, split.by = "disease") + umap_theme() + ggtitle('Neutrophils') + NoLegend()
Idents(Neutrophils) <- "RNA_snn_res.0.1"
FeaturePlot(Neutrophils, features = c("CD177"), max.cutoff = "3") + scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral"))) # pag walang ggtitle, yung name of gene lang lalabas
head(HCKDMISC)
# construct metacells in each GROUP
seurat_obj <- MetacellsByGroups(
seurat_obj = seurat_obj,
group.by = c("cellstates", "orig.ident"), # specify the columns in [email protected] to group by
k = 25, # nearest-neighbors parameter
max_shared = 10,
ident.group = 'cellstates') # set the Idents of the metacell seurat object
# normalize metacell expression matrix:
seurat_obj <- NormalizeMetacells(seurat_obj)
# Kapag nilagay mo "Neutrophils", 6 modules ang magenerate. Kapag "Condition", eh 2 modules.
seurat_obj <- SetDatExpr(
seurat_obj,
group_name = "Neutrophils", # the name of the group of interest in the group.by column
group.by='cellstates', # the metadata column containing the cell type info. This same column should have also been used in MetacellsBydiseases
assay = 'RNA', # using RNA assay
slot = 'data' # using normalized data
)
DimPlot(HCKDMISC)
# Test different soft powers:
seurat_obj <- TestSoftPowers(
seurat_obj,
networkType = 'signed') # you can also use "unsigned" or "signed hybrid"
# plot the results:
plot_list <- PlotSoftPowers(seurat_obj)
# assemble with patchwork
wrap_plots(plot_list, ncol=2)
power_table <- GetPowerTable(seurat_obj)
head(power_table, 10)
png("power-threshold.png", width=6, height=6, res=1200, units='in')
dev.off()
# construct co-expression network:
seurat_obj <- ConstructNetwork(
seurat_obj, soft_power=8,
setDatExpr=FALSE,
tom_name = 'Run1' # name of the topoligical overlap matrix written to disk
) # overwrite_tom = TRUE...maybe need new tom, try not overwriting?
png("hdWGCNA-dendrogram-neutrophil.png", width=7, height=3, res=1200, units='in')
PlotDendrogram(seurat_obj, main='Neutrophil hdWGCNA Dendrogram')
dev.off()
# need to run ScaleData first or else harmony throws an error:
seurat_obj <- ScaleData(seurat_obj, features=VariableFeatures(seurat_obj))
# compute all MEs in the full single-cell dataset
seurat_obj <- ModuleEigengenes(
seurat_obj,
group.by.vars="orig.ident"
)
# harmonized module eigengenes:
hMEs <- GetMEs(seurat_obj)
# module eigengenes:
MEs <- GetMEs(seurat_obj, harmonized=FALSE)
# compute eigengene-based connectivity (kME):
seurat_obj <- ModuleConnectivity(
seurat_obj,
group.by = 'cellstates', group_name = "Neutrophils")
# rename the modules
seurat_obj <- ResetModuleNames(
seurat_obj,
new_name = "N-M")
#We can visualize the genes in each module ranked by kME using the PlotKMEs function.
# plot genes ranked by kME for each module
PlotKMEs(seurat_obj, ncol=2)
PlotKMEs(seurat_obj, ncol=4)
PlotKMEs(seurat_obj, ncol=6)
png("hdWGCNA-PlotKME-neutrophil.png", width=6, height=6, res=1200, units='in')
dev.off()
png("hdWGCNA-PlotKME-neutrophil.png", width=8, height=2, res=1200, units='in')
dev.off()
# Get the module assignment table:
modules <- GetModules(seurat_obj)
# show the first 6 columns:
head(modules[ ,1:6]) #head(modules[,1:6])
delete <- head(modules, 50)
write.table(delete, file = "modules-top50.txt", sep = ",", quote = FALSE, row.names = F)
save(seurat_obj, file='hdWGCNA_object.Rdata')
saveRDS(seurat_obj, file='hdWGCNA_object.rds')
load(file='hdWGCNA_object.Rdata')
# # compute gene scoring for the top 25 hub genes by kME for each module
# # with Seurat method
seurat_obj <- ModuleExprScore(
seurat_obj,
n_genes = 25,
method='Seurat'
)
############# BASIC VISUALIZATION
# make a featureplot of hMEs for each module
plot_list <- ModuleFeaturePlot(
seurat_obj,
features='hMEs', # plot the hMEs
order=TRUE # order so the points with highest hMEs are on top
)
# stitch together with patchwork
wrap_plots(plot_list, ncol=2)
wrap_plots(plot_list, ncol=4)
wrap_plots(plot_list, ncol=5)
wrap_plots(plot_list, ncol=6)
png("hdWGCNA-FeaturePlot-hME-neutrophil.png", width=8, height=2, res=1200, units='in')
dev.off()
# make a featureplot of hub scores for each module
plot_list <- ModuleFeaturePlot(
seurat_obj,
features='scores', # plot the hub gene scores
order='shuffle', # order so cells are shuffled
ucell = FALSE # depending on Seurat vs UCell for gene scoring
)
# stitch together with patchwork ### Eto yung expression sa plot
wrap_plots(plot_list, ncol=2)
wrap_plots(plot_list, ncol=4)
wrap_plots(plot_list, ncol=6)
png("hdWGCNA-FeaturePlot-hubScores-neutrophil.png", width=8, height=2, res=1200, units='in')
dev.off()
# plot module correlagram
ModuleCorrelogram(seurat_obj)
png("hdWGCNA-ModuleCorrelogram-neutrophil.png", width=4, height=4, res=1200, units='in')
dev.off()
seurat_obj
# use GGally to investigate 6 selected modules: +1 MO SA TOTAL MODULES
GGally::ggpairs(GetMEs(seurat_obj)[,c(1:2,3:4)])
GGally::ggpairs(GetMEs(seurat_obj)[,c(1:3,4:5)])
GGally::ggpairs(GetMEs(seurat_obj)[,c(1:3,4:6)])
GGally::ggpairs(GetMEs(seurat_obj)[,c(1:4,5:6)])
# get hMEs from seurat object
MEs <- GetMEs(seurat_obj, harmonized=TRUE)
mods <- colnames(MEs); mods <- mods[mods != 'grey']
head(mods)
# add hMEs to Seurat meta-data:
[email protected] <- cbind([email protected], MEs)
head(seurat_obj)
#[email protected][26:29] <- NULL
# plot with Seurat's DotPlot function
# p <- DotPlot(seurat_obj, features = mods, group.by = 'cellstates')
p <- DotPlot(seurat_obj, features = mods, group.by = 'RNA_snn_res.0.1')
p + NoLegend()
# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- p +
coord_flip() +
RotatedAxis() +
scale_color_gradient2(high='red', mid='grey95', low='blue') + NoLegend()
p
png("DotPlot-NModules-per-Subset.png", width=3.25, height=1.75, res=1200, units='in')
p
dev.off()
# plot DotPlot output
p + theme(legend.text = element_text(family = 'Helvetica', size = 5), legend.key.height = unit(0.3, "cm"), legend.key.width = unit(0.3, "cm")) +
theme(legend.title = element_text(family = 'Helvetica', size = 5), legend.key.height = unit(0.3, "cm"), legend.key.width = unit(0.3, "cm"))
png("Module-Expression-DotPlot-per-Cluster.png", width=4, height=4, res=1200, units='in')
png("Module-Expression-DotPlot-per-Cluster.png", width=7, height=2.5, res=1200, units='in')
dev.off()
png("NM1-VlnPlot.png", width=5, height=5, res=1200, units='in')
png("NM2-VlnPlot.png", width=5, height=5, res=1200, units='in')
png("NM3-VlnPlot.png", width=5, height=5, res=1200, units='in')
png("NM4-VlnPlot.png", width=5, height=5, res=1200, units='in')
# Plot INH-M4 hME using Seurat VlnPlot function
p <- VlnPlot(seurat_obj, features = 'N-M4', group.by = 'RNA_snn_res.0.1', pt.size = 0)
p
p <- VlnPlot(seurat_obj, features = 'N-M4', group.by = 'Condition', pt.size = 0)
p
p <- VlnPlot(seurat_obj, features = 'N-M1', group.by = 'disease', pt.size = 0)
p
dev.off()
# add box-and-whisker plots on top:
p <- p + geom_boxplot(width=.25, fill='white')
# change axis labels and remove legend:
p <- p + xlab('') + ylab('hME') + NoLegend()
# plot output
p
Idents(HCKDMISC) <- "RNA_snn_res.0.1"
DimPlot(HCKDMISC)
############ NETWORK VISUALIZATION
############ NETWORK VISUALIZATION
############ NETWORK VISUALIZATION
# see FOLDER named "ModuleNetworks"
ModuleNetworkPlot(seurat_obj)
ModuleNetworkPlot(seurat_obj, edge.alpha = 10, vertex.label.cex =0.75)
library(RColorBrewer)
FeaturePlot(seurat_obj, "N-M2", split.by = "disease", pt.size = 0.4) & scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral")))
FeaturePlot(seurat_obj, "N-M4", split.by = "disease", pt.size = 0.4) & scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral")))
FeaturePlot(seurat_obj, "N-M2") + scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral")))
FeaturePlot(seurat_obj, "N-M3") + scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral")))
FeaturePlot(seurat_obj, "N-M4", max.cutoff = "q90", split.by = "disease", pt.size = 0.4) & scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral")))
FeaturePlot(seurat_obj, "N-M4")
# How to get hub genes
hub_df <- GetHubGenes(seurat_obj, n_hubs = 10)
hub_df <- GetHubGenes(seurat_obj, n_hubs = 25, mods=c("N-M1"))
M1 <- hub_df$gene_name
hub_df <- GetHubGenes(seurat_obj, n_hubs = 25, mods=c("N-M2"))
M2 <- hub_df$gene_name
hub_df <- GetHubGenes(seurat_obj, n_hubs = 25, mods=c("N-M3"))
M3 <- hub_df$gene_name
hub_df <- GetHubGenes(seurat_obj, n_hubs = 25, mods=c("N-M4"))
M4 <- hub_df$gene_name
HCKDMISC
head(HCKDMISC)
VlnPlot(HCKDMISC, "N-M1", split.by = "disease")
# hubgene network
HubGeneNetworkPlot(
seurat_obj,
n_hubs = 3, n_other=5,
edge_prop = 1,
mods = 'all')
p <- HubGeneNetworkPlot(seurat_obj, return_graph=TRUE)
p
png("Combined-Modules-NetworkPlot.png", width=5, height=5, res=1200, units='in')
dev.off()
par("mar")
par(mar=c(1,1,1,1))
# hubgene network
HubGeneNetworkPlot(
seurat_obj,
n_hubs = 3, n_other=5,
edge_prop = 2,
mods = 'all',
vertex.label.cex = 5,
edge.alpha = 10
)
seurat_obj <- RunModuleUMAP(
seurat_obj,
n_hubs = 5, # number of hub genes to include for the UMAP embedding
n_neighbors=10, # neighbors parameter for UMAP
min_dist=0.1 # min distance between points in UMAP space
)
# get the hub gene UMAP table from the seurat object
umap_df <- GetModuleUMAP(seurat_obj)
head(umap_df)
# plot with ggplot
ggplot(umap_df, aes(x=UMAP1, y=UMAP2)) +
geom_point(
color=umap_df$color, # color each point by WGCNA module
size=umap_df$kME*2 # size of each point based on intramodular connectivity
) +
umap_theme() #+ DarkTheme()
ModuleUMAPPlot(
seurat_obj,
edge.alpha=0.25,
sample_edges=TRUE,
edge_prop=0.1, # proportion of edges to sample (20% here)
label_hubs=10 ,# how many hub genes to plot per module?
keep_grey_edges=FALSE)
library(ggrepel)
ModuleUMAPPlot(
seurat_obj,
edge.alpha=0.25,
sample_edges=TRUE,
edge_prop=0.2, # proportion of edges to sample (20% here)
label_hubs=2 ,# how many hub genes to plot per module?
keep_grey_edges=FALSE)#+ labs(title = "geom_text_repel()")
g <- ModuleUMAPPlot(seurat_obj, return_graph=TRUE)
rm(p, p2)