单细胞实战之Roe、Augur、miloR——入门到进阶(中级篇4）

单细胞实战之Ro/e、Augur、miloR——入门到进阶(中级篇4）

该推文首发于公众号：单细胞天地

在上一讲中完成了inferCNV的分析流程，本讲将回顾学习单细胞分析中的三个工具：Ro/e、Augur、miloR，这三个工具分别能够判断细胞簇在组织分布倾向、扰动影响和细胞差异丰度。

本次内容涉及到的工程文件可通过网盘获得：中级篇2，链接: 提取码: yx93 。

此外，可以向“生信技能树”公众号发送关键词‘单细胞’，直接获取Seurat V5版本的完整代码。

Ro/e——判断细胞簇在组织分布倾向的分析工具

Ro/e分析又称为STARTRAC-distribution(张泽民院士团队开发的STARTRAC框架中的其中一个工具)的主要目的是通过计算组织中观测细胞数与期望细胞数（期望细胞数通过卡方检验计算得出）的比率。

与卡方检验中定义的 (观测值-期望细胞数)²/期望细胞数仅能反映观测值与期望细胞数的偏离程度不同，STARTRAC通过Ro/e 可进一步指示特定细胞簇在特定组织中是否富集或耗竭。如下图中是分析不同的T细胞簇/巨噬细胞簇在不同组织中的Ro/e富集情况。

Ro/e 指数的评分规则：

如果 Ro/e>1，则表明在该组织中，某细胞簇的细胞数量高于期望细胞数，即表现为富集；如果 Ro/e<1，则表明某细胞簇在该组织中的细胞数量低于期望细胞数，即表现为消耗。 +++, Ro/e > 1; ++, 0.8 < Ro/e ≤ 1; +, 0.2 ≤ Ro/e ≤ 0.8; +/−, 0 < Ro/e < 0.2; −, Ro/e = 0

个人认为期望细胞数也是基于现有数据进行测算的并且最后生成的值是比率，所以其实最终结果跟细胞数本身很有关系，同时细胞数量也不应太少，容易出现很大的偏倚。

Ro/e 指数分析步骤

1.导入

代码语言：javascript代码运行次数：0运行复制

rm(list = ls())
library(Startrac)
library(ggplot2)
library(tictoc)
library(ggpubr)
library(ComplexHeatmap)
library(RColorBrewer)
library(circlize)
library(tidyverse)
library(sscVis)
library(Seurat)
library(tidyverse)
library(readr)
library(qs)
library(BiocParallel)
register(MulticoreParam(workers = 8, progressbar = TRUE)) 

scRNA <- qread("./9-CD4+T/CD4+t_final.qs")
Idents(scRNA) <- "celltype"
DimPlot(scRNA)

2.数据预处理

代码语言：javascript代码运行次数：0运行复制

data  <- scRNA@meta.data
colnames(data)
#  [1] "orig.ident"         "nCount_RNA"         "nFeature_RNA"       "percent_mito"      
#  [5] "percent_ribo"       "percent_hb"         "RNA_snn_res.0.01"   "seurat_clusters"   
#  [9] "RNA_snn_res.0.05"   "RNA_snn_res.0.1"    "RNA_snn_res.0.2"    "RNA_snn_res.0.3"   
# [13] "RNA_snn_res.0.5"    "RNA_snn_res.0.8"    "RNA_snn_res.1"      "celltype"          
# [17] "location"           "pANN_0.25"          "DF.classifications" "contamination"     
# [21] "RNA_snn_res.0.4"    "RNA_snn_res.0.6"    "RNA_snn_res.0.7"    "RNA_snn_res.0.9"   
# [25] "RNA_snn_res.1.1"    "RNA_snn_res.1.2"    "RNA_snn_res.1.3"    "RNA_snn_res.1.4"   
# [29] "RNA_snn_res.1.5"    "RNA_snn_res.1.6"    "RNA_snn_res.1.7"    "RNA_snn_res.1.8"   
# [33] "RNA_snn_res.1.9"    "RNA_snn_res.2" 

# 需要得到clone.id(ID标识),patient(样本),majorCluster(细胞亚群),loc(不同组织/位置)列
# 并且需要重新命名
data <- data[,c(34,1,16,17)]
colnames(data) <- c("clone.id","patient","majorCluster","loc")
data$Cell_Name <- c("T cell")
data$patient <- as.character(data$patient)
data$clone.id <- as.character(data$clone.id)
data$majorCluster <- as.character(data$majorCluster)
data$loc <- as.character(data$loc)

3.Ro/e分析

代码语言：javascript代码运行次数：0运行复制

Roe <- calTissueDist(data,
         byPatient = F,
         colname.cluster = "majorCluster", # 不同细胞亚群
         colname.patient = "patient", # 不同样本
         colname.tissue = "loc", # 不同组织
         method = "chisq", # "chisq", "fisher", and "freq" 
         min.rowSum = 0) 
Roe
  #               left     right
  # ELK4+T    1.3286764 0.7309727
  # Naive T   0.6669732 1.2725882
  # Th1       0.8664169 1.1093401
  # Th17      1.2137669 0.8250281
  # Tm        0.9822346 1.0145413
  # Treg      1.3638576 0.7021762
  # ZNF793+T  1.0462073 0.9621785
  # ZSCAN12+T 0.5643888 1.3565553
  # 这是最关键的值，然后进行可视化即可

4.可视化

代码语言：javascript代码运行次数：0运行复制

table(data$majorCluster)
# ELK4+T   Naive T       Th1      Th17        Tm      Treg  ZNF793+T ZSCAN12+T 
#    204      1459       818       421      1045      1528       189       248 
col_fun <- colorRamp2(c(min(Roe, na.rm = TRUE), 1, max(Roe, na.rm = TRUE)), 
                      c("blue", "white", "red"))
Heatmap(as.matrix(Roe),
        show_heatmap_legend = TRUE, 
        cluster_rows = TRUE, 
        cluster_columns = TRUE,
        row_names_side = 'right', 
        show_column_names = TRUE,
        show_row_names = TRUE,
        col = col_fun,
        row_names_gp = gpar(fontsize = 10),
        column_names_gp = gpar(fontsize = 10),
        heatmap_legend_param = list(
          title = "Ro/e Index",  # 自定义图注名称
          at = seq(0.5, 1.5, by = 0.5), # 例刻度的位置/自己的数据必须修改一下！
          labels = seq(0.5, 1.5, by = 0.5) # 每个刻度的标签/自己的数据必须修改一下！
        ),
        cell_fun = function(j, i, x, y, width, height, fill) {
          grid.text(sprintf("%.2f", Roe[i, j]), x, y, gp = gpar(fontsize = 8, col = "black"))
        }
)
# cell_fun = function(j, i, x, y, width, height, fill)
# j 和 i：分别表示单元格的列和行索引。
# x 和 y：是单元格的中心位置坐标，用于确定文本的位置。
# width 和 height：表示单元格的宽度和高度，可以用来调整文本位置或大小。
# fill：表示单元格的背景颜色。

setwd("..")

从细胞数和Ro/e值可以看到左半结肠中存在较多数量的Treg细胞以及较高的Ro/e值，Naive T细胞则在右半结肠中更多和显著的Ro/e值。因此，这也符合现有的生物学知识，右半结肠癌通常具有更高的免疫原性而左半结肠癌则相反。

Augur——判断扰动影响的分析工具

Augur能够在单细胞数据中优先考虑对生物扰动最敏感的细胞类型，其采用机器学习框架来量化扰动和未扰动细胞在高维空间中的可分性(比如不同的疾病状态，药物刺激等)。Augur 使用的机器学习框架依赖于分类算法，如随机森林（Random Forests）或逻辑回归（Logistic Regression），这些算法非常适合处理具有高维特征空间的单细胞数据。

Augur的具体做法是：通过判断是否能够仅凭分子特征准确预测出每个细胞的实验处理标签（如处理组 vs 对照组），来量化这种可分离性。在实际操作中，它会为每种细胞类型单独训练一个机器学习模型，用于预测每个细胞来自哪个实验条件。接着，通过交叉验证评估每个细胞类型特定分类器的准确性，从而为细胞类型的优先级排序提供一个定量的依据。

Augur分析流程

1.导入

代码语言：javascript代码运行次数：0运行复制

rm(list = ls())
library(stringr)
library(Seurat)
library(Augur)
library(dplyr)
library(patchwork)
library(viridis)
library(qs)
library(BiocParallel)
register(MulticoreParam(workers = 8, progressbar = TRUE)) 

scRNA <- qread("./9-CD4+T/CD4+t_final.qs")
Idents(scRNA) <- "celltype"
DimPlot(scRNA)

dir.create("./12-STARTRAC-miloR-Augur")
setwd("./12-STARTRAC-miloR-Augur")

2.单一组别不同条件下augur分析

比如治疗前后,但是本次分析就把左右半结肠位置数据当做实验分组了。

代码语言：javascript代码运行次数：0运行复制

# 此函数的主要目的是评估并计算每种细胞类型在特定条件（如实验组别或时间点）下的表现，通常通过计算区域下曲线（AUC）来量化。它是单一条件下的细胞类型表现的度量，可以帮助识别在特定实验设置下哪些细胞类型表现最为显著或重要。
augur <- calculate_auc(scRNA,
                       cell_type_col = "celltype", # 细胞亚群数据
                       label_col = "location", # 实验分组数据
                       n_threads = 8)
qsave(augur,"augur.qs")

head(augur$AUC, 5)
# # A tibble: 5 × 2
#   cell_type   auc
#   <chr>     <dbl>
# 1 Naive T   0.628
# 2 Th17      0.621
# 3 Th1       0.604
# 4 ZSCAN12+T 0.588
# 5 ELK4+T    0.586

3.可视化

代码语言：javascript代码运行次数：0运行复制

# 明确哪个细胞群体的RNA水平变化最大
# plot_lollipop函数是基于ggplot框架的,所以很容易就可以美化图片
# 还可以提取原始函数进行修改
augur <- qread("augur.qs")
library(ggplot2)
plot_lollipop(augur) +
  geom_segment(aes(xend = cell_type, yend = 0.5), size = 1) +
  geom_point(size = 3, aes(color = cell_type)) +  # 增加点的大小和颜色映射
  scale_color_manual(values = c(
    "Naive T" = "#1f77b4",       # 蓝色
    "Th17" = "#ff7f0e",    # 橙色
    "Th1" = "#2ca02c",    # 绿色
    "ZSCAN12+T" = "#d62728",        # 红色
    "ELK4+T" = "#9467bd",   # 紫色
    "Treg" = "#8c564b",      # 棕色
    "Tm" = "#e377c2",     # 粉色
    "ZNF793+T" = "#7f7f7f"    # 灰色
  )) 

# 结合umap
# rank模式
plot_umap(augur,scRNA,
          mode = "rank",
          reduction = "umap",
          palette = "cividis",
          cell_type_col = "celltype")

# default模式
plot_umap(augur,scRNA,
          mode = "default",
          reduction = "umap",
          palette = "YlGnBu", #  "viridis", "plasma", "magma", "inferno"
          cell_type_col = "celltype")

setwd("..")

plot_lollipop展示了在不同处理情况下(group)中变化最显著的细胞(AUC值)，从结果来看Naive T细胞受到的扰动程度是最大的，ZNF793+T收到的扰动程度则是最小的。

mode中的rank和default模式的区别在于是否用AUC值进行颜色映射。

miloR——判断细胞差异丰度的分析工具

miloR能够随机定义细胞中的细胞节点，然后通过K最近邻法去识别所定义节点与周围的其他细胞节点之间的邻近关系，找到跟这些节点最近的(比如欧几里得法)其他细胞节点，将这些具有邻近关系的细胞小群看做不同的群组，此时可以想象一下，每个群组中都会存在具有疾病和正常信息的许多细胞，接着再将这些群组中疾病与否的细胞数分别进行比较和汇总，最后可以得到基于某一大类细胞中内部细胞变化趋势的结果(比如是否倾向与疾病的方向)。

有了这个结果之后，结合预先定义的细胞分群结果有助于使用者更好的去了解疾病与否、疾病状态、药物刺激等情况下的不同细胞亚群的变化状态，因此使用miloR可以有效的作为单细胞水平上细胞差异分析的补充工具。

miloR分析流程

1.导入

代码语言：javascript代码运行次数：0运行复制

rm(list = ls())
library(stringr)
library(Seurat)
library(miloR)
library(scater)
library(scran)
library(dplyr)
library(patchwork)
library(qs)
library(BiocParallel)
register(MulticoreParam(workers = 8, progressbar = TRUE)) 

scRNA <- qread("./9-CD4+T/CD4+t_final.qs")
Idents(scRNA) <- "celltype"
DimPlot(scRNA)

dir.create("./12-STARTRAC-miloR-Augur")
setwd("./12-STARTRAC-miloR-Augur")

2.数据预处理

需要一个seurat对象，建议提前预处理

代码语言：javascript代码运行次数：0运行复制

# check
table(scRNA$orig.ident)
 #    GSM5688706_WGC     GSM5688707_JCA GSM5688708_LS-CRC3 GSM5688709_RS-CRC1  GSM5688710_R_CRC3 
 #               738                542               1381               1879                807 
 # GSM5688711_R_CRC4 
 #               565
table(scRNA$celltype)
   # ELK4+T   Naive T       Th1      Th17        Tm      Treg  ZNF793+T ZSCAN12+T 
   #    204      1459       818       421      1045      1528       189       248

3.构建miloR对象

代码语言：javascript代码运行次数：0运行复制

# 由于需要输入SingleCellExperiment对象
# 因此先进行转化一下
scRNA_pre <- as.SingleCellExperiment(scRNA)
# 构建miloR对象
scRNA_milo <- Milo(scRNA_pre)
reducedDim(scRNA_milo,"UMAP") <- reducedDim(scRNA_pre,"UMAP")
scRNA_milo
# class: Milo 
# dim: 21308 5912 
# metadata(0):
# assays(3): counts logcounts scaledata
# rownames(21308): RP11-34P13.7 FO538757.2 ... GRIK1-AS1 AP001626.2
# rowData names(0):
# colnames(5912): GSM5688706_WGC_AACACACGTATCACGT-1 GSM5688706_WGC_AACACACTCAAGGTGG-1 ...
#   GSM5688711_R_CRC4_TTTGGTTTCTGGGATT-1 GSM5688711_R_CRC4_TTTGTTGGTAATGCTC-1
# colData names(35): orig.ident nCount_RNA ... RNA_snn_res.2 ident
# reducedDimNames(3): PCA HARMONY UMAP
# mainExpName: RNA
# altExpNames(0):
# nhoods dimensions(2): 1 1
# nhoodCounts dimensions(2): 1 1
# nhoodDistances dimension(1): 0
# graph names(0):
# nhoodIndex names(1): 0
# nhoodExpression dimension(2): 1 1
# nhoodReducedDim names(0):
# nhoodGraph names(0):
# nhoodAdjacency dimension(2): 1 1

4.构建K最邻图谱

代码语言：javascript代码运行次数：0运行复制

# d值与降维时选定的pca值一致
# k值与FindNeighbors时选定的值一致
scRNA_milo <- buildGraph(scRNA_milo, k = 30, d = 15,reduced.dim = "PCA")

5.在 KNN 图上定义具有代表性的邻域

代码语言：javascript代码运行次数：0运行复制

# 官网建议prop在0.1-0.2
# K和d和KNN图参数一致
scRNA_milo <- makeNhoods(scRNA_milo, prop = 0.2, k = 30, d = 15, 
                    refined = TRUE, reduced_dims = "PCA")
plotNhoodSizeHist(scRNA_milo)

开发者建议在50到100达到峰值，否则需要增加k/prop的值。

6.计算neighbourhoods中的细胞

计算每个样本中有多少细胞位于每个邻域中。

代码语言：javascript代码运行次数：0运行复制

scRNA_milo <- countCells(scRNA_milo, 
                    meta.data = as.data.frame(colData(scRNA_milo)),
                    sample="orig.ident") # sample代表样本

# colData(scRNA) 相当于代表了metadata的信息
# sample代表了样本数
head(nhoodCounts(scRNA_milo))
# 这个步骤能够给scRNA数据中增加一个n*m的matrix
# n代表不同的neighborhood数量,m代表实验样本
# 6 x 6 sparse Matrix of class "dgCMatrix"
#   GSM5688706_WGC GSM5688707_JCA GSM5688708_LS-CRC3 GSM5688709_RS-CRC1 GSM5688710_R_CRC3 GSM5688711_R_CRC4
# 1             19              7                 41                 18                 .                 .
# 2              1              .                  5                  3                 2                61
# 3              5              2                 84                  3                 .                 .
# 4             10             21                  4                 30                 .                 .
# 5             21             14                 59                 15                 .                 .
# 6             71              .                 19                  4                 .                 3

这为对象添加了一个矩阵，其中n是邻域的数量，

科技改变生活-雨落星辰 - 所有的伟大,都源于一个勇敢的开始

单细胞实战之Roe、Augur、miloR——入门到进阶(中级篇4）