内容

标志

1简介

在这里，我们描述了一个用于推断单细胞RNA-seq数据集的细胞周期位置的包。其中包括理论论证和基准(Zheng et al. 2022)．在我们手中，我们的方法(称为TriCycle)在各种数据模式中稳健地工作，包括跨物种(人类和小鼠)，细胞类型和检测技术(10X, Fluidigm C1);我们还没有遇到这种方法不起作用的数据集。主要输出是单元周期内相对时间的连续估计，表示为0到2pi之间的数字(我们称之为单元周期位置)。除了估计过程之外，我们还包括了一些方便的函数来可视化细胞周期时间，我们还提供了一个离散细胞周期阶段预测器的实现。

2先决条件

库(三轮车)

我们建议用户从a开始SingleCellExperiment对象。输出通常是SingleCellExperiment添加了新的信息。函数作用于矩阵或SummarizedExperiment对象，尽管输出会发生变化，因为这些类型的对象不具备存储输入对象和估算值的能力。

在软件包中，我们包含了一个示例SingleCellExperiment数据集，它是2个样本的小鼠神经球RNAseq的真实子集。从样本AX1和AX2中随机抽取200个细胞。此数据集与我们用于构建单元周期嵌入的数据相同。

数据(neurosphere_example)

重要的:请注意用户应标准化库大小之前放入三轮车的功能。库大小归一化可以通过normalizeCounts函数嘘包装或通过计算CPM值。

3.包功能概述

该方法基于获取一个新的数据集并将其投影到表示单元周期的嵌入中。这个嵌入是使用参考数据集构造的。也许令人惊讶的是，我们发现同样的嵌入适用于我们所研究的所有实验，包括跨物种、细胞类型和数据集。我们将这个嵌入空间作为包的一部分提供，并且我们不期望用户更改这个嵌入(尽管包中的函数支持其他嵌入)。

方法很简单:取一个新的数据集，将其投射到潜在空间中，并推断细胞周期时间。这里的关键功能是

下一步是验证成功预测了细胞周期时间。这包括观察一些有用的图表。这涉及到许多有用的可视化。例如，请注意我们使用的配色方案-因为单元格周期时间“环绕”\([0, 2π\]\)在间隔中，使用一个“环绕”的调色板是非常有用的。相关功能包括

我们还提供了一个单独的细胞周期阶段预测器，预测5个不同的阶段;estimate_Schwabe_stage ()．这个预测器是对作者提出的方法的一个小修改(Schwabe et al. 2020)．我们包含这个函数只是为了方便。
这不是三轮车方法的一部分!

最后，我们有一组函数来创建你自己的参考潜空间。

4将单细胞数据集投影到预学习的细胞周期空间

project_cycle_space ()将自动工程与名称化验logcounts在没有任何其他参数输入的情况下嵌入到单元格周期中。您可以指定物种(默认为鼠标)、基因ID、基因ID类型和AnnotationDb如果需要基因映射，则反对。指人(project_cycle_space)获取详细信息。

neurosphere_example <- project_cycle_space(neurosphere_example) #>不提供自定义参考投影矩阵。从鼠标Neuroshpere数据中学习到的ref将被使用。新数据中发现的投影基因数量为500。neurosphere_example #>类:singlecel实验#> dim: 1500 400 #>元数据(0):#> assays(2): counts logcounts #> rownames(1500): ENSMUSG00000056763 ENSMUSG00000025925…#> ENSMUSG00000044221 ENSMUSG00000040234 #>行数据名称(2):基因接入#> colnames(400): AX1_CTACGTCCAAACCCAT AX1_AACTCTTTCATTCACT…#> AX2_CGAGAAGCACTACAGT AX2_CATATTCGTCTGCGGT #> colData names(2): TotalUMIs sample #> reducedDimNames(1): tricycleEmbedding #> mainExpName: NULL #> altExpNames(0):

投影的细胞周期空间将被存储在reducedDims使用名称" tricycleEmbedding "(你可以设置其他嵌入名称)。

库(ggplot2)库(scattermore)库(scater) scater::plotReducedDim(neurosphere_example, dimred = "tricycleEmbedding") +实验室(x = "投影PC1"， y = "投影PC2") + ggtitle(sprintf("投影细胞周期空间(n=%d)"， ncol(neurosphere_example))) + theme_bw(base_size = 14))

5推断细胞周期位置

一旦新数据被投影到细胞周期嵌入中，就可以使用estimate_cycle_position ()．如果数据没有被投影，这个函数将为您进行投影。假设SingleCellExperiment作为输入，单元格周期位置将被添加到colData对象的名称tricyclePosition．

neurosphere_example <- estimate_cycle_position(neurosphere_example) names(colData(neurosphere_example)) #> [1] "TotalUMIs" "sample" "tricyclePosition"

估计的细胞周期位置在0和2之间。注意，我们努力获得高分辨率的细胞周期状态，我们认为连续位置在描述细胞周期时更合适。但是，为了便于用户理解位置变量，我们也注意到，用户可以将0.5pi近似为S阶段的开始，pi近似为G2M阶段的开始，1.5pi近似为M阶段的中间，1.75pi-0.25pi近似为G1/G0阶段。

6评估性能

我们有两种方法(快速)评估TriCycle是否有效。他们是

  1. 看看数据在细胞周期嵌入中的投影。
  2. 看看关键基因的表达作为细胞周期位置的函数。

如上所示，将数据投影到单元周期嵌入中。我们观察到，更深的序列数据将有一个更清晰的椭球图案与一个空的内部。随着排序深度的减小，椭球的半径减小，直到空的内部消失。因此，没有内部结构并不意味着这个方法不管用。

更重要的是检查几个基因作为细胞周期位置的函数。我们倾向于使用Top2a，它是高度表达的，因此在每个数据集中都是“可绘制的”。其他候选的例子是Smc2。为了绘制这些数据，我们提供了一个方便的函数fit_periodic_loess ()以黄土线之间的拟合为循环变量\θ(\ \)和其他响应变量。这种装配是通过制作来完成的theta.v3期(c(θ。V - 2 *。v,θ。V + 2 * pi))和重复y3次。只有与原值相对应的拟合值theta.v将被返回。在这个例子中，我们展示了如何很好地表达细胞周期标记基因Top2a改变一起\θ(\ \)．

top2a。idx <- where (rowData(neurosphere_example)$Gene == 'Top2a')匹配。l <- fit_periodic_黄土(neurosphere_example$tricyclePosition, assay(neurosphere_example， 'logcounts')[top2a.]idx，]， plot = TRUE, x_lab ="细胞周期位置\u03b8"， y_lab =" log2(Top2a)"， fig.title = paste0(" Top2a沿\u03b8的表达式(n="， ncol(neurosphere_example)， ")")) names(fit.l) #>[1] "拟合" "残距" "pred. "df”“黄土。O " "r平方" "无花果"适合。L $fig + theme_bw(base_size = 14)

对于Top2a，我们预计峰值表达在\π(\ \)．

7替代方案:推断细胞周期阶段

该方法是由Schwabe等人(2020)．我们做了一些小的修改来减少NA作业。但平均而言，性能与最初的实现非常相似Revelio包中。简而言之，我们计算z-数十个高表达阶段特异性细胞周期标记基因，并将细胞分配到最大的阶段z分数。

neurosphere_example <- estimate_Schwabe_stage(neurosphere_example, gname. sh)type = 'ENSEMBL'， species = 'mouse') #>这个函数是Schwabe等人2020的重新实现。如果你想使用三轮车方法，请运行estimate_cycle_position!#>不输入gname。将使用x的行名。没有指定注解db。org.Mm.eg.db将用于将小鼠ENSEMBL id映射到基因符号。#>批次1阶段G1。S基因:50 #> Batch 1 phase S基因:43 #> Batch 1 phase G2基因:61 #> Batch 1 phase G2M基因:67 #> Batch 1 phase M.G1基因:30 scater::plotReducedDim(neurosphere_example, dimred =" tricycleEmbedding"， colour_by =" CCStage") +实验室(x ="投影PC1"， y ="投影PC2"， title = paste0("投影细胞周期空间(n="， ncol(neurosphere_example)， ")) + theme_bw(base_size = 14))

8画出核密度

另一个有用的函数是plot_ccposition_den的核密度\θ(\ \)以冯·米塞斯分布的表现型为条件。输出图形以极坐标和笛卡尔坐标两种形式提供。这在比较不同的细胞类型、治疗方法或只是阶段时很有用。(因为我们在这里使用了一个非常小的数据集作为例子，我们设置了带宽，即冯米塞斯分布的浓度参数为10，以得到一条平滑的线。)

plot_ccposition_den(neurosphere_example$tricyclePosition, neurosphere_example$sample， 'sample'， bw = 10, fig.title = " \u03b8的内核密度")+ theme_bw(base_size = 14)

plot_ccposition_den(neurosphere_example$tricyclePosition, neurosphere_example$sample， 'sample'， type = "circular"， bw = 10, fig.title = " \u03b8内核密度")+ theme_bw(base_size = 14)

9绘制按细胞周期位置着色的嵌入散射图

以可视化细胞周期位置\θ(\ \)在任何嵌入上，我们都需要仔细选择一个循环的调色板。因此，我们包含这样的函数来绘制任意的嵌入SingleCellExperiment具有循环变量的对象。还有一个用于创建循环图例的辅助函数。

库(cowplot) p <- plot_emb_circle_scale(neurosphere_example, dimred = 1，点。大小= 3.5，点。Alpha = 0.9) + theme_bw(base_size = 14) legend <- circle_scale_legend(text. size = 14)plot_grid(p, legend, ncol = 2, rel_width = c(1, 0.4))

我们画出投影嵌入。在实践中，用户可以使用其他嵌入，如UMAP或t-SNE，并获得信息表示。

10创建一个新的引用

用户可以通过对细胞周期基因进行主成分分析来进行自己的参考，并在其他函数中使用学习到的旋转矩阵作为参考矩阵。这里有一个例子run_pca_cc_genes功能提取基因本体细胞周期基因(去:0007049)并运行PCA。通过用学习到的参考投影数据本身，这些投影等同于直接的PCA结果。但是您可以使用这个新学习到的参考来投影其他数据集。

gocc_sce set.seed(100)。o <- run_pca_cc_genes(neurosphere_example, express_values = "logcounts"， species = "mouse") #> No AnnotationDb指定。org.Mm.eg.db将用于将小鼠ENSEMBL id映射到基因符号。#>不输入gname。sce的行名。将使用O。#> 'select()'返回1:键和列之间的许多映射#>在你的数据中发现的3024个基因本体论细胞周期基因中的609个。新的。ref <- attr(reducedDim(gocc_sce。o， 'PCA')， 'rotation')[， seq_len(2)] head(new.ref) #> PC1 PC2 #> ENSMUSG00000040204 0.22243912 0.15499421 #> ENSMUSG00000023067 - 0.0295683 -0.16281149 #> ENSMUSG00000007872 -0.01837872 -0.02923750 #> ENSMUSG00000020649 0.18210797 0.14993867 #> ENSMUSG00000001403 0.18395862 -0.20999483 #> ENSMUSG00000060860 0.04130150 -0.07429779 new_sce <- estimate_cycle_position(neurosphere_example, ref.m = new。ref, dimred = 'tricycleEmbedding2') #>指定的dimred在singlecel实验或altexp中不存在。将运行project_cycle_space来计算embedding tricycleEmbedding2 #>在新数据中找到的投影基因的数量是500。

注意:如果用户要计算两个循环变量之间的相关性，如细胞周期位置，传统的pearson相关系数不会考虑周期性。用户可以使用(绝对的)循环相关值。(通过PCA重新学习参考时，PC1和PC2的符号是不确定的。如果PC1翻转，会有一个\π(\ \)转变。PC2也是如此。如果用户修复了引用，就不会有任何翻转。但考虑到周围的变化\ (0 \)或π\ [2 \ \)，仍应使用圆形相关系数代替pearson相关系数。)

cor(neurosphere_example$tricyclePosition, new_sce$tricyclePosition) #> [1] 0.3179227 CircStats::circ。cor(neurosphere_example$tricyclePosition, new_sce$tricyclePosition) #> r #> 1 -0.9698397 qplot(x = neurosphere_example$tricyclePosition,y = new_sce$tricyclePosition) +实验室(x ="原\u03b8"， y ="新\u03b8"， title = paste0("两个\u03b8的比较(n="， ncol(neurosphere_example)， ")) + theme_bw(base_size = 14))

11使用具有批处理效果的数据集创建新的引用

本节介绍如何使用具有批处理效果的数据集创建新的引用。仅建议在数据中识别批处理效果并希望使用该数据构建自定义引用的专家用户使用。理论上，用户可以使用其他方法来消除批效应。这里，我们使用修拉，用来构建我们的神经圈参考(Zheng et al. (2022))，作为例子。(不计算本节中的代码。)

#假设我们有一个计数矩阵包含所有批次的单元格;我们首先将矩阵子集化为GO细胞周期基因库(org. mm . e.g. .db)库(AnnotationDbi) cc.genes <- AnnotationDbi::select(org. mm . e.g. .db, keytype = "GOALL"， keys = "GO:0007049"， columns = "ENSEMBL")[， "ENSEMBL"] count_cc。M <- count.m[ensemble . M]id %in% cc.genes，] # ensembl。Ids是集合。每行计数的id。然后，我们使用子集矩阵构造一个Seurat对象，并设置批处理变量库(Seurat) Seurat。o <- CreateSeuratObject(counts = count_cc.m) seurat。O [["batch"]] <- batch。v #做一个Seurat列表，并为每批单独规范化#变量特征定义需要FindIntegrationAnchors函数Seurat。list <- lapply(SplitObject(seurat. list)o,分裂。by = "batch")， function(x) FindVariableFeatures(NormalizeData(x))) #查找锚点并合并数据seurat。- FindIntegrationAnchors(对象。list = seurat.list) seurat.integrated <- IntegrateData(anchorset = seurat.anchors) corrected。m <- seurat.integrated@assays$integrated@data #在批处理效果校正矩阵上运行PCA，并获得前2个PCs PCA的旋转分数。m <- scater::calculatePCA(已更正。M, ntop = 500) new。Ref <- attr(pca。M， 'rotation')[， seq_len(2)]

会话信息

#> 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。所以#> LAPACK: /home/biocbuild/bbs-3.16-bio /R/lib/libRlapack。so #> #> 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_TELEPHONE= c# > [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION= c# > #>附加基础包:#> [1]stats4 stats graphics grDevices utils datasets methods #>[8]基础#> #>其他附加包:#> [1] cowplot_1.1.1 scater_1.26.0 #> [3] sccuttle_1 .8.0 scattermore_0.8 #> [5] ggplot2_3.3.6 tricycle_1.6.0 # b> [7] SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0 #> [9] Biobase_2.58.0 genome icranges_1 .50.0 #> [11] GenomeInfoDb_1.34.0 IRanges_2.32.0 #> [13] S4Vectors_0.36.0 BiocGenerics_0.44.0 #> [15] MatrixGenerics_1.10.0 matrixStats_0.62.0 #> [17] BiocStyle_2.26.0 #> #>通过命名空间加载(并且没有附加):#> [1] CircStats_0.2-6 bitops_1.0-7 #> [3] bit64_4.0.5 rcolorbrewer_1 .3 #> b[5] httr_1.4.4 tools_4.2.1 #> [7] bslib_0.4.0 utf8_1.2.2 #> [9] R6_2.5.1 irlba_2.3.5.1 #> [11] vipor_0.4.5 DBI_1.1.3 #> [13] colorspace_2.0-3 withr_2.5.0 #> [15] gridExtra_2.3 tidyselect_1.2.0 #> [17] bit_4.0.4 compiler_4.2.1 #> [19] cli_3.4.1 BiocNeighbors_1.16.0 #> [23] bbbdelayedarray_0.24.0 labeling_0.4.2 #> [25] scales_1.2.1 mvtnorm_1.1-3 #> [27] string_1 .4.1 digest_0.6.30 #> [29]rmarkdown_2.17 XVector_0.38.0 # > [31] pkgconfig_2.0.3 htmltools_0.5.3 # > [33] sparseMatrixStats_1.10.0 highr_0.9 # > [35] fastmap_1.1.0 rlang_1.0.6 # > [37] RSQLite_2.2.18 circular_0.4 - 95 # > [39] DelayedMatrixStats_1.20.0 farver_2.1.1 # > [41] jquerylib_0.1.4 generics_0.1.3 # > [43] jsonlite_1.8.3 BiocParallel_1.32.0 # > [45] dplyr_1.0.10 rcurl_1.98 - 1.9 # > [47] magrittr_2.0.3 BiocSingular_1.14.0 # > [49] GenomeInfoDbData_1.2.9 Matrix_1.5-1 # > [51] ggbeeswarm_0.6.0 Rcpp_1.0.9 # > [53] munsell_0.5.0fansi_1.0.3 #> [55] viridis_0.6.2 ggng_0.1 -7 #> [61] zlibbioc_1.44.0 grid_4.2.1 #> [63] blob_1.2.3 ggrepel_0.9.1 #> [65] parallel_4.2.1 crayon_1.5.2 #> [67] lattice_0.20-45 Biostrings_2.66.0 #> [69] beachmat_2.14.0 KEGGREST_1.38.0 #> [71] magick_2.7.3 knitr_1.40 #> [73] pillar_1.8.1 boot_1.3-28 #> [75] codetools_0.2-18 ScaledMatrix_1.6.0 #> [79] BiocManager_1.30.19 png_0.1-7 #> [81]vctrs_0.5.0 org.Mm.eg.db_3.16.0 #> [83] gtable_0.3.1 assertthat_0.2.1 #> [85] cachem_1.0.6 xfun_0.34 #> [87] rsvd_1.0.5 mime_0.12 #> [89] viridisLite_0.4.1 tibble_3.1.8 #> [91] beeswarm_0.4.0 AnnotationDbi_1.60.0 #> [93] memoise_2.0.1

12参考文献

在包中，我们提供了一个参考，从小鼠神经圈RNAseq的完整数据集中学习。参考文献给出了500个细胞周期基因的权重及其id。虽然是从老鼠身上学来的，但它也适用于人类数据，通过基因符号来映射基因。

data(neuroRef) head(neuroRef) #> pc1。腐烂pc2。> ENSMUSG00000040204.6 -0.1989885 0.11245580 ENSMUSG00000040204 Pclaf Pclaf #> ENSMUSG00000020914.17 -0.1878494 0.04252486 ENSMUSG00000020914 Top2a Top2a #> ENSMUSG00000060860.8 -0.0567803 -0.07531686 ENSMUSG00000060860 Ube2s Ube2s #> ENSMUSG00000001403.13 -0.1662416 -0.15992459 ENSMUSG00000001403 Ube2c Ube2c #> ENSMUSG00000026605.14 -0.1721841 -0.11416861 ENSMUSG00000026605 Cenpf Cenpf #> ENSMUSG00000029177.9 -0.1357597 -0.19290262 ENSMUSG00000029177 Cenpa Cenpa

Schwabe, Daniel, Sara Formichetti, Jan Philipp Junker, Martin Falcke和Nikolaus Rajewsky, 2020。细胞周期中单细胞的转录组动力学分子系统生物学16 (11): e9946。https://doi.org/10.15252/msb.20209946．

郑世杰C.， Genevieve Stein-O 'Brien, Jonathan J. Augustin, Jared Slosberg, Giovanni A. Carosso, Briana Winer, Gloria Shin, Hans T. Bjornsson, Loyal A. Goff和Kasper D. Hansen. 2022。“基于迁移学习的细胞周期位置通用预测”基因组生物学23: 41。https://doi.org/10.1186/s13059-021-02581-y．