1简介

lineagpulse是一种用于单细胞RNA-seq (scRNA-seq)数据的差分表达算法。lineagpulse基于零膨胀负二项式噪声模型，可以捕获离散和连续的群体结构:离散群体结构是细胞组(例如实验条件或tSNE集群)。例如，连续群体结构可以是细胞的伪时间顺序或细胞的时间顺序。lineagpulse的主要用途和新颖性在于它能够很好地拟合细胞假时间顺序上的基因表达轨迹。请注意，lineagpulse并不推断假时间顺序，而是一个下游分析工具，用于分析给定假时间顺序上的基因表达轨迹(例如来自扩散假时间或monocle2)。

为了在scRNA-seq数据上运行LineagPulse，用户需要为包装器函数runlineagpulse使用最小的输入参数集，然后执行所有归一化、模型拟合和差分表达式分析步骤，而不需要任何更多的用户交互:

• 计数矩阵(基因x细胞)绝不能以任何方式归一化。一个有效的输入选项是预期从校准器计数。注意TPM或深度规范化预期计数不是计数数据。lineagpulse的统计框架基于matCounts包含计数数据的假设。计数矩阵也可以作为.mtx文件(稀疏矩阵格式文件)的路径提供，或者作为summarizeexperimental或singlecel实验对象提供。
• dfAnnotation数据帧，包含单元格注释。dfAnnotation的行名必须等于matCounts的列名，因为两者都对应单元格。注意，如果counts是summarizeexperiment或singlecel实验对象，如果dfAnnotation为NULL，注释数据帧将被取为colData(counts)。dfAnnotation必须包含
• 如果拟合一个连续模型，则一个名为“continuous”的列(例如，strMuModel是脉冲或样条，表达式是时间或伪时间坐标的函数拟合)
• 如果一个离散种群模型是拟合的(例如strMuModel是组，表达式是组赋值的函数，例如集群或实验条件)，则一个名为“groups”的列，其中包含将细胞赋值到这些组的字符串
• 描述批处理结构的列(如果有的话)。

此外，还可以提供:

• matPiConstPredictors退出率的基因特异性预测因子矩阵(基因x预测因子)。我们建议使用log平均表达(除非strDropModel为“logistic_ofMu”)和可能影响测序效率的潜在参数，如基因的GC含量。
• strMuModel用于差分表达式分析的替代模型的表达式模型类型:“脉冲”表示脉冲模型，“样条”表示自然三次样条模型。
• vecConfoundersMu一个字符串向量，对应于dfAnnotation中的列名，该列名描述了要更正的批处理结构。
• scaDFSplinesMu如果strMuModel设置为“splines”，则平均参数的基于样条模型的自由度。
• vecnormconst要使用的外部单元格归一化常数(例如测序深度校正因子)。元素的名称必须是matCounts(单元格)的列名。
• 设置并行化的进程数。
• boolVerbose在包装器函数运行时输出基本进度报告。
• boolSuperVerbose为包装器函数的每个步骤输出详细的进度报告。

最后，对lineagpulse的数学和算法基础有扎实掌握的有经验的用户可以更改这些高级输入选项的默认值:

• vecConfoundersDisp批处理变量用于纠正离散度(方差)。
• strdismodelfull用于完整模型的分散模型。
• strdismodelred用于简化模型的色散模型。
• strDropModel要使用的退出模型。
• strDropFitGroup接收一个退出模型参数化的单元格组。
• scaDFSplinesDisp离散参数的基于样条模型的自由度，如果strdismodel被设置为“splines”。
• boolEstimateNoiseBasedOnH0是否在空表达式模型或替代表达式模型上估计退出模型。注意，将其设置为FALSE会大大增加运行时间。
• scaMaxEstimationCycles退出和表达式模型估计迭代周期的最大数目。

2微分表达式分析

在这里，我们提出了一个纵向排序的微分表达式分析方案。差分表达式结果在一个数据帧中，可以通过properties($)这样的列表从输出对象访问该数据帧。差异表达分析的核心结果是差异表达的p值和假发现率校正p值，这是对非恒定表达模型(脉冲、样条或组)与恒定表达模型进行基因假设检验的结果。

library(LineagePulse) lsSimulatedData <- simulateContinuousDataSet(scaNCells = 100, scaNConst = 10, scaNLin = 10, scaNImp = 10, scaMumax = 100, scaSDMuAmplitude = 3, vecNormConstExternal=NULL, vecdisexternal =rep(20,30)， vecgenewisedropoutates =rep(0.1, 30))

画平均轨迹

##统一设置大小因子=1

##绘制色散

模拟负二项噪声

##模拟辍学

objLP <- runlineagpulse (counts = lsSimulatedData$counts, dfAnnotation = lsSimulatedData$annot)

计数数据的lineagpulse: v1.18.0

##——数据预处理

100个细胞中有0个没有连续的协变量，因此被排除在外。

30个基因中有0个不包含非零观察值，因此被排除在分析之外。

## # 100个单元格中有0个不包含非零观测值，因此被排除在分析之外。

计算归一化常量:

## #所有大小因子都设置为1。

H1和H0都适合ZINB模型。

## ### a)拟合ZINB模型a (H0: mu=常数disp=常数)与噪声模型。

## #。初始化:ll -26109.5891848408

## # 1。用ll -14357.6637411404迭代0.03分钟。

## # 2。在0.03分钟内迭代ll -14187.8787319392。

## # 3。在0.02分钟内迭代ll -14187.8782829273。

在0.11 min内完成了零膨胀负二项式模型A与噪声模型的拟合。

## ### b)拟合ZINB模型b (H1: mu=样条disp=常数)。

## #。初始化:ll -15775.1929206345

## # 1。在0.02分钟内迭代ll -13802.8519233628。

0.03分钟拟合零膨胀负二项式模型B

## ### c)拟合NB模型A (H0: mu=常数disp=常数)。

## #。初始化:ll -15467.7936881531

## # 1。在0.01分钟内迭代ll -15329.4995983373。

0.02分钟完成NB模型B的拟合。

## ### d)拟合NB模型B (H1: mu=样条disp=常数)。

## #。初始化:ll -15799.3286387389

## # 1。用ll -15255.2272440065迭代0.02分钟。

0.03分钟完成NB型B的拟合

ZINB装配时间:0.24分钟

运行差异表达式分析。

##完成运行lineagpulse()。

头(objLP dfResults美元)

## gene_1 gene_1 4.712952e-01 6.147329e-01 78.52002 4.712952e-01 0.994102 ## gene_3 gene_3 2.084452e-01 3.291240e-01 36.02987 2.084452e-01 0.994102 ## gene_5 gene_5 4.836890e -02 4.238977e-02 48.79000 1.836890e-02 0.994102 ## gene_6 gene_6 4.433910e-07 1.477970e-06 26.51645 4.433910e-07 0.994102 ## df_full_zinbdf_red_zinb df_full_nb df_red_nb loglik_full_zinb # # gene_1 7 2 7 2 -467.7836 # # gene_2 7 2 7 2 -451.6175 # # gene_3 7 2 7 2 -401.1858 # # gene_4 7 2 7 2 -471.6015 # # gene_5 7 2 7 2 -435.9502 # # gene_6 7 2 7 2 -365.0421 # # loglik_red_zinb loglik_full_nb loglik_red_nb allZero # #假# # gene_2 gene_1 -470.0659 -536.1737 -536.3946 -455.3063 -509.0278 -509.8222假# # gene_3 -404.7697 -450.1640 -452.2496假# # gene_4假# # gene_5 -477.4103 -527.5220 -528.3468 -442.7496 -487.4970 -488.0378错误## gene_6 -383.8676 - 41010.3024 -411.7151

除了原始p值之外，人们可能会对表达式模型的进一步细节感兴趣，例如作为伪时间函数的表达式平均值的形状，对数折叠变化(LFC)和作为伪时间函数的全局表达式趋势。我们将在下面分别讨论这些后续问题。请注意，所有这些后续问题都是基于适合计算差分表达式p值的模型来回答的。因此，一旦runlineagpulse()被调用一次，就不需要进一步的模型拟合。

#进一步检查结果#绘制基因轨迹

有多个选项可用于基因表达轨迹绘制:观察结果可以通过dropout的后验概率(boolcolorbydropout)进行着色。观测值可以基于替代表达式模型进行归一化，也可以作为散点图的原始观测值(boolH1NormCounts)。谱系轮廓可以添加到伪时间相关效应(boolLineageContour)中，以辅助非均匀人口密度的可视化解释。如果折叠变化很大，则可以显示日志计数而不是计数(boolLogPlot)。在任何情况下，基因表达轨迹分析器绘图函数plotGene的输出对象是一个ggplot2对象，然后可以对其进行打印或修改。

#绘制差异表达p值最低的基因gplotExprProfile <- plotGene(objLP = objLP, boolLogPlot = FALSE, strGeneID = objLP$dfResults[it .min(objLP$dfResults$p)，]$gene, boolLineageContour = FALSE) gplotExprProfile

函数plotGene还显示了负二项式噪声模型(“H1(NB)”)下的H1模型拟合，作为参考，以显示如果不考虑退出，模型拟合是什么样的。

#提取平均参数拟合p值最低的基因的每个细胞。matMeanParamFit <- getFitsMean(lsMuModel = lsMuModelH1(objLP)， vecGeneIDs = objLP$dfResults[it .min(objLP$dfResults$p)，]$gene) cat("最小拟合平均参数:"，round(min(matMeanParamFit)，1)))

最小拟合平均参数:66.7

cat("均值拟合均值参数:"，round(均值(matMeanParamFit)，1))

平均拟合平均参数:180.2

#首先，提取适合给定基因的模型vecMeanParamFit <- getFitsMean(lsMuModel = lsMuModelH1(objLP)， vecGeneIDs = objLP$dfResults[it .min(objLP$dfResults$p)，]$gene) #计算轨迹上从第一个到最后一个细胞的log2倍变化idxFirstCell <- where .min(dfAnnotationProc(objLP)$pseudotime) idxLastCell <- where .max(dfAnnotationProc(objLP)$pseudotime) cat("LFC轨迹上第一个到最后一个细胞:， round((log(vecMeanParamFit[idxLastCell]) - log(vecMeanParamFit[idxFirstCell])) / log(2)，1))

LFC轨迹上的第一个到最后一个单元格:

#计算表达式轨迹cat("LFC最小值到模型拟合的最大表达式值:"，round((log(max(vecMeanParamFit)) - log(min(vecMeanParamFit))) / log(2)，1))

## LFC最小到最大模型拟合表达式值:2.2

#创建所有差异表达基因的热图lsHeatmaps <- sortGeneTrajectories(vecIDs = objLP$dfResults[which(objLP$dfResults$padj < 0.01)，]$gene, lsMuModel = lsMuModelH1(objLP)， dirHeatmap=NULL) print(lsHeatmaps$hmGeneSorted)

3.会话信息

sessionInfo ()

## R版本4.2.1(2022-06-23)##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 20.04.5 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas。/home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack。所以## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_GB LC_COLLATE= c# # [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats graphics grDevices utils datasets methods基础## ##其他附加包:## [1]lineagpulse_1 .18.0 BiocStyle_2.26.0 ## ##通过命名空间加载(且未附加):# # # # [1] MatrixGenerics_1.10.0 Biobase_2.58.0 [3] sass_0.4.2 splines_4.2.1 # # [5] jsonlite_1.8.3 foreach_1.5.2 # # [7] gtools_3.9.3 bslib_0.4.0 # # [9] assertthat_0.2.1 highr_0.9 # # [11] BiocManager_1.30.19 stats4_4.2.1 # # [13] GenomeInfoDbData_1.2.9 yaml_2.3.6 # # [15] pillar_1.8.1 lattice_0.20-45 # # [17] glue_1.6.2 digest_0.6.30 # # [19] GenomicRanges_1.50.0 RColorBrewer_1.1-3 # # [21] XVector_0.38.0 colorspace_2.0-3 # # [23] htmltools_0.5.3 Matrix_1.5-1 # # [25] pkgconfig_2.0.3 GetoptLong_1.0.5 # #[27] magick_1 .7.3 bookdown_0.29 ## [29] zlibbioc_1.44.0 scales_1.2.1 ## [31] BiocParallel_1.32.0 tibble_3.1.8 ## [33] farver_2.1.1 generics_0.1.3 ## [35] IRanges_2.32.0 ggplot2_3.3.6 ## [37] cachem_1.0.6 SummarizedExperiment_1.28.0 ## [39] BiocGenerics_0.44.0 cli_3.4.1 ## [41] magrittr_2.0.3 crayon_1.5.2 ## [43] evaluate_0.17 fansi_1.0.3 ## [45] gplots_3.1.3 doParallel_1.0.17 ## [47] Cairo_1.6-0 tools_4.2.1 ## [49] GlobalOptions_0.1.2 lifecycle_1.0.3 ## [51] matrixStats_0.62.0ComplexHeatmap_2.14.0 # # [53] stringr_1.4.1 S4Vectors_0.36.0 # # [55] munsell_0.5.0 cluster_2.1.4 # # [57] DelayedArray_0.24.0 compiler_4.2.1 # # [59] jquerylib_0.1.4 GenomeInfoDb_1.34.0 # # [61] caTools_1.18.2 rlang_1.0.6 # # [63] grid_4.2.1 rcurl_1.98 - 1.9 # # [65] iterators_1.0.14 rjson_0.2.21 # # [67] SingleCellExperiment_1.20.0 circlize_0.4.15 # # [69] labeling_0.4.2 bitops_1.0-7 # # [71] rmarkdown_2.17 gtable_0.3.1 # # [73] codetools_0.2-18 DBI_1.1.3 # # [75] R6_2.5.1 knitr_1.40 # # [77] dplyr_1.0.10fastmap_1.1.0 ## [79] utf8_1.2.2 clue_0.3-62 ## [81] KernSmooth_2.23-20 shape_1.4.6 ## [83] stringi_1.7.8 Rcpp_1.0.9 ## [85] parallel_4.2.1 vctrs_0.5.0 ## [87] png_0.1-7 tidyselect_1.2.0 ## [89] xfun_0.34