前言

读《Cell type prioritization in single-cell data》，并学习利用AUC来分析相同细胞类型在两个处理组别的差异

目前在单细胞转录组中比较不同生物学状态之间差异的方式主要由两种：
（1）基于统计学假设检验方法比较实验组和对照组之间不同细胞类型比例是否由明显的差异
（2）基于统计学模型比较不同分组中分子水平的差异

因此，作者提出了一个问题，即哪一种cell type对单细胞数据中不同的处理（treat/control）最为敏感

原理

作者采用机器学习的方式，对某一cell type通过其在两个处理组别之间建立分类模型，然后训练该模型，比较预测的结果与实验的标签，通过AUC值来反应该cell type在这两个组别中的差异程度

比方说，对于某一种cell type，在WT/KO两个组中，会有不同的“表现”。因此在二分类模型中，AUC的值大小反应了该cell type在WT/KO两个组别之间的差异程度

AUC是表征分类模型性能的指标，AUC越高，分类性能越好，说明该cell type在两组中表达矩阵“特征”可以有效区分出WT/KO两个组，因此该cell type在两个组中差异明显（这里的差异指的是表达差异）
注：虽然同属于该cell type，但是在WT/KO组中基因的表达还是有一定差异的

input：

output：

AUC值表征该cell type在两个组别中差异的程度

代码

作者将某一cell type的两个组别的label作为响应变量（分类变量）;将该cell type在两个组别中的单细胞表达矩阵作为决策变量

# y 为响应变量（分类变量），cell_types[1]代表其中一种cell type
y = labels[cell_types == cell_types[1]]
# X 为决策变量；expr代表这一种cell type在两个组别中的单细胞基因表达矩阵
## cell_types[1]代表其中一种cell type
X = expr[, cell_types == cell_types[1]]

核心代码如下：
建立模型

# set up model
## 利用随机森林建模
if (classifier == "rf") {
  importance = T
  clf = rand_forest(trees = !!rf_params$trees,
                    mtry = !!rf_params$mtry,
                    min_n = !!rf_params$min_n,
                    mode = mode) %>%
    set_engine(rf_engine, seed = 1, importance = T, localImp = T)
} 
## 利用线性模型建模
else if (classifier == "lr") {
  family = ifelse(multiclass, 'multinomial', 'binomial')
  if (is.null(lr_params$penalty) || lr_params$penalty == 'auto') {
    # first, get optimized penalty for dataset
    lr_params$penalty = withCallingHandlers({
      glmnet::cv.glmnet(X0 %>%
                          ungroup() %>%
                          select(-label) %>%
                          as.matrix() %>%
                          extract(, setdiff(colnames(.), 'label')),
                        X0$label,
                        nfolds = folds,
                        family = family) %>%
        extract2("lambda.1se")
    }, warning = function(w) {
      if (grepl("dangerous ground", conditionMessage(w)))
        invokeRestart("muffleWarning")
    })
  }
  
  clf = logistic_reg(mixture = lr_params$mixture,
                     penalty = lr_params$penalty,
                     mode = 'classification') %>%
    set_engine('glmnet', family = family)
}

建模以后用测试集去评估分类的性能，计算AUC值，来评估差异