1简介

的censcyt包提供了功能的差异丰度分析，当一个协变量是受权利审查。扩展了的差异发现功能diffcyt包(Weber et al. 2019)通过包括细胞群体丰度和一个删减变量之间的关联测试，如生存时间。一般工作流程与困难的包装小插图在继续之前最好先熟悉一下。

作为一个简短的总结，diffcyt工作流程包括原始细胞分析数据的预处理，通过自动聚类和差异丰度分析(或差异状态分析)分配细胞类型(Weber et al. 2019)．如果与细胞术样本相关的时间到事件变量(例如生存时间)是可用的，人们可能会对细胞种群丰度和这个时间到事件变量之间的关系感兴趣。在实践中，时间到偶数变量的一个常见问题是，并非所有事件都被观察到。相反，可能只有一个最小的时间是已知的，换句话说，这些值是(正确的)审查。如果使用标准分析(如testDA_GLMM或testDA_edgeR从diffcyt)在参数估计中引入偏差。最简单的解决方法是排除所有不完整的样本，这被称为完整案例分析或按列表删除。但在某些情况下，例如高审查率或非随机审查，这种方法表现不佳，因为它会大大降低统计能力或引入偏差。

通常，当一个变量被删除时，可以使用经典的生存分析方法(例如Cox比例风险模型)进行关联测试。但是在差异丰度分析的背景下diffcyt删减变量不是响应，而是协变量(与经典的生存分析相反)。使用的方法censcyt克服这一障碍的方法是多重归责鲁宾(1988)这允许以增加运行时间为代价使用标准的差分分析方法。多重归责包括三个步骤归责，分析而且池这些都是在censcyt．

在第一步(归责)通过从预测分布中随机抽取值来替换缺失(或删减)值，从而生成多个完整的数据集(见表1)．然后在第二步(分析)，在每个估算数据集上拟合一个GLMM，其中细胞种群计数被建模为响应，而截尾(或更确切地说是估算的)变量(例如生存时间)作为协变量。在最后一步(池)第二步的结果结合使用鲁宾斯规则鲁宾(1988)这就考虑到了来自归因的额外不确定性。

在函数中主要可用三种插补方法testDA_censoredGLMM均列于表1．

如果最大的值被删减，生存函数就不会达到它的理论最小值0，这就会留下一些未定义的替换。在Kaplan-Meier imputation中实现了多个处理此问题的选项(表2)．

2数据

首先，加载必要的包。重要的是要注意加载censcyt也会载入diffcyt．

suppressPackageStartupMessages({library(censcyt) library(ggplot2) library(summarizeexperiment) library(tidyr)})

DA分析的主要输入包括细胞群体丰度的聚类乘以样本矩阵。如何获得这些计数在diffcyt装饰图案这里就不再进一步讨论了。相反，为了说明目的，我们模拟了一个数据集。在这种情况下，删减协变量是根据指数分布建模的。为了获得截尾，还模拟了一个附加变量截尾时间，实际测量值为这两个变量中的最小值。

set.seed(1234) nr_samples <- 50 nr_clusters <- 6 #协变量从指数分布中模拟:X_true <- rexp(nr_samples) #检查时间也从指数分布中采样:C <- rexp(nr_samples) #实际观测值是两者中的最小值:X_obs <- pmin(X_true,C) #此外，我们有事件指示器:Event_ind <- ifelse(X_true>C,0,1) #比例检查:(length(Event_ind)-sum(Event_ind))/length(Event_ind)

## [1] 0.52

接下来，细胞计数根据狄利克雷多项分布建模，其中一些集群与协变量有关联。这个函数simulate_multicluster可用于此。

#每个样本容量的单元数<- round(runif(nr_samples,1e3,1e4)) #狄利克雷多项分布的alpha参数。alphas <- runif(nr_clusters,1,10) #主要模拟函数simulation_output <- simulate_multicluster(alphas = alphas, sizes = sizes, covariate = X_obs, nr_diff = 2, return_summarized_experiment = TRUE, slope = list(0.9)) #作为一个summarizeexperiment对象d_counts <- simulation_output$counts

为了了解这个模拟数据的样子，我们可以绘制每个样本的相对细胞群体丰度与每个细胞群体的删减协变量。

is_diff_cluster <- ifelse(is.na(simulation_output$row_data$paired)，FALSE,TRUE) #第一次转换为比例:比例<- as.data.frame(t(apply(assay(d_counts)，2,function(x)x/sum(x)))) names(比例)<- paste0("Nr:"，names(比例)，" - "，ifelse(is_diff_cluster，"DA"，"non DA")) # then to long format for counts_long <- pivot_long (proportion, cols= seq_len(nr_clusters)， names_to = "cluster_id"， values_to = "proportion") # add observed(部分检查)变量counts_long$X_obs <- rep(X_obs,each=nr_clusters) # add event indicator counts_long$Event_ind <- factor(rep(Event_ind,each=nr_clusters)， levels = c(0,1)，labels=c(" deleted "，"observed")) ggplot(counts_long) + geom_point(aes(x=X_obs,y=proportion,color=Event_ind)) + facet_wrap(~cluster_id)

图1:相对细胞群丰度与存活时间(X_obs)
Event_ind表示X_obs是否被截尾或观察。

3.建立元数据

下一步是设置元数据，即测试中使用的协变量。这和在diffcyt但另外，必须给出截尾协变量的事件指示符，因为截尾变量由两个值表示，测量值本身和测量值(观察值(1)或截尾(0))时的指示符。

#所有的变量和样本id实验者信息<- data.frame(sample_id = seq_len(nr_samples)， X_obs = X_obs, Event_ind = Event_ind)

的createFormula函数与中使用的函数类似diffcyt但是除此之外event_indicator必须通过在中给出各自的列名来指定experiment_info．如果在参数中给出多个协变量cols_fixed然后event_indicator与第一个给定的元素有关cols_fixed．

da_formula <- createFormula(experiment_info = experiment_info, cols_fixed = "X_obs"， cols_random = "sample_id"， event_indicator = "Event_ind") #也创建对比矩阵测试对比<- matrix(c(0,1))

4微分测试

主要部分是使用函数进行差分测试testDA_censoredGLMM．它包含了与未审查版本相同的论点testDA_GLMM(从diffcyt)，并附带两个特定于多重imputation的附加参数。第一个附加参数是mi_reps(多次输入重复)，即执行多个输入步骤的次数。不幸的是，目前还没有明确的规则来确定需要多少次估算，只能凭经验来判断(Buuren 2018)．一般来说，如果审查率高，在需要高统计能力的应用中，需要更多的imputation(Buuren 2018)．一个“规则”是的二次规则Hippel (2018)哪一个可以显示如果详细的设置为真正的．第二个论点是imputation_method用于指定imputation方法。见表1和表2．

# test with 50 repeat with method risk set imputation (rs) da_out <- testDA_censoredGLMM(d_counts = d_counts, formula = da_formula, contrast = contrast, mi_rep = 30, imputation_method = "km"， verbose = TRUE, BPPARAM=BiocParallel::MulticoreParam(workers = 1))

50个值中的26个被审查

##建议的imputation数量:68

看看我们可以用的结果diffcyt的topTable功能:

topTable (da_out)

6行3列的数据帧## cluster_id p_val p_adj ## <字符> <数字> <数字> ## 33 0.000898221 0.00538933 ## 6 6 0.005527278 0.01658183 ## 1 1 0.807657070 0.95288779 ## 22 0.952887790 0.95288779 ## 4 4 0.813003970 0.95288779 ## 55 0.703777553 0.95288779

is.na(simulation_output$row_data$paired))

## [1] 3 6

5包装器函数

与其使用上面所示的单个函数，还可以使用包装器函数censcyt哪个类似于diffcyt包装。首先，模拟一些表达式数据，即一个比上面例子中更真实的起始位置。

#函数来创建表达式的数据(样本)fcs_sim < -函数(n = 2000,意味着= 0,sd = 1, ncol = 10,咖啡= 5){d <——矩阵(sinh (rnorm (n * ncol,意思是,sd)) *咖啡,ncol = ncol) #每一个标记在一个人口严重表示,即这应该结果#“ncol”集群(我在seq_len (ncol)) {d (seq (n / ncol *(张)+ 1,n / ncol * (i)),我]< - sinh (rnorm (n / ncol,意味着+ 2,sd)) *咖啡}colnames (d) < - paste0(“标记”,sprintf(“% 2 d”,1:ncol)) d} #用DA信号创建多个表达式数据样本comb_sim <- function(X, nr_samples = 10, mean = 0, sd = 1, ncol = 10, cof = 5){#创建随机数据(没有差分信号)d_input <- lapply(seq_len(nr_samples)， function(i) fcs_sim(mean = mean, sd = sd, ncol = ncol，咖啡(咖啡)#加微分丰富(DA)信号(我在seq_len (nr_samples)){#集群中的细胞数量1 n_da < -轮((plogis (X_true[我]))* 300)#将随机表达式d_input[[我]][seq_len (n_da)] < -矩阵(sinh (rnorm (n_da * ncol,意思是,sd)) *咖啡,ncol = ncol) #增加集群1 d_input表达式[[我]][seq_len (n_da), 1] < - sinh (rnorm (n_da,意味着+ 2,sd)) *咖啡}d_input} #设置参数和装病者d_input < - comb_sim (X_true nr_samples = 50)

的对象experiment_info，da_formula而且对比与上面指定的相同，但必须提供有关测量标记的额外信息。

markker_info <- data.frame(channel_name = paste0("channel"， sprintf("%03d"， seq_len(10)))， markker_name = paste0("marker"， sprintf("%02d"， seq_len(10)))， marker_class = factor(c(rep("type"， 10))， levels = c("type"， "state"， "none"))， stringsAsFactors = FALSE)

最后，可以运行包装器函数。论点mi_reps是指定重复次数在多重imputation和imputation_method允许设置imputation方法。为了加快计算速度，aBiocParallelParam对象(如。MulticoreParam，从包装BiocParallel)可以用于论证BPPARAM使部分计算并行化。

#运行包装器函数out_DA <- censcyt(d_input, experiment_info, marker_info, formula = da_formula, contrast = contrast, meta_clustering = TRUE, meta_k = 10, seed_clustering = 123, verbose = FALSE, mi_reps = 10, imputation_method = "km"， BPPARAM=BiocParallel::MulticoreParam(workers = 1))

## FlowSOM集群在5.9秒内完成

topTable(out_DA, format_vals = TRUE)

## 10行3列的数据帧## cluster_id p_val p_adj ## <字符> <数字> <数字> ## 8 8 0.00142 0.0142 ## 3 3 0.02110 0.1060 ## 7 7 0.31900 0.7970 ## 10 10 0.31900 0.7970 ## 11 0.90300 0.9030 ## 2 2 0.58500 0.9030 ## 4 4 0.76200 0.9030 ## 5 5 0.83400 0.9030 ## 9 9 0.85400 0.9030 ##

(out_DA, analysis_type = "DA"， sample_order = order(X_true))