植物空間轉(zhuǎn)錄組分析1：Seurat基本流程 - 簡書 (jianshu.com)
植物空間轉(zhuǎn)錄組分析2：STEEL+Seurat - 簡書 (jianshu.com)
植物空間轉(zhuǎn)錄組分析 3：hdWGCNA - 簡書 (jianshu.com)
前言
時(shí)隔半年的更新，這次的主題是單細(xì)胞hdWGCNA分析，文章剛剛被接收，目前在植物空間轉(zhuǎn)錄組已經(jīng)有文章使用了該方法勺卢，基本的內(nèi)容就是分module然后去看每個(gè)module的hub基因搁宾，或者做做GO富集之類的
image.png
R包安裝以及數(shù)據(jù)準(zhǔn)備
本次分析還是以蘭花的空間轉(zhuǎn)錄組為基礎(chǔ)稀颁，具體的數(shù)據(jù)參考以前的推文
# install BiocManager
install.packages("BiocManager")
# install Bioconductor core packages
BiocManager::install()
# install additional packages:
install.packages(c("Seurat", "WGCNA", "igraph", "devtools"))
devtools::install_github('smorabit/hdWGCNA', ref='dev')
大家可能會(huì)因?yàn)榫W(wǎng)絡(luò)問題無法下載表锻，我直接將壓縮包發(fā)給大家
## 數(shù)據(jù)準(zhǔn)備
## 將spot坐標(biāo)信息加到metadata中
# make a dataframe containing the image coordinates for each sample
image_df <- do.call(rbind, lapply(names(orc1_new@images), function(x){
orc1_new@images[[x]]@coordinates
}))
# merge the image_df with the Seurat metadata
new_meta <- merge(orc1_new@meta.data, image_df, by='row.names')
# fix the row ordering to match the original seurat object
rownames(new_meta) <- new_meta$Row.names
ix <- match(as.character(colnames(orc1_new)), as.character(rownames(new_meta)))
new_meta <- new_meta[ix,]
# add the new metadata to the seurat object
orc1_new@meta.data <- new_meta
構(gòu)建metaspots
## Here we set up the data for hdWGCNA and run MetaspotsByGroups.
subdata <- subset(orc1_new, idents = c(19,21,37,38,39,40))
subdata <- SCTransform(subdata, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)
subdata <- RunPCA(subdata, features = VariableFeatures(subdata))
subdata <- SetupForWGCNA(
subdata,
features = VariableFeatures(subdata),
wgcna_name = "SCT"
)
subdata <- MetacellsByGroups(
seurat_obj = subdata,
group.by = c("seurat_clusters"),
k = 5,
max_shared= 10,
min_cells = 6,
reduction = 'pca',
ident.group = 'seurat_clusters',
slot = 'scale.data',
assay = 'SCT'
)
m_obj <- GetMetacellObject(subdata)
m_obj
An object of class Seurat
14006 features across 180 samples within 1 assay
Active assay: SCT (14006 features, 0 variable features)
這里需要注意，一般空間分析推薦使用的是 MetaspotsByGroups函數(shù)们衙，我嘗試了一遍softpower只有1，module只有兩個(gè)碱呼，所以又使用了MetacellsByGroups函數(shù)嘗試蒙挑。
共表達(dá)網(wǎng)絡(luò)分析
# set up the expression matrix, set group.by and group_name to NULL to include all spots
subdata <- SetDatExpr(
subdata,
group.by=NULL,
group_name = NULL
)
# test different soft power thresholds
subdata <- TestSoftPowers(subdata)
plot_list <- PlotSoftPowers(subdata)
wrap_plots(plot_list, ncol=2)
image.png
# construct co-expression network:
subdata <- ConstructNetwork(
subdata,
soft_power = 1,
tom_name='test',
overwrite_tom=TRUE
)
# plot the dendrogram
PlotDendrogram(subdata, main='Spatial hdWGCNA dendrogram')
image.png
compute module eigengenes
subdata <- ModuleEigengenes(subdata)
subdata <- ModuleConnectivity(subdata)
modules <- GetModules(subdata) %>% subset(module != 'grey')
head(modules[,1:3])
gene_name module color
PAXXG297280 PAXXG297280 turquoise turquoise
PAXXG074500 PAXXG074500 blue blue
PAXXG239080 PAXXG239080 turquoise turquoise
PAXXG067750 PAXXG067750 turquoise turquoise
PAXXG327430 PAXXG327430 turquoise turquoise
PAXXG311150 PAXXG311150 blue blue
模塊數(shù)目比較少也是正常的
Data visualization
tom<-GetTOM(subdata)
##Data visualization
##Here we visualize module eigengenes using the Seurat functions DotPlot and SpatialFeaturePlot.
##For network visualization, please refer to the Network visualization tutorial.
# get module eigengenes and gene-module assignment tables
MEs <- GetMEs(subdata)
modules <- GetModules(subdata)
mods <- levels(modules$module)
write.table(file="modules.xls",modules,row.names = F,sep="\t")
# plot with Seurat's DotPlot function
p <- DotPlot(subdata, features=mods, group.by = 'seurat_clusters', dot.min=0.1)
# 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') +
xlab('') + ylab('')
p
p <- SpatialFeaturePlot(
subdata,
features = mods,
ncol = 4,
pt.size.factor = 3.5
)
p
image.png
image.png
module network plots
分為兩種，一種每個(gè)模塊的hubgene愚臀，一種是全部模塊的hubgene
## Individual module network plots
ModuleNetworkPlot(subdata,
plot_size = c(10, 10),
vertex.label.cex = 0.6
)
hub_df <- GetHubGenes(subdata, n_hubs = 25)
write.table(file="hub_df.xls",hub_df,row.names = F,sep="\t")
## hubgene network
HubGeneNetworkPlot(
subdata,
n_hubs = 4, n_other=6,
edge_prop = 0.75,
mods = 'all',
hub.vertex.size = 4,
other.vertex.size = 1,
edge.alpha = 0.5
)
image.png
總結(jié)
目前只是初步上手忆蚀，其實(shí)這個(gè)方法主要的目的還是根據(jù)表達(dá)對基因進(jìn)行聚類，再加個(gè)功能富集可能會(huì)更好，后續(xù)還有UMAP to co-expression networks流程馋袜，但因?yàn)槲以谧龅臅r(shí)候效果比較差所以沒放男旗。
轉(zhuǎn)載請注明>>>周小釗的博客, 打賞推薦博客
