单细胞轨迹分析（Trajectory Analysis）¶

# Seurat对象转Monocle3
library(Seurat)
library(monocle3)
library(SeuratWrappers)

# 假设已有seurat_obj
# 子集选取相关细胞
sub_obj <- subset(seurat_obj, idents = c("HSC", "Progenitor", "Monocyte", "DC"))

# 转换为monocle3 cds对象
cds <- as.cell_data_set(sub_obj)
cds <- cluster_cells(cds, resolution = 1e-3)

# 学习轨迹图
cds <- learn_graph(cds, use_partition = TRUE)

# 可视化轨迹
plot_cells(cds,
           color_cells_by = "cluster",
           label_groups_by_cluster = FALSE,
           label_leaves = FALSE,
           label_branch_points = FALSE)

# 指定root并排序（交互式或程序化）
# 方法1：交互式
cds <- order_cells(cds)

# 方法2：程序化指定root
# 找到特定cluster中最早分化状态的细胞
get_earliest_principal_node <- function(cds, cluster_name = "HSC"){
  cell_ids <- which(colData(cds)[, "seurat_clusters"] == cluster_name)
  closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
  closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
  root_pr_nodes <- igraph::V(principal_graph(cds)[["UMAP"]])$name
  root_pr_nodes <- closest_vertex[cell_ids, ]
  root_pr_nodes <- names(which.max(table(root_pr_nodes)))
  root_pr_nodes
}

cds <- order_cells(cds, root_pr_nodes = get_earliest_principal_node(cds))

# 可视化pseudotime
plot_cells(cds,
           color_cells_by = "pseudotime",
           label_cell_groups = FALSE,
           label_leaves = FALSE,
           label_branch_points = FALSE,
           graph_label_size = 1.5)

# 找沿轨迹变化基因
track_genes <- graph_test(cds, neighbor_graph = "principal_graph", cores = 4)
track_genes <- track_genes[order(track_genes$morans_I, decreasing = TRUE), ]
head(track_genes, 20)

library(slingshot)
library(SingleCellExperiment)

# 从Seurat转SCE
sce <- as.SingleCellExperiment(sub_obj)

# 运行slingshot
sce <- slingshot(sce,
                 clusterLabels = 'seurat_clusters',
                 reducedDim = 'PCA',
                 start.clus = "0")  # 指定起始cluster

# 获取pseudotime
pseudotime_mat <- slingPseudotime(sce)
head(pseudotime_mat)

# 获取曲线
curves <- slingCurves(sce)

# 可视化
library(ggplot2)
library(RColorBrewer)

# 在UMAP上画轨迹
umap_coords <- reducedDim(sce, "UMAP")
df <- data.frame(UMAP1 = umap_coords[,1],
                 UMAP2 = umap_coords[,2],
                 pseudotime = pseudotime_mat[,1],
                 cluster = colData(sce)$seurat_clusters)

ggplot(df, aes(x = UMAP1, y = UMAP2, color = pseudotime)) +
  geom_point(size = 0.5) +
  scale_color_viridis_c() +
  theme_minimal()

# 找沿轨迹差异表达基因（tradeSeq）
library(tradeSeq)
# 拟合GAM模型
sce_gam <- fitGAM(counts = counts(sce),
                   sds = SlingshotDataSet(sce),
                   nknots = 6,
                   verbose = TRUE)

# 关联检验
assoc_res <- associationTest(sce_gam)
assoc_res <- assoc_res[order(assoc_res$pvalue), ]
head(assoc_res, 20)

# 1. 用velocyto从BAM生成loom文件
velocyto run10x /path/to/cellranger_output /path/to/genes.gtf

# 或使用STARsolo直接输出spliced/unspliced
STAR --runMode alignReads \
     --soloType CB_UMI_Simple \
     --soloFeatures Gene Velocyto \
     --outSAMtype BAM SortedByCoordinate \
     ...

# 2. Python: scVelo分析
import scvelo as scv
import scanpy as sc
import anndata as ad

# 读取数据
adata = sc.read_h5ad("processed_adata.h5ad")
ldata = scv.read("velocyto_output.loom", cache=True)

# 合并spliced/unspliced信息
adata = scv.utils.merge(adata, ldata)

# 预处理
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# 3. 动态模型（推荐）
scv.tl.recover_dynamics(adata, n_jobs=8)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)

# 4. 可视化
scv.pl.velocity_embedding_stream(adata, basis='umap', color='celltype')
scv.pl.velocity_embedding_grid(adata, basis='umap', color='celltype')

# 5. Latent time（类似pseudotime）
scv.tl.latent_time(adata)
scv.pl.scatter(adata, color='latent_time', cmap='gnuplot')

# 6. 找driver基因
top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index[:20]
scv.pl.scatter(adata, basis=top_genes[:5], ncols=5, frameon=False)

# CytoTRACE（R版本）
# 安装
# devtools::install_github("digitalcytometry/cytotrace")
library(CytoTRACE)

# 输入：基因×细胞的表达矩阵（raw counts或normalized均可）
expr_matrix <- as.matrix(GetAssayData(sub_obj, slot = "counts"))

# 运行CytoTRACE
results <- CytoTRACE(expr_matrix, ncores = 8)

# 结果包含：
# CytoTRACE score: 每个细胞的分化潜能（0=分化完全, 1=最干）
# CytoTRACErank: 排名

# 添加到Seurat对象
sub_obj$CytoTRACE <- results$CytoTRACE[colnames(sub_obj)]

# 可视化
FeaturePlot(sub_obj, features = "CytoTRACE") +
  scale_color_viridis_c(option = "plasma")

# 按cluster查看分化潜能
VlnPlot(sub_obj, features = "CytoTRACE", group.by = "seurat_clusters")

# CytoTRACE 2（Python版本）
import cytotrace2 as ct2

# 输入AnnData
results_adata = ct2.cytotrace2(adata)

# 结果存储在 adata.obs['CytoTRACE2_Score']
sc.pl.umap(results_adata, color='CytoTRACE2_Score', cmap='RdYlBu_r')

# Monocle3: 模块化基因表达
# 已有track_genes结果
sig_genes <- track_genes %>%
  filter(q_value < 0.05) %>%
  arrange(desc(morans_I)) %>%
  head(200)

# 基因模块分析
gene_module_df <- find_gene_modules(cds[rownames(cds) %in% sig_genes$gene_short_name, ],
                                     resolution = 1e-2)

# 热图展示基因模块沿轨迹的表达
plot_cells(cds,
           genes = gene_module_df %>% filter(module %in% c(1, 2, 3)),
           label_cell_groups = FALSE,
           show_trajectory_graph = FALSE)

# tradeSeq: 沿轨迹的基因表达趋势
library(tradeSeq)
plotSmoothers(sce_gam, counts(sce), gene = "CD34")
plotSmoothers(sce_gam, counts(sce), gene = "GATA1")

# 在分支点比较两条lineage
branch_res <- earlyDETest(sce_gam, knots = c(3, 4))
branch_res <- branch_res[order(branch_res$pvalue), ]

# ===== 环境安装 =====
# R包安装
# install.packages("BiocManager")
# BiocManager::install(c("monocle3", "slingshot", "tradeSeq"))
# devtools::install_github("digitalcytometry/cytotrace")

# Python包安装
pip install scvelo scanpy velocyto anndata

# velocyto安装（需要samtools）
pip install velocyto
conda install -c bioconda samtools

# ===== velocyto生成loom =====
velocyto run10x \
  /data/cellranger/sample1/outs \
  /ref/refdata-gex-GRCh38-2020-A/genes/genes.gtf \
  -o /data/velocity/

# ===== 完整Monocle3流程 =====
# 见上文R代码

# ===== scVelo完整流程 =====
# 见上文Python代码