PseudotimeDE在单细胞伪时间中的建模分析
前言
PseudotimeDE的安装
该R包尚未在Bioconductor 或 CRAN 上公布，因此在安装的时候需要利用devtools
install.packages("devtools")
library(devtools)
# 由于我缺少两个依赖包，故需要先安装
BiocManager::install('slingshot')
BiocManager::install('BiocStyle')
devtools::install_github("SONGDONGYUAN1994/PseudotimeDE",dependencies = T)
广义加性模型（GAM）简介
当响应变量和决策变量之间的关系不明确时，即不是能用线性回归或者多项式回归明确关系的时候，往往采用的是广义加性模型来解决此类问题，
对比线性回归，多项式回归和一般加性模型（从上到下依次是线性回归，多项式回归和广义加性模型）：
对于广义加性模型来说，f(x)为未指明的函数，需要进行非参数估计，广义加性模型允许我们在预先不知晓因变量与自变量之间关系的情况下，使用非线性平滑项来拟合模型。而非线性平滑项的做法是，当在全局角度无法知晓响应变量和决策变量之间的关系时，通常将自变量划分未多个区间，每一个区间都用单独的线性函数或非线性的低阶多项式函数来拟合。这种方法被称为样条函数（Spline），也成为样条插值法，属于数值分析里面的内容。其生成的回归线为平稳、光滑的曲线，因此经样条函数变化后的自变量也被称作非线性平滑项
那么R包mgcv中的gam()函数可以完成广义加性模型的建立：
# 仅有非线性平滑项，三次样条
fit <- gam(y ~ s(x, k = 3, bs=’cr’), data = data)
## 其中s(x, k = 3, bs=’cr’)表示的是非线性平滑项，cr为三次样条
## k表示平滑程度，越小越平滑，反之越摇摆
# 线性项和非线性平滑项，三次样条
fit <- gam(y ~ x1 + s(x2, k = 3, bs=’cr’), data = data)
## 其中x1表示线性项
## 其中s(x2, k = 3, bs=’cr’)内表示的是非线性平滑项，cr为三次样条
## k表示平滑程度，越小越平滑，反之越摇摆
# 仅有非线性平滑项，自然样条
fit <- gam(y ~ ns(x, k = 3, df=6), data = data)
## 其中ns(x, k = 3, df=6)表示的是非线性平滑项，ns()代表自由样条
## k表示平滑程度，越小越平滑，反之越摇摆
关于广义加性模型的学习资料目前国内还没有比较好的，主要可以参考：
- https://cran.r-project.org/web/packages/gam/gam.pdf
- https://www.andrew.cmu.edu/user/achoulde/95791/lectures/lecture02/lecture02-95791.pdf
- https://www.r-bloggers.com/2017/07/generalized-addictive-models-and-mixed-effects-in-agriculture/
- https://m-clark.github.io/generalized-additive-models/building_gam.html
- https://s3.amazonaws.com/assets.datacamp.com/production/course_6413/slides/chapter2.pdf
- http://environmentalcomputing.net/intro-to-gams/
关于R包mgcv中的gam()函数的结果解读：
其中：
- parametric coefficient：Intercept代表截距项的参数估计和显著性p_val
- Approximate significance of smooth terms：(1).edf 为 1 等价于线性关系；(2).edf > 1 且 edf ≤ 2 是弱非线性关系；(3). edf > 2 表示高度非线性关系
- R-sq 代表R方的值
- Deviance explained 代表方差解释率
依据广义加性模型的原理，上图表示为
expv = 9.5944935 - 1.0480144 * s(pseudotime).1 + 0.2672279 * s(pseudotime).2 - 1.7394072 * s(pseudotime).3 + 2.2212207 * s(pseudotime).4 + 1.3692433 * s(pseudotime).5
而 s(pseudotime).1， s(pseudotime).2，s(pseudotime).3， s(pseudotime).4， s(pseudotime).5代表关于pseudotime不同的样条插值函数，这里一共有5个样条插值函数是因为我们选取 k = 6，6个节点分出5个区段，每一个区段对应一个拟合程度最好的样条插值函数s(pseudotime)（s(pseudotime).1， s(pseudotime).2，s(pseudotime).3， s(pseudotime).4， s(pseudotime).5），然后5个样条插值函数乘上对应的系数相加，即为最终的加性模型方程
以三次样条插值为例，s(pseudotime).2 代表在区间2中，s(pseudotime).2这个三次多项式拟合效果最好
如何理解样条插值函数：
我们举一个例子：
x <- seq(0, 8, 0.1)
sin_x <- sin(x)
y <- sin_x + rnorm(n = length(x), mean = 0, sd = sd(sin_x / 2))
Sample_data <- data.frame(y,x)
gam_y <- gam(y ~ s(x,k = 6, bs = "cr"))
model_matrix <- predict(gam_y, type = "lpmatrix")
plot(y ~ x)
abline(h = 0)
## 画出s(x).1和s(x).2
lines(x, model_matrix[, "s(x).1"], type = "l", lty = 2)
lines(x, model_matrix[, "s(x).2"], type = "l", lty = 2)
表一
model_matrix 表示的是每一个样条插值函数的值是以x为自变量预测的y值，那么上图每一列（2~6列）对应的数值是样条函数s()的预测值（y值）
下图的纵坐标即为上表（表一）每个s()函数（2-6列）对应的值，每一个值对应图中的一个点
上图表示的是两个样条函数s(x).1和样条函数s(x).2的曲线图
我们画出所有的样条函数：
plot(y ~ x)
abline(h = 0)
x_new <- seq(0, max(x), length.out = 100)
y_pred <- predict(gam_y, data.frame(x = x_new))
## 所有的样条函数s()的曲线
matplot(x, model_matrix[,-1], type = "l", lty = 2, add = T)
## 最终的拟合曲线
lines(y_pred ~ x_new, col = "red", lwd = 2)
粗红色曲线代表最终拟合的曲线，该例子的加性模型方程为：
y = 0.15944433 + 0.86656365 * s(x).1 - 0.27461751 * s(x).2 - 1.14086262 * s(x).3 + 0.04573517 * s(x).4 + 1.01969190 * s(x).5
也就是每一个非线性多项式s(x)乘对应系数，然后累加起来即为最后的拟合曲线
相比于普通的回归模型拟合，广义加性模型的优点是拟合程度更高，更精确
PseudotimeDE代码分析
数据预处理准备：
data(LPS_ori_tbl, package = "PseudotimeDE")
data(LPS_sub_tbl, package = "PseudotimeDE")
# 以CCL5为例
gene = "CCL5"
ori.tbl = LPS_ori_tbl
sub.tbl = LPS_sub_tbl[1:50]
sce = LPS_sce
model = "nb"
k = 6
knots = c(0:5/5)
fix.weight = TRUE
aicdiff = 10
seed = 123
LPS_ori_tbl
LPS_ori_tbl表示的是每一个细胞对应的pseudotime
我们对PseudotimeDE的核心函数进行改写如下：
set.seed(seed)
## Avoid R checking warning
splits <- time.res <- NULL
# 建立基因表达量和pseudotime之间的非线性平滑项，方法是三次样条
fit.formula <- stats::as.formula(paste0("expv ~ s(pseudotime, k = ", k, ",bs = 'cr')"))
# 获取细胞数
num_cell <- length(ori.tbl$cell)
# 获取细胞对应的pseudotime
pseudotime <- ori.tbl$pseudotime
# 获取总细胞数
num_total_cell <- dim(sce)[2]
# 获取每个细胞基因"CCL5"的表达量（count值）
expv <- assays(sce)$counts[gene, ori.tbl$cell]
count.v <- expv
# 整合pseudotime和expv
dat <- cbind(pseudotime, expv) %>% as_tibble()
# 对每个细胞基因"CCL5"的表达量（count值）进行标准化
expv.quantile <- quantile(log(count.v + 1))
# 对每个细胞基因"CCL5"的表达量（count值）计算均值
expv.mean <- mean(log(count.v + 1))
# 对每个细胞基因"CCL5"的表达量（count值）统计0表达情况
expv.zero <- sum(count.v == 0)/length(count.v)
其中dat表示每一个细胞对应相应的pseudotime和基因表达量（基因CCL5）
dat
如果以 model = "nb" 为例：
# 数据输入准备
dat = dat
zinf = FALSE
use_weights = FALSE
k = 6
knots = c(0:5/5)
nthreads = 1
cellWeights <- rep(1, dim(dat)[1])
fit.formula <- stats::as.formula(paste0("expv ~ s(pseudotime, k = ", k, ",bs = 'cr')"))
## fit.formula = expv ~ s(pseudotime, k = 6, bs = "cr")
fit <- gam(fit.formula, family = nb(link = "log"),
data = dat,
knots = list(pseudotime = knots), control = list(nthreads = nthreads))
## 其中knots的用处是控制节点
# 我们可以利用model.matrix()函数来查看模型对应的响应变量值
Xp <- model.matrix(fit)
怎么理解广义加性模型在pseudotime中的应用
上述曲线即为每个细胞中CCL5表达量与pseudotime的一个曲线关系，而如何描述这种曲线关系是否是显著的，我们主要看下图的指标
该例子中，如果截距和非线性平滑项是显著的，那么我们就认为模型是显著的，因此就认为该基因的表达量和pseudotime之间满足如下函数关系式：
- expv = 9.5944935 - 1.0480144 * s(pseudotime).1 + 0.2672279 * s(pseudotime).2 - 1.7394072 * s(pseudotime).3 + 2.2212207 * s(pseudotime).4 + 1.3692433 * s(pseudotime).5
- 其中s(pseudotime)表示的是非线性平滑项，可以理解为一个区间内拟合得最好的多项式
由图片上来看，如果拟合曲线越抖动，那么该基因表达量随pseudotime的变化就越明显，应用于基因表达量和时间序列的分析中可以解释为该基因的表达量随时间波动明显
p值矫正方式
那么PseudotimeDE除了去发现一些表达量随pseudotime变化相关的基因，但是往往gam模型显著的事实上并不显著，而在前言中，作者就阐明了PseudotimeDE中一种更为严格p值矫正方法
PseudotimeDE不仅保证产生的p值有严格的统计学意义，还实现了更高的检验效力（power）和更好的对错误发现率的控制（FDR control）。总而言之，PseudotimeDE 对利用单细胞测序数据解析细胞动态有着重要意义。
那么经过最严格的p值筛选后：
data(LPS_ori_tbl, package = "PseudotimeDE")
data(LPS_sub_tbl, package = "PseudotimeDE")
res <- PseudotimeDE::runPseudotimeDE(gene.vec = "CXCL10",
ori.tbl = LPS_ori_tbl,
sub.tbl = LPS_sub_tbl[1:50],
sce = LPS_sce,
model = "nb",
mc.cores = 1)
那么挑选与gam模型显著的基因重点看para.pv，以该值作为显著性去做筛选即可
