基于Seurat的空转单样本数据分析流程学习(四)-hdWGCNA

rm(list = ls())
library(Seurat)
library(tidyverse)
library(cowplot)
library(patchwork)
library(WGCNA)
library(hdWGCNA)
enableWGCNAThreads(nThreads = 8)

dir.create("0-hdWGCNA")   
setwd("0-hdWGCNA")   

getwd()
data <- Read10X_h5("./GSE281978/GSM8633891/filtered_feature_bc_matrix.h5")
dim(data)
object <- CreateSeuratObject(counts = data, 
                             assay = "Spatial", 
                             min.cells = 3, #过滤在少于3个细胞中表达的基因，以减少低表达基因的干扰
                             project = "GSM8633891")
object

# 再读取
image <- Read10X_Image(image.dir = "./GSE281978/GSM8633891/spatial/", 
                       image.name = "tissue_hires_image.png",
                       filter.matrix = TRUE)
image
  
dim(image)
image <- image[Cells(x = object)]# 筛选图像对象中包含的SPOT
DefaultAssay(object = image) <- "Spatial" # 设置图像对象的默认assay为"Spatial"

# 添加图片到object中
object[["GSM8633891"]] <- image
object@images[["GSM8633891"]]@scale.factors$lowres <- object@images[["GSM8633891"]]@scale.factors$hires

# scale.factors
object@images[["GSM8633891"]]@scale.factors

seurat_obj <- object
## 简单探索一下数据结构
as.data.frame(seurat_obj[["Spatial"]]$counts[1:4,1:4])
as.data.frame(LayerData(seurat_obj, 
                        assay = "Spatial", 
                        layer = "counts")[1:5,1:5])
head(seurat_obj@meta.data)
table(seurat_obj$orig.ident)

Layers(seurat_obj)
Assays(seurat_obj)

# 添加坐标信息
# 增加barcode列信息给Seurat obj 
seurat_obj$barcode <- colnames(seurat_obj)

# 读取样本路径
tissue_positions <- read.csv('./GSE281978/GSM8633891/tissue_positions_list.csv',header = T)
tissue_positions <- subset(tissue_positions, barcode %in% seurat_obj$barcode)

# 将image_df与Seurat的元数据合并
new_meta <- dplyr::left_join(seurat_obj@meta.data, tissue_positions, by='barcode')

# 将新的元数据添加到Seurat对象seurat_vhd
seurat_obj$row <- new_meta$array_row 
seurat_obj$imagerow <- new_meta$pxl_row_in_fullres
seurat_obj$col <- new_meta$array_col
seurat_obj$imagecol <- new_meta$pxl_col_in_fullres

library(qs)
qsave(seurat_obj, file="./GSE281978/GSM8633891/GSM8633891.qs")

# load the samples
dat1 <- qread("./GSE281978/GSM8633891/GSM8633891.qs")
dat1$region <- 'GSM8633891'
dat2 <- qread("./GSE281978/GSM8633892/GSM8633892.qs")
dat2$region <- 'GSM8633892'

# merge
seurat_obj <- merge(dat1, dat2)
seurat_obj$region <- factor(as.character(seurat_obj$region), 
                            levels=c('GSM8633891', 'GSM8633892'))

#joinlayers
seurat_obj <- JoinLayers(seurat_obj)

qsave(seurat_obj,file = "seurat_obj.qs")

# normalization, feature selection, scaling, and PCA
seurat_obj <- seurat_obj %>%
  NormalizeData() %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA()

# Louvain clustering and umap
seurat_obj <- FindNeighbors(seurat_obj, dims = 1:30)
seurat_obj <- FindClusters(seurat_obj,verbose = TRUE)
seurat_obj <- RunUMAP(seurat_obj, dims = 1:30)

# 设定样本分组
seurat_obj$region <- factor(as.character(seurat_obj$orig.ident),
                             levels=c('GSM8633891', 'GSM8633892'))

# show the UMAP
p1 <- DimPlot(seurat_obj, 
              label=TRUE, 
              reduction = "umap", 
              group.by = "seurat_clusters") + NoLegend()
p1

p2 <- SpatialDimPlot(seurat_obj, 
                     label = TRUE, 
                     label.size = 3)
p2

# 把准备好的注释添加到Seurat对象中
annotations <- read.csv('annotations.csv',sep = ",") 
head(annotations)
#   seurat_clusters  annotation
# 1               0 fibroblasts
# 2               1 fibroblasts
# 3               2 endothelial
# 4               3   malignant
# 5               4   dendritic
# 6               5   malignant
ix <- match(seurat_obj$seurat_clusters, annotations$seurat_clusters)
seurat_obj$annotation <- annotations$annotation[ix]
a <- seurat_obj@meta.data
colnames(seurat_obj@meta.data)
# 设定注释
Idents(seurat_obj) <- seurat_obj$annotation

p3 <- SpatialDimPlot(seurat_obj, label = TRUE, label.size = 3)
p3 + NoLegend()

seurat_obj <- SetupForWGCNA(
  seurat_obj,
  gene_select = "fraction",
  fraction = 0.05,
  wgcna_name = "vis"
)

seurat_obj <- MetaspotsByGroups(
  seurat_obj,
  group.by = c("region"),
  ident.group = "region",
  assay = 'Spatial'
)
seurat_obj  <- NormalizeMetacells(seurat_obj)

# set up the expression matrix, set group.by and group_name to NULL to include all spots
seurat_obj  <- SetDatExpr(
  seurat_obj,
  group.by=NULL,
  group_name = NULL
)

# test different soft power thresholds
library(doParallel)
registerDoParallel(1)  # 单线程
seurat_obj <- TestSoftPowers(seurat_obj)
plot_list <- PlotSoftPowers(seurat_obj)

wrap_plots(plot_list, ncol=2)

# construct co-expression network:
seurat_obj <- ConstructNetwork(
  seurat_obj,
  tom_name='test',
  overwrite_tom=TRUE
)

# plot the dendrogram
PlotDendrogram(seurat_obj, main='Spatial hdWGCNA dendrogram')

seurat_obj <- ModuleEigengenes(seurat_obj)
seurat_obj <- ModuleConnectivity(seurat_obj)

seurat_obj <- ResetModuleNames(
  seurat_obj,
  new_name = "SM"
)

# get module eigengenes and gene-module assignment tables
MEs <- GetMEs(seurat_obj)
modules <- GetModules(seurat_obj)
mods <- levels(modules$module); mods <- mods[mods != 'grey']

# add the MEs to the seurat metadata so we can plot it with Seurat functions
seurat_obj@meta.data <- cbind(seurat_obj@meta.data, MEs)

# plot with Seurat's DotPlot function
p <- DotPlot(seurat_obj, features=mods, group.by = 'annotation', 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

# 图太大了
#p1 <- SpatialFeaturePlot(
#  seurat_obj,
#  features = mods,
#  alpha = c(0.1, 1),
#  ncol = 8
#)

#p1

options(future.globals.maxSize = 50 * 1024^3)  
# perform UMAP embedding on the co-expression network
seurat_obj <- RunModuleUMAP(
  seurat_obj,
  n_hubs = 5,
  n_neighbors=15,
  min_dist=0.3,
  spread=1
)

# make the network plot
ModuleUMAPPlot(
  seurat_obj,
  edge.alpha=0.5,
  sample_edges=TRUE,
  keep_grey_edges=FALSE,
  edge_prop=0.075, 
  label_hubs=5 
)