掌握单细胞动态：手把手教你RNA速率分析（附代码）

# 1. 加载必要包 --------------------------------------------------------------
suppressPackageStartupMessages({
  library(Seurat)       # 单细胞数据分析基石
  library(SeuratDisk)   # 用于处理 h5ad 和 loom 文件
  library(SeuratWrappers) # Seurat 的辅助包，提供与其他工具的接口
  library(velocyto.R)   # RNA 速率分析的核心算法包
  library(ggplot2)      # 强大的绘图工具
  library(Matrix)       # 处理稀疏矩阵
  library(loomR)        # 直接读取 loom 文件
  library(patchwork)    
})

# 2. 读取loom文件并创建Seurat对象 --------------------------------------------
process_loom <- function(loom_path) {
# 尝试使用SeuratDisk直接读取
if (requireNamespace("SeuratDisk", quietly = TRUE)) {
    tryCatch({
      cat("Attempting to read loom file with SeuratDisk...\n")
      seurat_obj <- LoadLoom(loom_path,
                             assay = "RNA",
                             layers = c("spliced", "unspliced"),
                             verbose = FALSE)
      cat("Successfully loaded loom with SeuratDisk\n")
      # 确保 "counts" 层设置为 "spliced" 数据
      if ("spliced" %in% Layers(seurat_obj[["RNA"]])) {
        LayerData(seurat_obj, assay = "RNA", layer = "counts") <- LayerData(seurat_obj, assay = "RNA", layer = "spliced")
      }
      # 检查 "unspliced" 层是否存在
      if ("unspliced" %in% Layers(seurat_obj[["RNA"]])) {
        cat("Unspliced layer loaded by LoadLoom\n")
      } else {
        stop("Unspliced layer not found after LoadLoom")
      }
      return(seurat_obj)
    }, error = function(e) {
      cat("SeuratDisk loading failed, using manual method...\n")
    })
  }
# 手动读取方法（当SeuratDisk无法加载时备用）
  cat("Reading loom file manually...\n")
  loom <- connect(loom_path, mode = "r", skip.validate = TRUE)
  spliced <- t(loom[["layers/spliced"]][, ])
  unspliced <- t(loom[["layers/unspliced"]][, ])
  cell_ids <- loom[["col_attrs/CellID"]][]
  gene_names <- loom[["row_attrs/Gene"]][]
  colnames(spliced) <- colnames(unspliced) <- cell_ids
  rownames(spliced) <- rownames(unspliced) <- gene_names
  loom$close_all()
  cat("Creating Seurat object...\n")
  seurat_obj <- CreateSeuratObject(
    counts = spliced,
    assay = "RNA",
    project = "RNA_velocity",
    min.cells = 3,
    min.features = 200
  )
# 添加 "unspliced" 作为 "RNA" assay 的一个层
  LayerData(seurat_obj, assay = "RNA", layer = "unspliced") <- unspliced
return(seurat_obj)
}

# 3. 自定义函数：绘制剪接与未剪接计数散点图 ---------------------------------
plot_spliced_unspliced <- function(seurat_obj, gene) {
# 从 "RNA" assay 中提取 spliced 和 unspliced 数据
  spliced <- GetAssayData(seurat_obj, assay = "RNA", layer = "counts")[gene, ]
  unspliced <- GetAssayData(seurat_obj, assay = "RNA", layer = "unspliced")[gene, ]
# 创建数据框
  df <- data.frame(spliced = spliced, unspliced = unspliced)
# 使用 ggplot2 创建散点图，应用 log 变换以处理稀疏数据
  p <- ggplot(df, aes(x = spliced + 1, y = unspliced + 1)) +
    geom_point(alpha = 0.5) +
    scale_x_log10() +
    scale_y_log10() +
    labs(title = gene, x = "Spliced + 1", y = "Unspliced + 1") +
    theme_minimal()
return(p)
}

# 4. 数据预处理 --------------------------------------------------------------
preprocess_data <- function(seurat_obj) {
  cat("Preprocessing data...\n")
  DefaultAssay(seurat_obj) <- "RNA"
# 保存 unspliced 数据以防被覆盖
if ("unspliced" %in% Layers(seurat_obj[["RNA"]])) {
    unspliced_data <- LayerData(seurat_obj, assay = "RNA", layer = "unspliced")
  } else {
    stop("Unspliced layer missing before preprocessing")
  }
# 基础过滤
  seurat_obj <- subset(seurat_obj, subset = nFeature_RNA > 200 & nCount_RNA < 50000)
# 标准化和特征选择
  seurat_obj <- NormalizeData(seurat_obj)
  seurat_obj <- FindVariableFeatures(seurat_obj, selection.method = "vst", nfeatures = 3000)
  seurat_obj <- ScaleData(seurat_obj, features = rownames(seurat_obj))
# 降维与聚类
  seurat_obj <- RunPCA(seurat_obj, npcs = 50, verbose = FALSE)
  seurat_obj <- FindNeighbors(seurat_obj, dims = 1:30)
  seurat_obj <- FindClusters(seurat_obj, resolution = 0.8)
  seurat_obj <- RunUMAP(seurat_obj, dims = 1:30, n.neighbors = 30, min.dist = 0.3)
# 恢复 unspliced 层（重要！确保 unspliced 数据在处理后依然存在）
  LayerData(seurat_obj, assay = "RNA", layer = "unspliced") <- unspliced_data
  cat("Unspliced layer restored after preprocessing\n")
return(seurat_obj)
}

# 5. RNA速率分析 --------------------------------------------------------------
run_velocity <- function(seurat_obj) {
  cat("Running RNA velocity analysis...\n")
# 打印可用层名以便调试
  cat("Layers in RNA assay:", Layers(seurat_obj[["RNA"]]), "\n")
# 从 "RNA" assay 中获取 spliced 和 unspliced 数据
  emat <- as.matrix(GetAssayData(seurat_obj, assay = "RNA", layer = "counts"))      # 已剪接数据
  nmat <- as.matrix(GetAssayData(seurat_obj, assay = "RNA", layer = "unspliced"))  # 未剪接数据
# UMAP 坐标（嵌入），用于后续可视化
  umap_coords <- Embeddings(seurat_obj, "umap")
# 计算 PCA 坐标上的欧几里得距离矩阵，用于构建细胞间的邻近关系
  pca_coords <- Embeddings(seurat_obj, "pca")[, 1:30]  # 使用前30个主成分
  cell_dist <- dist(pca_coords, method = "euclidean")
# 执行速度分析
  rvel <- gene.relative.velocity.estimates(
    emat = emat,
    nmat = nmat,
    deltaT = 1,          # 时间步长，通常设为1
    kCells = 25,         # 用于局部估计的邻近细胞数
    fit.quantile = 0.02, # 拟合模型时使用的分位数
    cell.dist = cell_dist, # 细胞距离矩阵
    verbose = TRUE
  )
# 将结果写入 Seurat 对象的 @misc 槽，方便后续调用
  seurat_obj@misc[["RunVelocity"]] <- rvel
return(seurat_obj)
}

# 6. 可视化与保存结果 --------------------------------------------------------
visualize_results <- function(seurat_obj, output_dir = "results") {
  cat("Visualizing results...\n")
if (!dir.exists(output_dir)) dir.create(output_dir)
# 保存 Seurat 对象，包含所有分析结果
  saveRDS(seurat_obj, file.path(output_dir, "velocity_seurat.rds"))
# 1. UMAP 聚类图：展示细胞群体分布
  pdf(file.path(output_dir, "umap_clusters.pdf"), width = 8, height = 6)
print(DimPlot(seurat_obj, reduction = "umap", label = TRUE) +
          ggtitle("Cell Clusters"))
  dev.off()
# 2. RNA 速率矢量图：细胞命运的动态轨迹
  pdf(file.path(output_dir, "rna_velocity.pdf"), width = 10, height = 8)
  ident.colors <- scales::hue_pal()(length(levels(seurat_obj)))
  names(ident.colors) <- levels(seurat_obj)
  cell.colors <- ident.colors[Idents(seurat_obj)]
  names(cell.colors) <- colnames(seurat_obj)
  show.velocity.on.embedding.cor(
    emb = Embeddings(seurat_obj, "umap"),
    vel = seurat_obj@misc[["RunVelocity"]],
    n = 300,             # 在图中显示的细胞子集数量（为避免过于拥挤）
    scale = "sqrt",      # 矢量长度的缩放方式
    cell.colors = ac(cell.colors, alpha = 0.5), # 细胞颜色及透明度
    cex = 0.8,           # 细胞点的大小
    arrow.scale = 3,     # 箭头长度的缩放因子
    show.grid.flow = TRUE, # 是否显示网格流线
    grid.n = 50,         # 网格密度
    main = "RNA Velocity"
  )
  dev.off()
# 3. 剪接/未剪接基因图：展示部分高变基因的剪接状态
  pdf(file.path(output_dir, "spliced_vs_unspliced.pdf"), width = 9, height = 7)
  var_genes <- VariableFeatures(seurat_obj)
  genes_to_plot <- if (length(var_genes) > 20) sample(var_genes, 9) else var_genes
# 使用自定义函数生成并排列散点图
  plots <- lapply(genes_to_plot, function(gene) plot_spliced_unspliced(seurat_obj, gene))
  p <- wrap_plots(plots, ncol = 3) # 使用 patchwork 组合多图
print(p)
  dev.off()
  cat("Results saved to:", output_dir, "\n")
}