内容

1安装

bayNorm已经提交给Bioconductor，一旦被接受，可以通过以下方式进行安装:

库(BiocManager) BiocManager::安装(“bayNorm”)

目前bayNorm可以通过github安装:

库(devtools) devtools: install_github(“WT215 / bayNorm”)

2快速入门:用于一组或多组单元格

主要功能是bayNorm它是先验参数估计和规格化数组或矩阵生成的包装函数。

运行基本参数bayNorm是:

输入参数BETA_vec，条件(如果指定了的话),UMI_sffl(如果指定了的话),Prior_type，FIX_MU，BB_SIZE而且GR存储在列表中input_params应该从哪个输出bayNorm．对象先知先觉在一起input_params输出bayNorm应该输入bayNorm_sup用于在规格化计数矩阵的三维阵列、模态或平均版本输出之间进行转换。

2.1捕获效率的估计

除了原始数据外，用户可能需要提供的另一个参数是平均捕获效率\(β> < \ \)因此β\ (\ \)可以用任何其他方法估计的比例因子进一步计算。默认的β\ (\ \)计算的总计数归一化为6%。或者你可以用BetaFun在bayNorm中估计β\ (\ \)：

数据(EXAMPLE_DATA_list) #假设输入数据是SummarizedExperiment对象:#这里30细胞用于说明证交所< - SummarizedExperiment:: SummarizedExperiment(化验= SimpleList(数量= EXAMPLE_DATA_list inputdata美元[,seq (30)])) # SingleCellExperiment bayNorm对象也可以输入:#证交所< - SingleCellExperiment:: SingleCellExperiment(化验=列表(数量= EXAMPLE_DATA_list inputdata美元))BETA_est < -BetaFun (data =证交所,MeanBETA = 0.06)总结(BETA_estβ美元)
最小第一曲，中位数，平均第三曲，最大值。## 0.02299 0.03859 0.05255 0.06000 0.07473 0.11168

这个函数BetaFun选择一个基因子集，排除细胞中异常基因和高比例零的基因。然后将每个细胞中基因子集的总数归一化为\(β> < \ \)．另一种方法是将R包估计的比例因子归一化食物来\(β> < \ \)．

2．2运行bayNorm

BETA_vec可以设置为零，则β\ (\ \)估计总数将正常化至6%

#返回三维数组归一化数据，从后验分布中抽取20个样本:bayNorm_3D<-bayNorm(data =rse, BETA_vec = NULL, mode_version=FALSE, mean_version =FALSE, S=20,verbose =FALSE, parallel = TRUE)
# # | | | 0 % | |==== | 5 % | |======= | 10 % | |========== | 15 % | |============== | 20 % | |================== | 25 % | |===================== | 30 % | |======================== | 35 % | |============================ | 40 % | |================================ | 45 % | |=================================== | 50 % | |====================================== | 55 % | |========================================== | 60 % | |============================================== | 65% ||================================================= | 70年 % | |==================================================== | 75年 % | |======================================================== | 80年 % | |============================================================ | 85年 % | |=============================================================== | 90年 % | |================================================================== | 95年 % | |======================================================================| 100%
#返回二维矩阵归一化数据(后验的MAP): #简单设置mode_version=TRUE，但保留mean_version=FALSE #返回二维矩阵归一化数据(后验的均值):#简单设置mean_version=TRUE，但保留mode_version=FALSE

2．3非umi scRNAseq数据集

bayNorm的数学模型适用于UMI数据集。但它也可以应用于非umi数据集。在bayNorm时，需要指定以下参数:*UMI_sffl: bayNorm也可以应用于非umi数据集。但是，用户需要提供一个比例的数字。原始数据将被缩放数分割，并对四舍五入的缩放数据应用bayNorm。通过这样做，Dropout与Mean表达式图将类似于UMI数据集。

2.4利用现有估计的先验参数生成三维数组或二维矩阵。

如果你之前在数据集上运行过bayNorm，但想要输出另一种类型的数据(3D数组或2D矩阵)，你可以使用这个函数bayNorm_sup．中指定以下参数来输入现有的估计参数是很重要的bayNorm_sup：

  • BETA_vec:如果条件之前已经指定，然后输入unlist (bayNorm_outputβ美元)
  • 先知先觉：输入bayNorm_output PRIORS_LIST美元
  • 条件:请确保指定相同的参数条件像以前一样。您可以从bayNorm函数之前的输出中找到这两个对象，bayNorm函数是一个列表。
#现在如果你想生成2D矩阵(MAP)使用相同的先验#估计之前生成:bayNorm_2D<-bayNorm_sup(数据=rse, PRIORS=bayNorm_3D$PRIORS, input_params=bayNorm_3D$input_params, mode_version=TRUE, mean_version =FALSE, verbose =FALSE) #或者你可能想生成2D矩阵#(后验平均值)使用相同的先验#估计之前生成:bayNorm_2D<-bayNorm_sup(Data=rse, PRIORS=bayNorm_3D$PRIORS, input_params=bayNorm_3D$input_params, mode_version=FALSE, mean_version = TRUE, verbose =FALSE)

3.其他功能

  1. Prior_fun的包装函数EstPrior，BB_Fun而且AdjustSIZE_fun：
    • EstPrior:使用MME方法估计先验。
    • BB_Fun:通过最大化两者的边际似然函数来估计先验φ\ (\ \)或两个\μ(\ \)而且φ\ (\ \)．
    • AdjustSIZE_fun:调整φ\ (\ \)估计从BB_Fun．

  2. noisy_gene_detection:这是一个比较高级的包装器函数bayNorm，bayNorm_sup而且SyntheticControl．有关此函数背后的原理的详细信息，请参阅本文中的方法部分2．

  3. SyntheticControl给出一个真实的scRNA-seq数据β\ (\ \)，它将估计\μ(\ \)而且φ\ (\ \)并从中模拟控制数据的使用泊松分布。参见论文中的方法部分2欲知详情。

#下游分析:bayNorm论文中用到的DE基因检测DE函数，由SCnorm作者提供。

library(MAST) library(reshape2) ########SCnorm_runMAST3########## SCnorm_runMAST3 <- function(Data, NumCells) {if(length(dim(Data))==2){resultss<-SCnorm_runMAST(Data, NumCells)} else if(length(dim(Data))==3){library(foreach) resultss<- foreach(sampleind=1:dim(Data)[3]，.combine=cbind)%do%{print(sampleind) qq<-SCnorm_runMAST(Data[，，sampleind]， NumCells) return(qq$adjpval)}} return(resultss)} SCnorm_runMAST <- function(Data，NumCells) {NA_ind <——(rowSums(数据)= = 0)Data = as.matrix (log2(数据+ 1))G1 <拼,seq [(1, NumCells [1])] G2 <拼,seq [(1, NumCells [1])] qq_temp < - rowMeans (G2) -rowMeans (G1) qq_temp [NA_ind] = NA numGenes =暗(数据)[1]datalong =融化(数据)气孔导度= c(代表(c1, NumCells [1] * numGenes),代表(c2, NumCells [2] * numGenes))按日计算工资= cbind (datalong电导率)colnames(包含日期的)= c(“基因”,“细胞”、“价值”、“电导率”)按日计算工资$基因=因素(包含日期的基因美元)按日计算工资$细胞=因素(包含日期的单元格美元)vdata = FromFlatDF (dataframe =按日计算工资,idvars = "cell"， primerid = "gene"， measurement = "value"， id = numeric(0)， cellvars = "Cond"， featurevars = NULL, phenovars = NULL) zlm。zlm(~ Cond, vdata, method='glm'， ebayes=TRUE)lr = lrTest(zlm。output， 'Cond') gpval = zlm.lr[，，'Pr(>Chisq)'] adjpval = p.adjust(gpval[，1]，"BH") ##只使用连续部分的pvalue。adjpval=adjpval [rownames(Data)] return(list(adjpval=adjpval, logFC_re=qq_temp))}
#现在检测两组细胞之间的DE基因，每组15个细胞#对于3D数组DE_out<-SCnorm_runMAST3(Data=bayNorm_3D$Bay_out, NumCells=c(15,15)) med_pvals<-apply(DE_out,1,median) #DE基因的调用阈值为0.05:DE_genes<-names(which(med_pvals<0.05)) #对于2D数组DE_out<-SCnorm_runMAST3(Data=bayNorm_2D$Bay_out, NumCells=c(15,15)) DE_genes<-names(which(DE_out$adjpval<0.05))

4方法

4．1bayNorm的基本原理

scRNA-seq数据集通常用维度矩阵表示问\ \ (P \倍)，其中P表示观察到的基因总数，Q表示研究的细胞总数。的元素\(间{ij} \)（\(i \in \{1,2，\ldots, P\}\)而且\(j \in \{1,2，\ldots, Q\}\))，表示报告的转录本数目\(^文本{th}} {\ \)基因\ (j ^文本{th}} {\ \)细胞。对于非umi方案，这等于映射到该细胞中该基因的测序读取的总数。对于基于UMI的协议，这等于映射到每个基因的单个UMI的数量。该基质可以包括来自不同组或批次细胞的数据，代表不同的生物条件。这可以表示为单元格组或条件的标签向量(\ (C_j \))．bayNorm为每个细胞(j)中的每个基因(i)生成原始表达计数(\ \(间{ij} ^ 0))，考虑到该基因的scRNA-seq读取计数(\(间{ij} \))．

标准化scRNA-seq数据的常用方法是使用全局缩放因子(\ (s_j \))，忽略任何特定基因的偏差。规范化的数据\ \(波浪号{x} _ {ij} \)是通过划分每个单元格的原始数据得到的吗\ (j \)通过它的整体比例因子\ (s_j \)：

{方程}\ \[\开始标签{装备:缩放}\波浪号{x} _ {ij} = \压裂{间{ij}} {s_j} \{方程}结束\]

在bayNorm中，我们使用贝叶斯方法实现全局扩展。我们假设给定细胞中转录本的原始数量(\ \(间{ij} ^ 0))，所观察到的转录本数目(\(间{ij} \))遵循具有概率的二项式模型\ \ beta_j \ ()，我们称之为捕获效率，它表示细胞中原始转录本被观察到的概率。另外，我们假设原数还是真数\(^文本{th}} {\ \)基因\ (j ^文本{th}} {\ \)电池(\ \(间{ij} ^ 0))服从负二项分布，参数均值为(\μ(\ \))，尺寸(或色散参数，φ\ (\ \))，以便:\[\公关(x ^ 0 _ {ij} = n | \ phi_i \ mu_i) = \压裂{\γ(n + \ phi_i)}{\γ(\ phi_i) n !}(\压裂{\ phi_i} {\ mu_i + \ phi_i}) ^ {\ phi_i}(\压裂{\ mu_i} {\ mu_i + \ phi_i}) ^ {n} \]

总的来说，我们有以下模型:

\[开始\{方程}{模型:湾}\ \标签开始{分裂}间{Binom} {ij} \ sim \文本(间{ij} ^ 0,文本{概率}= \ \ beta_j) \ \间{ij} ^ 0 \ sim \文本{NB}({意味着}= \ mu_i \文本,文本{大小}= \ \ phi_i)结束\{分裂}\{方程}结束\]

利用贝叶斯规则，我们得到了每个细胞中每个基因的mrna原始数量的后验分布:

方程\[开始\ {}\ underbrace{\公关(间{ij} ^ 0 |间{ij}, \ beta_j \ mu_i \ phi_i)} _{后}= \ \文本dfrac {\ overbrace{\公关(间{ij} |间{ij} ^ 0 \ beta_j)} ^ \文本{可能性}\ * \ overbrace{\公关(间{ij} ^ 0 | \ mu_i \ phi_i)} ^ \文本之前}{}{\ underbrace{\公关(间{ij} | \ mu_i \ phi_i \ beta_j)} _ \文本{边际似然}}\{方程}结束\]

先验参数\μ(\ \)而且φ\ (\ \)的补充资料中详细讨论了使用经验贝叶斯方法估计每个基因的1．

有关bayNorm基本原理的更多详细信息，请参见1．

4.2捕获效率的估计

细胞特异性捕获效率\ \ beta_j \ ()而全局比例因子(\ (s_j \))是密切相关的。我们可以将不同方法估计的比例因子(见下文)转换为\ \ beta_j \ ()取值为:

\[开始\{方程}\ beta_j = (s_j / \酒吧{年代})\酒吧{β\}\{方程}结束\]

\(\酒吧{β\}\)或\(β> < \ \)，一个标量，是所有单元的全局平均捕获效率的估计值，范围在0到1之间。

有两种不同的估算方法\(\酒吧{β\}\)而且\ \ beta_j \ ()：

  1. 如果spike-ins或smFISH数据可用，它们可用于估计捕获效率。我们可以将每个单元中观察到的spike-ins总数除以输入的spike-ins总数，或者我们可以拟合线性回归来估计单元特异性\ \ beta_j \ ()．如果有smFISH数据，我们可以在原始数据的平均表达式(响应变量)和smFISH数据的平均表达式(解释变量)之间拟合线性回归。解释变量的系数可用作\(\酒吧{β\}\)．
  2. 原始数据本身可直接用于估计单元特定的整体比例因子(\ (s_j \))．然后估计\(\酒吧{β\}\)可以用来估计吗\ \ beta_j \ ()．有不同的方法可用于估计全局比例因子。有些是针对大量RNA-seq数据开发的，有些是针对scRNA-seq数据开发的。的价值\(\酒吧{β\}\)取决于所使用的协议，可以是批处理相关的。例如，对于基于Droplet的协议，大约是0.06 (inDrop)或0.12 (Drop-seq)。\(\酒吧{β\}\)也可以通过如上所述的spike-ins或smFISH数据进行估计。

我们最后注意到，上面讨论的捕获效率估计将假设细胞具有相似的原始转录内容。因此，bayNorm输出典型细胞的原始转录计数估计值，并根据细胞大小和转录内容的变化进行校正。这通常是下游分析(如DE检测)所需要的。然而，如果对绝对原始计数感兴趣，并且有额外的信息，如细胞大小或每个细胞的总转录内容，则可以为此目的适当地重新调整捕获效率。

4.３先验参数估计

4.3.1边际分配最大化

使用经验贝叶斯方法，可以使用跨细胞观察计数的边际似然分布的最大化来估计先验参数。让\ (M_i \)的边际似然函数\ (^ {th} \ \文本)跨细胞的基因。假设单元间独立，则函数的对数边际分布\ (^ {th} \ \文本)基因是

\[开始\{方程}\ log M_i = \ sum_ {j = 1} ^问\ log \公关(间{ij} | \ mu_i \ phi_i \ beta_j),结束\{方程}\]在哪里\ \公关(间{ij} | \ mu_i \ phi_i \ beta_j) \)为负二项式(见方法)。使这个方程最大化就得到了这个对\ ((\ mu_i \ phi_i) \)．

以上的优化需要对每一个P基因进行。我们参考φ\ (\ \)和/或\μ(\ \)为了方便起见，bayNorm利用了R包中名为“BB”的光谱投影梯度方法(spg)，通过最大化边际分布作为BB估计来估计。优化两者的边际分布\μ(\ \)而且φ\ (\ \)(2D优化)是计算密集型的。如果我们有一个好的估计\μ(\ \)，那么我们就可以优化关于的边际分布φ\ (\ \)单独行动，这样效率更高。

4.4矩量法

一种启发式的估计方法\ \ mu_i \ ()而且\ \ phi_i \ ()是通过矩量法的一种变体。第一步是对原始数据进行简单的规范化，在给定单元特定的捕获效率(\ \ beta_j \ ())．简单规范化计数\(间{ij} ^ \)计算公式如下:

\[开始\{方程}间{ij} ^ s =间{ij} \压裂{\ langle \ sum_ {i = 1} ^ Px_ {ij} / \ beta_j \ rangle_j} {\ sum_ {i = 1} ^ P{间{ij}}},结束\{方程}\]的比例因子的分子在哪里\(间{ij} \)是通过对单元格的缩放总数取平均值得到的。

基于简单的标准化数据，我们能够估计先验参数\μ(\ \)而且φ\ (\ \)使用矩估计方法(MME)，它简单地等同于理论和经验矩的负二项分布。这种估计方法是快速的，仿真表明它提供了很好的估计\μ(\ \)但缺点是估计φ\ (\ \)显示系统偏倚(见附录图S24 a-b)1)．

4．5组合方法(bayNorm中用于估计先验的默认设置)

基于模拟研究(补充图S24)1)，最稳健和有效的估计\μ(\ \)而且φ\ (\ \)可以通过以下组合方式获得，这是bayNorm中的默认设置:

  1. 在简单归一化数据的基础上，对每个基因使用MME方法得到MME估计值\μ(\ \)而且φ\ (\ \)．
  2. 尽管BB估计φ\ (\ \)更接近真实吗φ\ (\ \)时，许多估值位于搜索空间的上边界(补充图S24 c-d in1)．因此，我们找到了调整估计的MMEφ\ (\ \)通过拟合估计的MME之间的线性回归来估计一个因子φ\ (\ \)和BB估计φ\ (\ \)效果最佳(补充图S24 c-d1)．这是调整后的MME估计值φ\ (\ \)连同MME估算\μ(\ \)估计\ \ beta_j \ ()可用于近似每个细胞中每个基因的后验分布。

根据细胞特定的属性(\ (C_j \))．无论实验条件如何，都可以对所有细胞进行先验估计。我们把这个过程称为“全局的”。或者，假设数据集中有多组细胞，我们有理由相信每组细胞的行为可能不同。然后我们可以估计先验参数"\μ(\ \)而且φ\ (\ \)，在每组内分别(组内具有相同的\ (C_j \)值)。我们将此过程称为“本地”过程。基于“全局”过程估计特定单元组的先验参数可以消除潜在的批处理效应。基于“局部”程序的多组归一化可以放大组间的差异，同时减轻组内的可变性，适用于DE检测。

5会话信息

sessionInfo ()
## R版本4.2.0 RC (2022-04-19 r82224) ##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 20.04.4 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas。/home/biocbuild/bbs-3.15-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]stats4 stats graphics grDevices utils datasets methods ##[8]基础## ##其他附加包:[1] SingleCellExperiment_1.18.0 SummarizedExperiment_1.26.0 ## [3] Biobase_2.56.0 genomicranges_1 .32.0 IRanges_2.30.0 ## [7] S4Vectors_0.34.0 BiocGenerics_0.42.0 ## [9] MatrixGenerics_1.8.0 matrixStats_0.62.0 ## [11] bayNorm_1.14.0 knitr_1.38 ## [13] BiocStyle_2.24.0 ## ##通过命名空间加载(并且没有附加):## [28] htmltools_0.5.2 snow_0.4-4 yaml_2.3.5 ## [10] survival_3.3-1 rlang_1.0.2 jquerylib_0.1.4 ## [13] fitdistrplus_1.1-8 BiocParallel_1.30.0 GenomeInfoDbData_1.2.8 ## [19] codetools_0.2-18 evaluate_0.15 fastmap_1.1.0 ## [22] parallel_4.2.0 Rcpp_1.0.8.3 BiocManager_1.30.17 ## [25] DelayedArray_0.22.0 jsonlite_1.8.0 XVector_0.36.0 ## [28][34] bitops_1.0-7 magrittr_2.0.3 sass_0.4.1 ## [37] RCurl_1.98-1.6 MASS_7.3-57 Matrix_1.4-1 ## [40] rmarkdown_2.14 iterators_1.0.14 R6_2.5.1 ## [43] compiler_4.2.0

6参考文献