##装载所需的包:针织机包:基础知识
作者:Catalina Vallejos (cnvallej@uc.cl)及艾玲(eling@ebi.ac.uk）
编辑日期:2018-03-18
单细胞mRNA测序可以揭示看似同质的细胞群体中基因表达水平的新细胞间异质性。然而，这些实验容易出现高水平的技术噪音，为识别在研究细胞组中显示真正异质表达的基因带来了新的挑战。
基础知识(Bayesian一个分析的如果角-Cl形的年代排序数据)是一个集成的贝叶斯层次模型，它通过同时执行数据规范化(全局缩放)、技术噪声量化和两种类型的统计不确定性来传播统计不确定性监督下游分析:
对于一组细胞[1]: BASiCS提供了一个识别组内高(和低)变量基因的标准。
用于两组(或更多组)细胞[2]: BASiCS允许识别组间差异表达基因。与传统的差分表达式工具一样，BASiCS可以揭示组间平均表达式的变化。除此之外，BASiCS还可以发现over-dispersion-在考虑技术噪声后观察到的剩余细胞间变异的度量。例如，这一特征导致了在衰老过程中免疫细胞的新见解。
在这两种情况下，都提供了一个概率输出，通过预期错误发现率(EFDR)[4]校准后验概率阈值。
目前，BASiCS依赖于的使用激增的基因-人工引入到每个细胞的裂解液中-来进行这些分析。
本文档的“方法论”部分提供了在BASiCS中实现的统计模型的简要描述。
重要的BASiCS是在监督实验的背景下设计的，其中所研究的细胞组(例如实验条件，细胞类型)是先验的(例如病例对照研究)。因此，我们不建议在无监督的环境中使用BASiCS，其目的是通过聚类来揭示细胞的亚群。
basic的输入数据集必须存储为SingleCellExperiment对象(如SingleCellExperiment包)。
的newBASiCS_Data函数可以根据以下信息创建输入数据对象:
计数:原始表达式的矩阵以\(q\)乘以\(n\)维数计数。在这个矩阵中，\(q_0\)行必须对应生物基因，\(q-q_0\)行必须对应技术尖峰基因。基因名称必须存储为rownames(计数)．
科技的向量真正的/假长度为q的元素。如果Tech[i] = FALSE的基因我是生物;否则基因就会被插入。此载体必须以与中相同的基因顺序指定计数矩阵。
SpikeInfo:一个data.frame有\(q-q_0\)行。第一列必须包含与尖峰基因相关的名称(如rownames(计数))．第二列必须包含插入基因的输入分子数量(每个细胞的数量)。
BatchInfo(可选):长度\(n\)向量，表示批处理结构(当单元格使用多批处理时)。
例如，下面的代码模拟了一个包含50个基因(40个生物基因和10个插入基因)和40个细胞的数据集。
set.seed(1) Counts = Counts = matrix(rpois(50* 40,2)， ncol = 40) rownames(Counts) <- c(paste0("Gene"， 1:40)， paste0("Spike"， 1:10)) Tech = c(rep(FALSE,40)，rep(TRUE,10)) set.seed(2) SpikeInput = rgamma(10,1,1) SpikeInfo <- data.frame("SpikeID" = paste0("Spike"， 1:10)， "SpikeInput" = SpikeInput) #无批处理结构DataExample = newBASiCS_Data(Counts, Tech, SpikeInfo) #有批处理结构DataExample = newBASiCS_Data(Counts, Tech, SpikeInfo)， BatchInfo = rep(c(1,2)， each = 20))注:scRNA-seq数据集在执行分析之前通常需要质量控制过滤。这是为了去除表达量非常低的细胞和/或转录本。这个函数BASiCS_Filter可用于执行此任务。有关示例，请参见帮助(BASiCS_Filter)．
注意:嘘软件包为scRNA-seq数据集的预处理提供了增强的功能。
要转换现有的SingleCellExperiment对象(数据)转换为一个可以在BASiCS中使用的对象，元信息必须存储在对象中。
isSpike(数据、“ERCC”)=科技:指示生物/技术基因(见上文)的逻辑载体必须储存在int_metadata槽通过isSpike函数。
元数据(数据):SpikeInfo而且BatchInfo对象存储在元数据槽位号SingleCellExperiment对象:metadata(Data)=list(SpikeInput = SpikeInfo[，2]， BatchInfo = BatchInfo)．一旦包含了附加信息，就可以在BASiCS中使用对象。
参数估计使用BASiCS_MCMC函数。运行该算法的基本参数是:
N:总迭代次数薄:减薄周期的长度(即仅每薄的输出中存储迭代BASiCS_MCMC）燃烧:老化期长度(即初始燃烧类的输出中将丢弃的迭代BASiCS_MCMC）如果可选参数PrintProgress设置为真正的， R控制台将显示MCMC算法的进度。其他可选参数请参见帮助(BASiCS_MCMC)．
在这里，我们说明的用法BASiCS_MCMC使用内置的合成数据集。
Data <- makeExampleBASiCS_Data() Chain <- BASiCS_MCMC(Data = Data, N = 1000, Thin = 10, Burn = 500, PrintProgress = FALSE)## ------------------------------------------------------------- ## 密度采样器已经开始:1000次迭代。## ------------------------------------------------------------- ## ------------------------------------------------------------- ## 老化期结束。## ------------------------------------------------------------- ## ## ------------------------------------------------------------- ## ------------------------------------------------------------- ## 所有1000名获得迭代已经完成。## ------------------------------------------------------------- ## ------------------------------------------------------------- ## ## ------------------------------------------------------------- ## 请参阅下面的总结整体录取率。## ------------------------------------------------------------- ## ## mu[i] s中的最小接受率:0.476 ## mu[i] s中的平均接受率:0.70152 ## mu[i] s中的最大接受率:0.896 ## ## delta[i] s中的最小接受率:0.506 ## delta[i] s中的平均接受率:0.58424 ## delta[i] s中的最大接受率:0.706 ## ## phi(关节)的接受率:0.924 ## ## nu中的最小接受率:0.462 # #之间的平均录取率ν[j]的年代:0.5504 # #最大ν[j]的接受率:0.71 # # # #θ最低接受率[k]的年代:0.792 # #平均录取率在θ[k]的年代:0.792 # #θ[k]最大的接受率的年代:0.792 ## ## ------------------------------------------------------------- ##重要讲话:
的控制台输出中显示的接受率BASiCS_MCMC在0.44左右。如果它们离这个值太远，就应该增加N而且燃烧．
它是至关重要的来评估MCMC算法的收敛性之前执行下游分析。有关此步骤的指导，请参阅本插图的“收敛评估”部分
通常,设置N = 20000，瘦= 20而且燃烧= 10000导致稳定的结果。
我们在分析[5]中提供的单细胞样本时，使用[2]中获得的MCMC链的一小部分提取物来说明这一分析。这包括在基础知识随着ChainSC数据集。
数据(ChainSC)下面的代码是用来识别高变基因(HVG)而且低变基因在这些细胞内。的VarThreshold参数为分配给生物成分的可变性比例设置较低的阈值(σ)．在下面的例子中:
对于每个基因，这些函数返回后验概率作为HVG/LVG证据的度量。通过控制EFDR(默认选项:EviThreshold定义使EFDR = 0.10)。
par(mfrow = c(2,2)) HVG <- BASiCS_DetectHVG(ChainSC, VarThreshold = 0.6, Plot = TRUE) LVG <- BASiCS_DetectLVG(ChainSC, VarThreshold = 0.2, Plot = TRUE)
要访问这些测试的结果，请使用。
头(美元HVG表)21 21 Map3k11 10.044337 2.344024 0.7599326 1.00 TRUE ## 71 71 Lefty1 3.574049 2.972835 0.7657068 0.96 TRUE ## 48 48 Ctgf 4.296390 2.455164 0.7370870 0.91 TRUE ## 88 88 Naa11 4.997926 2.503697 0.7438627 0.91 TRUE ## 166 166 Pnma2 2.773063 2.781834 0.7270004 0.89 TRUE ## 185 185 Alg8 3.931154 2.489690 0.7379904 0.88 TRUE头美元举债经营(LVG表)## GeneIndex GeneName Mu Delta Sigma Prob LVG ## 16 16 Gm10653 1165.3430 0.04023528 0.05845448 1 TRUE ## 63 63 Luc7l2 1942.1656 0.07219708 0.10128552 1 TRUE ## 65 65 Atp5g2 338.0637 0.06562415 0.09368431 1 TRUE ## 89 89 Rpl14 1348.3208 0.02466278 0.03720835 1 TRUE ## 90 Rpl11 1256.5145 0.02305986 0.03439021 1 TRUE ## 92 92 Rcc2 294.2660 0.06411164 0.09191962 1 TRUESummarySC <- Summary(ChainSC) plot(SummarySC, Param = "mu"， Param2 = "delta"， log = "xy") with(HVG$Table[HVG$Table$HVG == TRUE，]， points(mu, delta)) with(LVG$Table[LVG$Table$LVG == TRUE，]， points(mu, delta))
请注意:相对于BASiCS的原始版本，这个阈值标准已经发生了变化EviThreshold定义为EFDR = EFNR)。然而，新的选择更稳定(有时，不可能找到一个阈值，使EFDR = EFNR)。
为了说明在两个细胞群之间使用差分平均表达和差分过分散检验，我们在分析[5]数据集(单细胞vs池-分裂样本)时使用了从[2]中获得的MCMC链的提取。这些都是通过独立运行BASiCS_MCMC每组单元格的函数。
数据(ChainSC)数据(ChainRNA)Test <- BASiCS_TestDE(Chain1 = ChainSC, Chain2 = ChainRNA, GroupLabel1 = "SC"， GroupLabel2 = "PaS"， EpsilonM = log2(1.5)， EpsilonD = log2(1.5)， EFDR_M = 0.10, EFDR_D = 0.10, Offset = TRUE, OffsetPlot = TRUE, Plot = TRUE)
在BASiCS_TestDE，EpsilonM在表达式(\(\mu\))和中设置log2 fold变化(log2FC)EpsilonDlog2FC过色散(\(\delta\))。作为默认选项:ε = ε = log2(1.5)(即增加50%)。为了调整总体RNA含量的差异，当抵消= TRUE．这是推荐的默认值。
生成的输出列表可以使用
头(测试TableMean美元)# # GeneName MeanOverall Mean1非常刻薄的MeanFC MeanLog2FC ProbDiffMean # # 47 BC018473 1502.198 2948.429 94.024 31.251 4.966 1 # # 71 Lefty1 10.756 5.170 16.195 0.328 -1.607 1 # # 329二488.756 354.929 619.062 0.573 -0.802 1 # # 418 Zfp937 149.071 221.552 78.498 2.808 1.490 1 # # 437 Snora33 6.428 1.943 10.795 0.176 -2.510 1 # # 445 Zfp71-rs1 126.059 62.323 188.118 0.329 -1.605 1 # # ResultDiffMean # # 47 SC + # # 71不是+ # # 329不是+ # # 418 SC + # # 437不是445不是+ + # #头(测试TableDisp美元)# # GeneName MeanOverall DispOverall Disp1 Disp2 DispFC DispLog2FC # # 49 Gm5643 1821.217 0.062 0.105 0.020 5.171 - 2.370 # # 58 Luc7l2 2809.375 0.045 0.072 0.018 4.076 - 2.027 # # 87 Hist2h2ab 70.985 0.346 0.591 0.107 5.626 2.492 # # 94 Gm16287 658.777 0.081 0.141 0.022 6.669 - 2.737 # # 102 Zyg11b 1657.344 0.063 0.108 0.019 5.772 2.529 # # 117 Stxbp2 615.573 0.064 0.105 0.024 4.349 - 2.121 # # ProbDiffDisp ResultDiffDisp # # 49 1 SC + # # 58 1 SC + # # 87 SC + # # 94 SC + # # 102 1 SC SC + + # # 117 1注意:由于通常在scRNA-seq数据集中观察到的均值和过度分散之间的混淆，我们只评估那些组间均值没有变化的基因的过度分散变化。使用= 0作为保守的选择:
Test <- BASiCS_TestDE(Chain1 = ChainSC, Chain2 = ChainRNA, GroupLabel1 = "SC"， GroupLabel2 = "PaS"， EpsilonM = 0, EpsilonD = log2(1.5)， EFDR_M = 0.10, EFDR_D = 0.10, Offset = TRUE, OffsetPlot = TRUE, Plot = TRUE)
的输出在外部存储BASiCS_MCMC(推荐)，附加参数StoreChains，StoreDir而且RunName是必需的。例如:
Data <- makeExampleBASiCS_Data() Chain <- BASiCS_MCMC(Data, N = 1000, Thin = 10, Burn = 500, PrintProgress = FALSE, StoreChains = TRUE, StoreDir = tempdir()， RunName = "Example")## ------------------------------------------------------------- ## 密度采样器已经开始:1000次迭代。## ------------------------------------------------------------- ## ------------------------------------------------------------- ## 老化期结束。## ------------------------------------------------------------- ## ## ------------------------------------------------------------- ## ------------------------------------------------------------- ## 所有1000名获得迭代已经完成。## ------------------------------------------------------------- ## ------------------------------------------------------------- ## ## ------------------------------------------------------------- ## 请参阅下面的总结整体录取率。## ------------------------------------------------------------- ## ## mu[i] s中的最小接受率:0.476 ## mu[i] s中的平均接受率:0.70152 ## mu[i] s中的最大接受率:0.896 ## ## delta[i] s中的最小接受率:0.506 ## delta[i] s中的平均接受率:0.58424 ## delta[i] s中的最大接受率:0.706 ## ## phi(关节)的接受率:0.924 ## ## nu中的最小接受率:0.462 # #之间的平均录取率ν[j]的年代:0.5504 # #最大ν[j]的接受率:0.71 # # # #θ最低接受率[k]的年代:0.792 # #平均录取率在θ[k]的年代:0.792 # #θ[k]最大的接受率的年代:0.792 ## ## ------------------------------------------------------------- ##在本例中，的输出BASiCS_MCMC将被存储为BASiCS_Chain“chain_Example”文件中的对象。Rds”，在tempdir ()目录中。
为了加载预先计算的MCMC链，
Chain <- BASiCS_LoadChain("Example"， StoreDir = tempdir())要评估链的收敛性，包提供的收敛诊断coda可以使用。此外，可以直观地检查链条。例如:
plot(Chain, Param = "mu"， Gene = 1, log = "y")
plot(Chain, Param = "phi"， Cell = 1)
在上图中:
acf ?）要访问与单个参数关联的MCMC链，请使用该函数displayChainBASiCS．例如,
displayChainBASiCS(Chain, Param = "mu")[1:5,1:5]## Gene1 Gene2 Gene3 Gene4 Gene5 ## [1，] 8.516936 4.502855 3.826309 5.657261 19.23518 ## [2，] 5.548291 3.961501 2.843601 3.589559 19.51600 ## [3，] 10.131468 5.688056 5.369110 4.545069 19.64434 ## [4，] 5.664402 5.120068 4.006658 6.431824 22.93474 ## [5，] 7.068831 4.759326 3.862751 7.077321 26.41031作为后验分布的总结，函数总结计算每个模型参数的后验中位数和高后验密度(HPD)区间。作为默认选项，HPD间隔包含0.95概率。
ChainSummary <- Summary(Chain)这个函数displaySummaryBASiCS提取个体参数的后验总结。例如
(displaySummaryBASiCS(ChainSummary, Param = "mu"))## Mu lower upper ## Gene1 7.694477 4.331570 10.680111 ## Gene2 4.977089 2.588200 7.964820 ## Gene3 4.109454 2.509790 6.490364 ## Gene4 4.843340 2.742174 8.544155 ## Gene5 18.517451 10.747134 23.065117 ## Gene6 8.766424 5.511452 13.280344下图显示了基因特定参数\(\mu_i\)(均值)和\(\delta_i\)(过分散)的后验中值和相应的HPD 95%区间
par(mfrow = c(2,2)) plot(ChainSummary, Param = "mu"， main = "所有基因"，log = "y") plot(ChainSummary, Param = "mu"，基因= 1:10,main = "前10个基因")plot(ChainSummary, Param = "delta"， main = "所有基因")plot(ChainSummary, Param = "delta"，基因= c(2,5,10,50)， main = "5个定制基因")
也可以得到归一化常数\(\phi_j\)和\(s_j\)的类似摘要。
图(ChainSummary, Param = "s"， Cells = 1:5)
最后，也可以为基因特异性参数的后验估计创建散点图。通常情况下，该图将显示在平均和过度分散之间观察到的混淆效应。
par(mfrow = c(1,2)) plot(ChainSummary, Param = "mu"， Param2 = "delta"， log = "x"， SmoothPlot = FALSE)
的选项SmoothPlot = TRUE通常推荐，因为在分析实际数据集时，该图将包含数千个基因。
也有可能产生一个矩阵归一化和去噪的表达式计数，其中技术变化的影响被删除。为此，我们实现了这个函数BASiCS_DenoisedCounts．对于每个基因\(i\)和细胞\(j\)，该函数返回
\ [x ^ * _ {ij} = \压裂{间{ij}}{{φ\}\帽子帽子_j \{\ν}_j}, \]
其中\(x_{ij}\)是观察到的基因\(i\)在细胞\(j\)中的表达量，\(\hat{\phi}_j\)表示\(\phi_j\)的后中值，\(\hat{\nu}_j\)是\(\nu_j\)的后中值。
DenoisedCounts = BASiCS_DenoisedCounts(Data = Data, Chain = Chain) DenoisedCounts[1:5, 1:5]##[，1][，2][，3][，4][，5] ##基因1 0.000000 34.355728 47.81306 6.127459 35.141310 ##基因2 0.000000 0.000000 0.00000 12.254919 0.000000 ##基因3 0.000000 4.580764 0.00000 6.127459 0.000000 ##基因4 4.915316 2.290382 23.90653 0.000000 25.770294 ##基因5 4.915316 6.871146 0.00000 55.147134 9.371016或者，用户可以计算归一化和去噪的表达率，这些表达率是跨细胞的所有基因表达的基础BASiCS_DenoisedRates．这个函数的输出由
\[\Lambda_{ij} = \hat{\mu_i} \hat{\rho}_{ij}， \]
其中\(\hat{\mu_i}\)表示\(\mu_j\)的后置中值，\(\hat{\rho}_{ij}\)由其后置均值(基于所有模型参数的MCMC样本的蒙特卡罗估计)给出。
降噪率<- basics_降噪率(数据=数据，链=链，倾向= FALSE)降噪率[1:5,1:5]## [，1] [，2] [，3] [，4] [，5] ## [1，] 3.069806 30.5252453 18.407656 7.091326 30.9276224 ## [2，] 1.656280 0.9535089 3.420656 9.734775 0.9691411 ## [3，] 2.481371 4.4320480 3.584329 4.830411 1.7525373 ## [4，] 4.864397 2.9188147 9.113103 2.222938 20.4012677 ## [5，] 9.262084 9.1816807 12.2837785 43.926025 11.3747969另外，也可以提取去噪的表达倾向\(\hat{\rho}_{ij}\)
DenoisedProp = BASiCS_DenoisedRates(数据=数据，链=链，倾向= TRUE) DenoisedProp[1:5, 1:5]## [，1] [，2] [，3] [，4] [，5] ## [1，] 0.3989623 3.9671630 2.3923205 0.9216124 4.0194573 ## [2，] 0.3327809 0.1915797 0.6872806 1.9559176 0.1947205 ## [3，] 0.6038202 1.0785005 0.8722154 1.1754387 0.4264648 ## [4，] 1.0043476 0.6026450 1.8815741 0.4589680 4.2122313 ## [5，] 0.5001814 0.4958393 0.6633626 2.3721421 0.6142744我们首先描述了[1]中介绍的模型，它与一组细胞有关。
在整个过程中，我们考虑\(q\)基因的表达计数，其中\(q_0\)在被研究的细胞群体中表达(生物基因)，其余的\(q-q_0\)是外部的插入(技术)基因。设\(X_{ij}\)为一个随机变量，表示一个基因\(i\)在细胞\(j\)中的表达量(\(i=1，\ldots,q\);\ (j = 1, \ ldots n \))。BASiCS基于以下层次模型:\[X_{ij} \大| \mu_i， \phi_j， \nu_j， \rho_{ij} \sim \left\{\begin{array}{ll} \mbox{Poisson}(\phi_j \nu_j \mu_i \ ij})， \mbox{for}i=1，\ldots,q_0, j=1，\ldots,n \mbox{for}i=q_0+1，\ldots,q, j=1，\ldots,n， \end{array} \right.\]
其中\(\nu_j\)和\(\rho_{ij}\)是相互独立的随机效应，使得\(\nu_j|s_j，\theta \sim \mbox{Gamma}(1/\theta,1/ (s_j \theta))\和\(\rho_{ij} | \delta_i \sim \mbox{Gamma}(1/\ delta_i,1/\delta_i)\)[footnoteGamma]。
该模型的图形表示如下所示。这是基于2个基因(\(i\):生物和\(i'\):技术)在2个细胞(\(j\)和\(j'\))的表达计数。方形和圆形节点分别表示已知的观察量(观察到的表达计数和添加的spike-in mRNA分子数量)和未知元素。黑色圆形节点表示在我们的层次结构中起中间作用的随机效应，红色圆形节点与我们模型中层次结构顶层的未知模型参数有关。蓝色、绿色和灰色区域分别突出了生物基因、技术基因或细胞中共享的元素。
\中心线{\ includegraphics(高度= 4){。/ BASiCS_DAG.jpg}}
在此设置中，用于下游分析的关键参数是:
\(\mu_i\):基因\(i\)在被研究细胞组中的平均表达参数。对于插入技术基因，假设\(\mu_i\)是已知的，并且等于相应插入技术基因的输入分子数)。
\(\delta_i\):基因\(i\)的过色散参数，在考虑技术噪声(关于泊松采样)后观察到的过度变异的度量
附加(有害)参数解释如下:
\(\phi_j\):与mRNA含量差异相关的细胞特异性归一化参数(可识别性约束:\(\sum_{j=1}^n \phi_j = n\))。
\(s_j\):与技术单元格特定偏差相关的单元格特定归一化参数(有关此解释的更多详细信息，请参阅[6])。
\(\theta\):技术过色散参数，控制细胞间技术变异性的强度。
当来自同一组的细胞在多个测序批次中处理时，该模型被扩展，因此技术上的过分散参数\(\theta\)是批次特定的。该扩展允许为每批单元推断不同强度的技术噪声。
［footnoteGamma]:我们将Gamma分布参数化，如果\(X \sim \mbox{Gamma}(a,b)\)，则\(\mbox{E}(X)=a/b\)和\(\mbox{var}(X)=a/b^2\)。
在[2]中，该模型已扩展到研究多组细胞的情况。这是通过假设基因特异性参数也是群体特异性参数来实现的。基于此设置，差异表达的证据通过组间基因特异性参数(均值和过度分散)的log2倍变化进行量化。
关于模型设置、先前的规范和实现的更多细节在[1]和[2]中描述。
我们感谢马里奥尼实验室(EMBL-EBI;CRUK-CI)在整个R库开发过程中提供支持和讨论。我们特别感谢Aaron Lun (@LTLA, CRUK-CI)在准备Bioconductor提交期间提供的建议和支持。
我们也感谢来自(括号内提供的Github别名)的反馈和贡献:Ben Dulken (@bdulken)， Chang Xu (@xuchang116)， Danilo Horta (@Horta)， Dmitriy Zhukov (@dvzhukov)， Jens Preußner (@jenzopr)， Joanna Dreux (@Joannacodes)， Kevin Rue-Albrecht (@kevinrue)， Luke Zappia (@lazappi)， Simon Anders (@s-andrews)， Yongchao Ge和Yuan Cao (@yuancao90)，等等。
这项工作得到了MRC生物统计部门的资助(MRC批准号:MRC_MC_UP_0801/1;Catalina Vallejos和Sylvia Richardson)， EMBL欧洲生物信息学研究所(核心欧洲分子生物学实验室资助;Catalina Vallejos, Nils Eling和John Marioni)， CRUK剑桥研究所(CRUK核心资金;约翰·马里奥尼)和艾伦·图灵研究所(EPSRC批准no。EP / N510129/1;卡特琳娜Vallejos)。
[1] Vallejos CA, Marioni JCM和Richardson S(2015)基础:单细胞测序数据的贝叶斯分析。PLoS计算生物学11 (6)， e1004333。
[2] Vallejos CA, Richardson S和Marioni JCM(2016)超越手段的比较:了解单细胞水平上基因表达的变化。基因组生物学17(1)， 1-14。
[3] Martinez-Jimenez CP, Eling N, Chen H, Vallejos CA, Kolodziejczyk AA, Connor F, Stojic L, Rayner TF, Stubbington MJT, Teichmann SA, de la Roche M, Marioni JC和Odom DT(2017)免疫刺激后衰老增加细胞间转录变性。科学355(6332)， 1433-1436。
[4]牛顿MA, Noueiry A, Sarkar D, Ahlquist P(2004)用半参数分层混合方法检测差异基因表达。生物统计学5(2)， 155-76。
[5] Grün D, Kester L, van Oudenaarden A(2014)单细胞转录组学噪声模型的验证。自然方法11(6)， 637-40。
[6] Vallejos CA, Risso D, Scialdone A, Dudoit S和Marioni JCM(2017)规范化单细胞rna测序数据:挑战和机遇。自然方法14日,565 - 571。
[7]罗伯茨GO和罗森塔尔JS(2009)。自适应MCMC的例子。计算与图形统计杂志18: 349 - 367。
sessionInfo ()## R版本3.4.3(2017-11-30)##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 16.04.4 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.6-bioc/R/lib/libRblas。/home/biocbuild/bbs-3.6-bioc/R/lib/libRlapack。所以## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_US。UTF-8 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统计图形grDevices utils数据集##[8]方法基础## ##其他附加包:## [1] BASiCS_1.0.1 SingleCellExperiment_1.0.0 ## [3] SummarizedExperiment_1.8.1 DelayedArray_0.4.1 ## [5] matrixStats_0.53.1 Biobase_2.38.0 ## [7] GenomicRanges_1.30.3 GenomeInfoDb_1.14.0 ## [9] IRanges_2.12.0 S4Vectors_0.16.0 ## [11] BiocGenerics_0.24.0 BiocStyle_2.6.1 ## [13] knitr_1.20 ## ##通过命名空间加载(并且没有附加):[1] bitops_1.0-6 bit64_0.9-7 progress_1.1.2 ## [4] httr_1.3.1 rprojroot_1.3-2 dynamicTreeCut_1.63-1 ## [7] tools_3.4.3 backports_1.1.2 DT_0.4 ## [10] R6_2.2.2 KernSmooth_2.23-15 vipor_0.4.5 ## [13] DBI_0.8 lazyeval_0.2.1 colorspace_1.3-2 ## [16] gridExtra_2.3 prettyunits_1.0.2 bit_1.1-12 ## [19] compiler_3.4.3 scales_0.5.0 string_1 .3.0 ## [22] digest_0.6.15 rmarkdown_1.9 XVector_0.18.0 ## [25] scater_1.6.3 pkgconfig_2.0.1 htmltools_0.3.6 ## [28] highr_0.6 limma_3.34.9 ## [31]rlang_0.2.0 RSQLite_2.0 FNN_1.1 ## [34] shiny_1.0.5 bindr_0.1.1 zoo_1.8-1 ## [37] BiocParallel_1.12.0 dplyr_0.7.4 RCurl_1.95-4.10 ## [40] magrittr_1.5 GenomeInfoDbData_1.0.0 Matrix_1.2-12 ## [43] Rcpp_0.12.16 ggbeeswarm_0.6.0 munsell_0.4.3 ## [46] viridis_0.5.0 stringi_1.1.7 yaml_2.1.18 ## [49] edgeR_3.20.9 zlibbioc_1.24.0 rhdf5_2.22.0 ## [52] plyr_1.8.4 grid_3.4.3 blob_1.1.0 ## [55] shinydashboard_0.6.1 lattice_0.20-35 locfit_1.5-9.1 ## [58] pillar_1.2.1 igraph_1.2.1 rjson_0.2.15 ## [61] reshape2_1.4.3 biomaRt_2.34.2 XML_3.98-1.10 ## [64] glue_1.2.0 evaluate_0.10.1 scran_1.6.9 ## [67] data.table_1.10.4-3 httpuv_1.3.6.2 testthat_2.0.0 ## [70] gtable_0.2.0 assertthat_0.2.0 ggplot2_2.2.1 ## [73] mime_0.5 xtable_1.8-2 coda_0.19-1 ## [76] viridisLite_0.3.0 tibble_1.4.2 AnnotationDbi_1.40.0 ## [79] beeswarm_0.2.3 memoise_1.1.0 tximport_1.6.0 ## [82] bindrcpp_0.2 statmod_1.4.30