382_单细胞时间序列分析¶

library(Seurat)      # 单细胞分析
library(slingshot)   # 拟时序分析
library(SingleCellExperiment)  # SCE格式
library(ggplot2)

# 将 Seurat 对象转换为 SingleCellExperiment 格式
sce <- as.SingleCellExperiment(seurat_obj)

# 提取 UMAP 坐标作为降维结果
reducedDims(sce)$UMAP <- Embeddings(seurat_obj, "umap")

# 运行 Slingshot（指定起始细胞类型）
sce <- slingshot(
  sce,
  clusterLabels = "cell_type",     # 细胞类型标签列
  reducedDim = "UMAP",             # 使用 UMAP 坐标拟合轨迹
  start.clus = "Stem_cell",        # 起始细胞类型（祖细胞/干细胞）
  end.clus = c("Mature_cell_A", "Mature_cell_B")  # 终点（可选）
)

# 提取拟时序值
pseudotime_matrix <- slingPseudotime(sce)   # 每条轨迹一列
print(head(pseudotime_matrix))

# 取每个细胞的最小拟时序（如果有多条轨迹，取主轨迹时序）
seurat_obj$Pseudotime <- apply(pseudotime_matrix, 1, function(x) {
  x[!is.na(x)][which.min(abs(x[!is.na(x)]))]   # 取最近的非NA值
})

# 可视化：UMAP + 拟时序着色
FeaturePlot(
  seurat_obj,
  features = "Pseudotime",
  reduction = "umap",
  order = TRUE
) + scale_color_viridis_c(option = "B") +    # 黄-橙-紫配色
  ggtitle("Slingshot 拟时序分析")

# 分析沿拟时序的基因表达变化
# 找随拟时序显著变化的基因（用 tradeSeq 包）
library(tradeSeq)

# 过滤掉NA（未处于轨迹上的细胞）
keep_cells <- !is.na(seurat_obj$Pseudotime)
seurat_sub <- seurat_obj[, keep_cells]

# 拟合 GAM（广义加法模型）
sce_sub <- as.SingleCellExperiment(seurat_sub)
sce_sub <- fitGAM(
  counts = assay(sce_sub, "counts"),  # 使用raw counts
  pseudotime = data.frame(
    Lineage1 = seurat_sub$Pseudotime
  ),
  nknots = 6    # 样条曲线节点数
)

# 检验哪些基因随时间显著变化
at_test <- associationTest(sce_sub)    # 关联检验
sig_genes <- at_test[at_test$pvalue < 0.05, ]
sig_genes <- sig_genes[order(sig_genes$waldStat, decreasing = TRUE), ]
cat("随拟时序显著变化的基因数:", nrow(sig_genes), "\n")

import scvelo as scv    # RNA速度分析（pip install scvelo）
import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt

# 设置绘图参数
scv.settings.verbosity = 2    # 显示关键信息
scv.settings.presenter_view = True

# ============================================================
# 前置：需要 STARsolo/Alevin 输出的 spliced + unspliced counts
# ============================================================
# 读取 velocyto 输出的 loom 文件
loom_file = "sample.loom"   # Alevin-fry 或 STARsolo 输出
adata = scv.read_loom(loom_file)

# 或者读取预计算好的 h5ad（已包含 spliced/unspliced 层）
adata = sc.read_h5ad("rna_velocity_data.h5ad")

# 检查必要的层是否存在
print(adata.layers.keys())   # 应该包含 'spliced' 和 'unspliced'

# 步骤1：预处理（过滤、归一化）
scv.pp.filter_and_normalize(
    adata,
    min_shared_counts=20,    # 基因在至少20个细胞中有共享reads
    n_top_genes=2000         # 保留2000个高变基因
)

# 步骤2：PCA 和 UMAP（如果还没做）
sc.pp.pca(adata, n_comps=50)
sc.pp.neighbors(adata, n_pcs=30)
sc.tl.umap(adata)

# 步骤3：计算 RNA 速度
scv.pp.moments(
    adata,
    n_pcs=30,   # 使用前30个PC计算kNN
    n_neighbors=30   # kNN邻居数
)

# 动态模型（推荐，比随机模型更准）
scv.tl.recover_dynamics(adata, n_jobs=4)   # 4线程，计算较慢
scv.tl.velocity(adata, mode="dynamical")   # 动态模型

# 构建速度图
scv.tl.velocity_graph(adata)

# 步骤4：可视化 RNA 速度
scv.pl.velocity_embedding_stream(
    adata,
    basis="umap",               # 在 UMAP 上展示速度流
    color="cell_type",          # 按细胞类型着色
    save="velocity_stream.pdf", # 保存图片
    title="RNA Velocity"
)

# 查看速度置信度（confidence低的细胞轨迹预测不可信）
scv.tl.velocity_confidence(adata)
scv.pl.scatter(
    adata,
    c="velocity_confidence",   # 用置信度着色
    cmap="Greens",
    perc=[5, 95]               # 显示5%-95%范围（去极端值）
)

# 步骤5：细胞命运预测（latent time）
scv.tl.latent_time(adata)
scv.pl.scatter(
    adata,
    color="latent_time",       # 潜在时间（类似拟时序）
    color_map="gnuplot",
    size=80,
    title="Latent Time"
)