hello，大家好，今天我们来分享一个新的分析点，跨样本的基因表达动态。文章在A statistical framework for differential pseudotime analysis with multiple single-cell RNA-seq samples，很有意义的一个分析点。

看看简介

基于单细胞 RNA-seq (scRNA-seq) 数据的伪时间分析已被广泛用于研究沿连续生物过程的动态基因调控程序，如细胞分化、免疫反应和疾病发展。 现有的伪时间分析方法主要解决重建细胞伪时间轨迹和推断一个生物样本中沿重建轨迹的基因表达变化的问题。 随着越来越多地对多个患者样本进行 scRNA-seq 研究，比较不同样本的基因表达动态已成为一种新的需求，但缺乏系统的分析解决方案。 （个人翻译，总的来说，其实我们更加关注相对于normal，疾病样本在轨迹分化的过程中的基因动态）。

其实软件Lamian就是对多样本轨迹分析的补充，考虑到单细胞数据样本包含多种生物学条件（例如，年龄、性别、样本类型、疾病状态等），这个分析软件框架有以下功能。

• (1) 构建细胞伪时间轨迹，评估轨迹分支结构的不确定性。
• (2) 评估与样本生物学因素相关的分支结构的变化。
• (3) 确定沿伪时间的细胞丰度和基因表达的变化。
• (4) 表征样本生物学差异如何修改细胞丰度和基因表达的伪时间动态。

重要的是，在识别与伪时间和样本生物学差异相关的细胞丰度或基因表达变化时，Lamian 考虑了其他现有伪时间方法未考虑的生物样本的变异性。 因此，在分析多样本数据时，Lamian 能够更恰当地控制错误发现率 (FDR)，这是任何其他方法都无法提供的特性。 （这个可能确实需要关注）。

接下来我们详细看看分析过程(使用示例数据)

options(warn=-1)
suppressMessages(library(Lamian))

Module 1: tree variability

Lamian 的模块 1 设计用于检测伪时间树形结构中分支的稳定性。 自动随机选取伪时间树结构中的分支，然后测试它们的检测率（这算是检验分支的稳定性）。 在每次cell-bootstrapping后，重建树结构并重新识别所有分支。 同时应用 Jaccard 统计量和重叠系数作为评估引导设置中的任何分支是否与原始分支之一匹配的统计量。 分支的检测率定义为分支发现它匹配的bootstrap settings的百分比。 模块 1 还报告和测试每个分支中的样本比例。 （有点费劲）。

数据准备

需要一个基因的细胞表达矩阵、细胞的低维表示和细胞注释（每个细胞属于哪个样本）。这里我们就用示例数据分析看看

data(man_tree_data)

man_tree_data 是一个包含基因按细胞表达矩阵、细胞的低维表示（PCA）和细胞注释的列表，其中第一列是细胞名称，第二列是样本名称，行名称是细胞 名称。

推断树形结构

set.seed(12345)
res = infer_tree_structure(pca = man_tree_data[['pca']], 
                           expression = man_tree_data[['expression']], 
                           cellanno = man_tree_data[['cellanno']], 
                           origin.marker = c('CD34'), 
                           number.cluster = 5,
                           xlab='Principal component 1', 
                           ylab = 'Principal component 2')

这个软件看来是按照pca的成分来构建树形结构，我想我们还是可以替换成tSNE或者UMAP，还要指定起始位点的基因特征，分支数量，看来这个软件，需要的前期知识有点多，这里的cluster我们就对应我们聚类的结果或者细胞注释的结果。

绘图

plotmclust(res, cell_point_size = 0.5, 
           x.lab = 'Principal component 1', 
           y.lab = 'Principal component 2')

图片.png

这里我们还是真的需要UMAP的降维结果最好。

估计分支的不稳定性

We can call the function evaluate_uncertainty() to evaluate the tree topology uncertainty. We suggested that users set n.permute = 100 or more to ensure enough randomness to construct the null distribution, but here we set n.permute = 5 in order to provide a simplified example.

result <- evaluate_uncertainty(res, n.permute=5)

这里作为简单的示例，n.permute 设置为5，在实际应用中设置 n.permute = 100 或更多时，结果会更有意义。

Module 2: Evaluate differential topoloty

Lamian 的模块 2 首先识别样本之间树拓扑的变化，然后评估是否存在与样本生物学差异相关的差异拓扑变化。 对于每个样本，Lamian 计算每个树分支中的cell比例，称为分支cell比例。 由于零或低比例可以反映分支的缺失或耗尽，因此可以使用分支cell比例变化来描述树拓扑的变化。 对于多个样本，Lamian 通过估计跨样本的分支cell比例的方差来表征每个分支的跨样本变异。 此外，可以拟合回归模型来测试分支单元格比例是否与样本生物学差异相关。 这允许识别不同条件之间的树拓扑变化，例如在病例相对于normal。

data = result[[2]]
rownames(data) <- c('HSC->lymphocyte', 'HSC->myeloid', 'HSC->erythroid')
design = data.frame(sample= paste0('BM',1,8), sex = c(0, rep(1,4), rep(0,3)))
branchPropTest(data = data, design = design)

Module 3: Trajectory differential tests about gene expression

XDE test(这个大家自行查看，就不多做解释了)。

Res <- lamian.test(expr = expdata$expr, 
                   cellanno = expdata$cellanno, 
                   pseudotime = expdata$pseudotime, 
                   design = expdata$design, 
                   test.type = 'variable', 
                   permuiter = 5)

downstream analysis 1: visualize and cluster XDE genes based on their multiple-sample temporal patterns across sample covariates

stat <- Res$statistics
stat <- stat[order(stat[,1], -stat[,3]), ]
diffgene <- rownames(stat[stat[, grep('^fdr.*overall$', colnames(stat))] < 0.05,])
Res$populationFit <- getPopulationFit(Res, gene = diffgene, type = 'variable')
## clustering
Res$covariateGroupDiff <- getCovariateGroupDiff(testobj = Res, gene = diffgene)
Res$cluster <- clusterGene(Res, gene = diffgene, type = 'variable', k=5)
colnames(Res$populationFit[[1]]) <- colnames(Res$populationFit[[2]]) <- colnames(Res$expr) 
plotXDEHm(Res, cellWidthTotal = 180, cellHeightTotal = 350, subsampleCell = FALSE, sep = ':.*')

图片.png

图片还是不错的。

We can plot the cluster mean and difference.

## plotClusterMeanAndDiff
plotClusterMeanAndDiff(Res)

We can also plot the cluster difference seperately by calling plotClusterDiff(testobj=Res, gene = diffgene).

downstream analysis 2: Gene ontology (GO) enrichment analysis

We can further check out whether the DDG in each cluster have any enriched genome functions by applying Gene Ontology (GO) analysis library(topGO); goRes <- GOEnrich(testobj = Res, type = 'variable'). Please note that there are no enriched GO terms in any of the clusters in this simplified example since we only subset 1,000 genes. In practice, if there are significant GO terms in the clusters, we can generate a heatmap of the top GO terms using the command plotGOEnrich(goRes = goRes).

TDE test(这个大家也自行查看)

Res <- lamian.test(expr = expdata$expr, 
                   cellanno = expdata$cellanno, 
                   pseudotime = expdata$pseudotime, 
                   design = expdata$design, 
                   test.type = 'time', 
                   permuiter = 5)
diffgene <- rownames(Res$statistics)[Res$statistics[,1] < 0.05]
Res$populationFit <- getPopulationFit(Res, gene = diffgene, type = 'time')
Res$cluster <- clusterGene(Res, gene = diffgene, type = 'time', k=3)
plotTDEHm(
  Res,
  subsampleCell  = FALSE,
  showCluster = TRUE,
  type = 'time',
  cellWidthTotal = 200,
  cellHeightTotal = 200
)

图片.png

plotClusterMean(testobj=Res, cluster = Res$cluster, type = 'time')

图片.png

downstream analysis 2: GO analysis for TDE genes

Similar to XDE test, we can further proceed to identify the enriched GO terms for each gene cluster using commands library(topGO); goRes <- GOEnrich(testobj = Res, type = 'time', sep = ':.*'). Please note that there are no enriched GO terms in any of the clusters in this simplified example since we only subset 1,000 genes. In practice, if there are enriched GO terms, we can plot them using commands plotGOEnrich(goRes = goRes).

Module 4: trajectory differential test based on cell composition(这部分跟上面的类似，就不多介绍了)

总之，软件有其特点，可以借鉴，取其长处，为己所用。

生活很好，有你更好