空间转录组对比：Squidpy/Giotto/STutility¶

# === Squidpy 标准流程 ===
import squidpy as sq
import scanpy as sc

# 读取Visium数据
adata = sq.datasets.visium_hne_adata()     # 示例数据
# 或从10x输出读取:
# adata = sc.read_visium("spaceranger_output/")

# 标准Scanpy预处理
sc.pp.normalize_total(adata)               # 标准化
sc.pp.log1p(adata)                         # log转换
sc.pp.highly_variable_genes(adata)         # 高变基因
sc.tl.pca(adata)                           # PCA
sc.pp.neighbors(adata)                     # KNN
sc.tl.leiden(adata, resolution=0.5)        # 聚类
sc.tl.umap(adata)                          # UMAP

# ★空间分析（Squidpy特有）
# 空间邻域图
sq.gr.spatial_neighbors(adata,
    coord_type="grid",                     # 网格型（Visium）
    n_neighs=6                             # 六角形邻域
)

# 空间自相关（Moran's I）
sq.gr.spatial_autocorr(adata,
    mode="moran",                          # Moran's I统计量
    genes=adata.var_names[:100]            # 前100个基因
)
# 结果在 adata.uns["moranI"]
# I值高→空间聚集 | I≈0→随机分布 | I值低→空间分散

# 邻域富集分析
sq.gr.nhood_enrichment(adata,
    cluster_key="leiden"                   # 聚类结果
)
sq.pl.nhood_enrichment(adata, cluster_key="leiden")  # 可视化

# 共现分析
sq.gr.co_occurrence(adata,
    cluster_key="leiden"                   # 哪些细胞类型倾向共同出现
)

# 空间可视化
sq.pl.spatial_scatter(adata,
    color="leiden",                        # 空间聚类图
    size=1.5
)
sq.pl.spatial_scatter(adata,
    color=["CD3D", "CD68"],               # 基因空间表达
    size=1.5
)

# === Giotto 标准流程 ===
library(Giotto)

# 创建Giotto对象（Visium数据）
my_giotto <- createGiottoObject(
  expression = expr_matrix,                # 表达矩阵
  spatial_locs = spatial_coords             # 空间坐标
)
# 或从10x Visium输出:
# my_giotto <- createGiottoVisiumObject(
#   visium_dir = "spaceranger_output/",
#   h5_visium_path = "filtered_feature_bc_matrix.h5"
# )

# 预处理
my_giotto <- normalizeGiotto(my_giotto)    # 标准化
my_giotto <- calculateHVG(my_giotto)       # 高变基因
my_giotto <- runPCA(my_giotto)             # PCA
my_giotto <- createNearestNetwork(my_giotto)  # KNN
my_giotto <- doLeidenCluster(my_giotto)    # Leiden聚类
my_giotto <- runUMAP(my_giotto)            # UMAP

# ★空间分析
# 空间网络
my_giotto <- createSpatialNetwork(my_giotto,
  method = "kNN",                          # 方法: kNN/Delaunay
  k = 6                                   # 邻居数
)

# 空间基因检测
spatial_genes <- detectSpatialCorGenes(
  my_giotto,
  method = "network",                      # 基于空间网络
  spatial_network_name = "kNN_network"
)

# 空间域检测（HMRF）
my_giotto <- doHMRF(my_giotto,
  spatial_genes = spatial_genes$gene,      # 空间基因
  spatial_network_name = "kNN_network",
  k = 5                                   # 空间域数
)

# 可视化
spatPlot(my_giotto,
  cell_color = "leiden_clus",              # 聚类着色
  point_size = 2
)

spatFeatPlot2D(my_giotto,
  feats = c("CD3D", "CD68"),              # 基因空间图
  point_size = 2
)

# === STutility（Seurat扩展） ===
library(STutility)
library(Seurat)

# 从Visium读取（基于Seurat）
infoTable <- data.frame(
  samples = "spaceranger_output/filtered_feature_bc_matrix.h5",
  spotfiles = "spaceranger_output/spatial/tissue_positions_list.csv",
  imgs = "spaceranger_output/spatial/tissue_hires_image.png",
  json = "spaceranger_output/spatial/scalefactors_json.json"
)

se <- InputFromTable(infotable = infoTable)  # 创建Seurat对象+空间信息

# Seurat标准流程
se <- se %>%
  SCTransform() %>%                        # SCT标准化
  RunPCA() %>%
  FindNeighbors(dims=1:30) %>%
  FindClusters(resolution=0.5) %>%
  RunUMAP(dims=1:30)

# ★空间可视化（STutility优势：H&E叠加）
ST.FeaturePlot(se,
  features = c("CD3D", "CD68"),            # 基因表达叠加到H&E上
  ncol = 2
)

ST.DimPlot(se,
  group.by = "seurat_clusters",            # 聚类叠加到H&E
  ncol = 1
)

# === 空间分析速查 ===
# Squidpy (Python)
import squidpy as sq
sq.gr.spatial_neighbors(adata)             # 空间邻域
sq.gr.spatial_autocorr(adata, mode="moran")  # 空间自相关
sq.gr.nhood_enrichment(adata, cluster_key="leiden")  # 邻域富集

# Giotto (R)
createSpatialNetwork(obj, method="kNN", k=6)
detectSpatialCorGenes(obj, method="network")

# STutility (R, 基于Seurat)
ST.FeaturePlot(se, features="gene")        # H&E叠加

# 选择: Python→Squidpy | R+简单→STutility | 全面→Giotto