把差异分析结果拆分到不同单细胞亚群

把差异分析结果拆分到不同单细胞亚群 by 生信技能树

差异分析大家都没有问题了，无论是表达量芯片还是转录组测序后的counts矩阵，都有各自的R包，很容易就拿到上下调的基因列表。如果是GEO数据库里面的表达量芯片数据处理，主要的难点是表达量矩阵获取和探针的基因名字转换，搞定后只需要一定的生物学背景对数据进行合理的分组后就是标准的差异分析，富集分析。主要是参考我八年前的笔记：

其他形式的矩阵也是可以差异分析的，比如蛋白质组学后的表达量矩阵，以及代谢组，等等。又或者是甲基化芯片后的矩阵，甲基化测序的 WGBS和RRBS，还有芯片是最高频的甲基化技术，其中甲基化芯片数据处理我是有视频课程的，首先需要阅读我在生信技能树的甲基化系列教程，目录如下：

但是因为传统的表达量矩阵都是bulk时代的产物，而且单细胞费用一直是居高不下，所以即使到2023大家的首选仍然是bulk测序，普通的bulk测序的转录组才四五百块钱，但是单细胞仍然是接近两万块钱一个样品。

我们针对bulk时代的传统的表达量矩阵进行分组后的差异分析，拿到的上下调基因这个时候有两个不同的属性：

不同分组本身不同单细胞亚群比例有变化，所以其单细胞亚群特异性基因会差异
不同分组的多个单细胞亚群都有同样的通路上下调导致bulk水平的统计学显著的上下调基因

这样的两个情况就需要区分一下，如果我们怀疑是不同分组本身不同单细胞亚群比例有变化，其实可以做一个去卷积的单细胞亚群比例推断，就可以验证。如果我们仅仅是想看自己的bulk水平的统计学显著的上下调基因是否有单细胞亚群特异性，一个R包：TOAST (TOols for the Analysis of heterogeneouS Tissues) 推荐给大家：

使用起来非常方便，就一个csTest函数即可，示范代码是：

Design_out <- makeDesign(design, Prop)
fitted_model <- fitModel(Design_out, Y_raw)
fitted_model$all_coefs # list all phenotype names
fitted_model$all_cell_types # list all cell type names
# coef should be one of above listed phenotypes
# cell_type should be one of above listed cell types
res_table <- csTest(fitted_model, coef = "age", 
                    cell_type = "Neuron", contrast_matrix = NULL)
head(res_table)

可以看到，运行这个csTest函数的流程需要起码3个元素：

Y _ raw 的表达量芯片矩阵或 DNA 甲基化信号值矩阵
一个称为 Design 的样本信息数据框
以及一个称为 prop 的细胞组成表(即混合比例)

Y _ raw 也可以是 SummarizedLaboratory 对象格式，而不是数据矩阵格式，不过这个不重要。我们的bulk时代的表达量芯片或者甲基化芯片很容易拿到矩阵，然后矩阵里面的每个样品都有自己的属性，比如疾病或者对象，处理和对照，性别，这样的信息，就是称为 Design 的样本信息数据框。比较麻烦的是一个称为 prop 的细胞组成表(即混合比例)，如果是甲基化芯片矩阵，可以使用EpiDISH-根据甲基化信号值推断样品的细胞成分，如果是表达量芯片矩阵，可以使用使用MuSiC来根据单细胞亚群特征基因集对bulk转录组数据推断细胞成分。

示例数据

加载这个包，就可以得到3种测试数据，两个甲基化芯片，一个普通表达量芯片矩阵，都有各自的样品分组以及各自的细胞比例矩阵。所以很容易运行这个csTest函数的流程。

## ----loadData-----------------------------------------------------------------
library(TOAST)
data("RA_100samples")

# 第一个例子数据集是450K DNA 甲基化数据。基于 GSE42861提供的原始数据，
# 我们获取并处理了该数据集。这是类风湿关节炎患者和对照组的 DNA 甲基化450K 数据。
# 原始数据集有485577个特征和689个样本。
# 我们已经将随机选择的50名 RA 患者和50名对照者的数据集减少到3000个 CpG。

Y_raw <- RA_100samples$Y_raw
Pheno <- RA_100samples$Pheno
Blood_ref <- RA_100samples$Blood_ref

## ----checkData----------------------------------------------------------------
dim(Y_raw) 
Y_raw[1:4,1:4]

## ----checkPheno---------------------------------------------------------------
dim(Pheno)
head(Pheno, 3)
table(Pheno$disease)

我们这里选择了第一个是甲基化示例数据，所以使用EpiDISH-根据甲基化信号值推断样品的细胞成分，代码如下所示：

refinx <- findRefinx(Y_raw, nmarker = 1000) 
Y <- Y_raw[refinx,]
Ref <- as.matrix(Blood_ref[refinx,])
library(EpiDISH)
outT <- epidish(beta.m = Y, ref.m = Ref, method = "RPC")
estProp_RB <- outT$estF
estProp_RB[1:4,1:4]

library(TOAST) 
K=6
set.seed(1234)
outRF1 <- csDeconv(Y_raw, K, TotalIter = 30, bound_negative = TRUE)  
estProp_RF_improved <- assignCellType(input=outRF1$estProp,
                                      reference=estProp_RB) 
mean(diag(cor(estProp_RF_improved, estProp_RB)))
Prop <- estProp_RF_improved
colnames(Prop) <- colnames(Ref) 
head(Prop)

如果是表达量芯片矩阵，需要使用MuSiC来根据单细胞亚群特征基因集对bulk转录组数据推断细胞成分。

 "CD8T"  "CD4T"  "NK"    "Bcell" "Mono"  "Gran"

这个时候，甲基化信号值矩阵把血液bulk矩阵拆分成为了常见的淋巴细胞以及髓系免疫细胞的比例。

有了甲基化信号矩阵，样品的各个单细胞亚群比例矩阵，以及样品分组，就很容易运行这个csTest函数的流程。


Design_out <- makeDesign(design, Prop) 
fitted_model <- fitModel(Design_out, Y_raw)
# print all the cell type names
fitted_model$all_cell_types
# print all phenotypes
fitted_model$all_coefs
## -----------------------------------------------------------------------------
res_table <- csTest(fitted_model, 
                    coef = "disease", 
                    cell_type = "Gran")
head(res_table, 3)
Disease_Gran_res <- res_table

其中 coef 参数可以控制我们想做的表型差异情况，示例数据里面的是类风湿关节炎患者和对照组这样的疾病状态。然后 cell_type 参数里面控制我们想看的单细胞亚群的差异情况。最后得到res_table就是差异分析结果啦，而且是针对粒细胞这个免疫细胞亚群的。

这样的话，我们就把常规的差异分析结果给细化到了单细胞亚群的水平，虽然我们并没有做指定的单细胞层面的数据。

学徒作业

使用包里面的第二个示例数据，GSE65133里面的芯片表达量矩阵，是PBMC，所以很容易使用PBMC3K的单细胞数据来去卷积，看 BCells ， CD8T ， CD4T ， NK cells ，Monocytes 这些单细胞亚群的比例的推断结果。就使用MuSiC来根据单细胞亚群特征基因集对bulk转录组数据推断细胞成分。

ixxmu / mp_duty