从 scVI 到 Total VI（二）| 单细胞转录组中的变分推断

我们在上一篇推送中介绍了scVI的模型基础和设计思路，先回顾一下。

一般的神经网络的整个前向和后向传播过程都是确定性的，所以之前实在很难想象神经网络是如何和概率分布联系在一起的。

• 自从Kingma的VAE提出以后， reparameterization 就搭建了借用加noise做sampling并且还能允许后向传播的的神经网络构建分布。
• 如果不做sampling，而是只是推断概率密度函数的参数，则又回到了确定性问题。
• 如果是这样，那么对于形式确定的分布，我们就可以显式地写出概率密度 ，进而对整个dataset 计算似然 。

scVI 的前半部分通过四个 decoder分别编码测序深度 和隐变量 的均值和标准差。

scVI 后半部分的两个 decoder 就是在做ZINB分布的参数推断， 第一个参数 是NB的均值， 第二个参数是 ， 第三个表示dispersion的参数 则直接指定为可更新的参数， 在pytorch 里就可以通过实例化nn.Parameter对象并设置requires_grad=True来实现。

这一篇我们来看看我们可以用scVI来做些什么？

scVI 隐变量有大用处

👆图中其实可以看到，在模型的隐变量 🥷 和输出 📤 各有一个箭头指向了它们的下游应用：

• 隐变量 🥷 的作用：聚类、 可视化、 batch 矫正
• 输出 📤 的作用： 数据补全 、 差异表达分析

有了数据， 定义了模型。 fit完之后得到一个从 到 的东西，自己猜自己？到底有啥用？这个在【链接】 中提过。 隐变量模型的核心在于得到可以概括数据的低维隐变量 。

batch correction 和 聚类

scVI 的 是不包含 batch 信息 的。 大家可以再仔细看一下模型图， 经过编码器之后得到那四个参数。 但是在生成的过程中， 又重新放在  的边上一起输入解码器。这么做就能让编码器在编码的时候意识到batch的不同，然后把 s 从表达矩阵 中解耦。

因为隐变量浓缩了原始表达数据的规律。 所以把低维数据 来代表一个细胞， 我们可以发现用它构建的UMAP和用表达矩阵构建的UMAP 很相似。 cell type的信息还是在隐空间中保留了下来（A）。而热图中细胞clustering的结果也和SIMLR的结果相差不大， 主要的大类都被保留下来了。 最后是第三个数据集 RETINA 中，很明显地SIMLR没有办法融合两个batch。 而scVI的两个batch 在UMAP中融合地很均匀，说明 中尽量保存了生物学信息而去掉了batch的信息 。

基于隐变量的可视化

可视化这个过程如何应用隐变量不像上面两个功能那么直观。其实基于隐变量的这三个都是紧密结合在一起的。UMAP或者t-SNE这两个可视化降维是基于K -nearest neighbor 的，那么用什么数据来构建 KNN 就决定了最后单细胞在图中的散布。

基于原始的表达矩阵（raw count）计算KNN是单细胞分析的基本做法。但是raw count收到许多实验变量的影响，用它计算出来的KNN会出现保留下这些变量带来的偏差。scVI 的隐变量 得益于其 解耦联机制，能把batch的信息给提取出来，然后 只保留其生物学特性，，因此用低维隐变量 构建的KNN可以既覆盖原来的基因表达特点又可以去掉这些实验变量的影响。之后再把这个graph 作为输入传入UMAP和t-SNE就能得到一个更加纯净的可视化结果了。

基于scVI的后验估计

校正观测数据

隐变量很有用🤓， 它能以更低的维概括高维的数据， 它的某个维度甚至可能含有具体的生物学含义。那scVI的高维输出有什么用呢🤔？按理来说它只是重建模型的输入而已， 考虑到重建还有误差， 我为什么不直接看输入呢 ❓

回到最基础的AE，模型的输出和输入的维度是相同的， 基本上就是重建输入。但是到了VAE， 由于重参数化这个过程引入了不确定性，为了让输出和输入更像， decoder只能想办法更贴近数据的分布 。最后到了scVI， 我们不仅仅满足于让隐变量 满足某种分布，我们还想让 满足某种形式（例如上面提到的ZINB），因此不仅增加了不确定， 而且模型的输出 📤 成为有统计意义的参数。因为scVI的 输出实际上在近似表达矩阵的后验 ，所以我们称之为细胞n和基因 的后验估计。

scVI 里对 的后验估计。注意这里还以 和 为底做了积分。 这个过程可以直接理解为对多次的输出求个平均，因为re-parameterization有不确定性所以单独一次的输出其实也是不准的。

👆说了一大堆， 还是没有说我们为什么需要后验估计。

对每一个观测到的细胞，基因A的reads count有可能是和其他处于同一个cell state的细胞的同一个基因的reads count相差很远。其误差可以来源于细胞之间的生物学差异， 也可能是测序过程引入的技术差异。Bulk-seq的数据由于把所有细胞合在一起，会平均化这种差异。 scRNA-seq则会把这个variation暴露出来。如果该细胞的其他基因的表达量和同cell state的其他细胞的表达量差不多，但是唯独基因A的表达差很多。那我们大概会认为这个基因的read count是不准的。scVI 的对该基因A的后验估计由其隐变量决定，而隐变量基于所有的基因计算的。因此我们可以认为scVI的输出📤的每个维度都是考虑了整体表达之后对输入的观测数据的校正。

用基因表达举证作为例子还是太抽象了（狗头）， 以柴犬为例吧。假设你现在有一堆柴犬的图片， 你通过scVI也学了柴犬图片的概率分布。 这时你遇到一哪哪都像柴犬， 唯独头是大橘的图片。 那么经过scVI 之后， 模型会返回以一个更符合柴犬分布的图像。

填补缺失值

scRNA-seq的数据是高度稀疏的，这也是常用使用ZINB来描述其raw count matrix 的原因。表达量较低的基因有更大的可能会因为drop-out effect而观测不到任何的reads。这个时候我们就可以用上面提到的方法用后验估计来填补这些drop-out的count。

Supplementary Figure 4-5 记录了scVI用于在CORTEX dataset上填补缺失值的表现。

我们可以看到一共三个baseline model拿来作为比较。 ZIFA和ZINB-WaVE也都属于用ZINB来模拟RNA expression的模型。 最后一个模型MAGIC的表现就显得那么好。 右图是binomial corrupted 的dataset。

至于模型输出的第二个应用—— 用scVI进行差异表达分析， 我们就放到下一次再分享啦。

TotalVI ：joint-representation的scVI

**scVI**就是个框架， 它仅用简单的设计就很好地建模了scRNA-seq数据，它还兼容了很多对表达数据统计分析的方法。在框架搭建好之后，我们将它拓展到新的单细胞数据类型：

现在，我们拓展到CITE-seq数据上。

TotalVI的模型构架如👇所示， 出乎我意料的是， 它并没有将X（RNA 表达）和Y（蛋白水平）分开并单独用自己的encoder。好处是显而易见的，那就是简单且只有一个隐变量。利用这个联合隐空间我们将处于两个omics-type空间的raw count用统一且连续的低维表征来概括。

那现在问题来了， 蛋白质水平服从什么分布呢？

非特异性抗体结合带来测量误差

充分了解数据是构建准确模型的前提，CITE-seq的蛋白质read count是怎么来的呢？

👆图中我们看到细胞表面的蛋白质和DNA-tagged antibody（DTA）结合了， 这个抗体-寡核苷酸聚合物（后面简称DTA了）含有PCR handle， 所以在洗掉未结合的抗体之后， 我们利用就可以通过PCR扩增出结合的片段，再利用Ab-Barcode鉴定出这段read对应的蛋白。

由于抗体结合存在非特异性， 所以DTA也可能结合到别的蛋白质上，当然非特异性结合的亲和力不如原靶点高。在一个原靶点蛋白不表达的细胞， 该蛋白在CITE-seq数据里却可能是有reads count的🙀。 非特异性结合的特点导致了这个reads count不会很高。 如果把所有细胞的该蛋白水平都画出来，我们就会发现histogram大致呈现几个峰。低read count的峰被称之为background（BG），一般认为是技术误差导致的； 高read count的峰称之为foreground（FG）， 更可能代表特异性结合。

了解了蛋白分布的特性， 我们就能设计出更好的模型来表述这种双峰的分布。 一种简单直白的办法就是设置一个阈值， 把低于该阈值的read count都过滤掉。另一种做法是用GMM来给每一个蛋白都构建一个一个mixture model， 给FG和BG用不同的component来描述。

但是问题就出现了，GMM应该设置多少个component呢？FG和BG两个？如果对所有细胞都是用同一个GMM， 那么我们实际上假设该dataset里的所有cell-type都服从一个样FG和BG的分布。同样的问题也发生在hard-thresholding上， 如果我们给不同的cell-type都设置同一个阈值合理吗？如果想要考虑cell type的差异，那我们应该设置多少个component或者要设置多少个阈值呢？

用细胞特异的混合模型来描述蛋白水平

TotalVI的解决方法是给每个细胞都设置单独的mixture model‼️ 听起来也太夸张了，这不就妥妥地over-parameterize了吗 ？不然，TotalVI还是假设只有FG和BG两个component，每个component都服从negative-binomial。只是分布的参数 ( ) 不再是对所有的细胞都一致，而是取决于

能高度概括该细胞的隐变量 还有实验的批次 。具体怎么做呢？

• 计算两个component的比例，
• 计算Foreground的调节项，
• 分配component，
• 利用Gamma-posion抽样来获得蛋白水平

  • 从gamma获得均值 co
  • 抽样

以上就是TotalVI如何构建蛋白质部分的分布。 和 就是下图所标的绿色和黄色的模块。而且在实际运行的时候， 我们只需要将 的激活层设置为sigmoid就能很简单地完成 的采样。这两个网络的输入是 ， 这样就实现了综合考虑RNA和蛋白水平之后对每个细胞设置FG的水平和FG-BG的比例。

RNA表达：从ZINB到NB

知道了蛋白部分如何建模，再用**scVI** 把RNA部分也建模了，完整的模型不就来了吗。

细心的你肯定发现了， 原本RNA项的zero-inflation不见了❔（发现了吧狗头）

我们在上一篇推送介绍scVI的时候提到， zero-inflation 的目的是为了建模单细胞表达数据的high drop-out rate。而实际上，Gamma-Poisson中的Poisson本身也是可以描述低概率事件的，所以领域内对于需不需要zero-inflation也有争议。

故事还得从2020年说起， 当时Svensson在NBT上直接抛出他的结论。

To the Editor
Potential users of single-cell RNA-sequencing (scRNA-seq) often encounter a choice between high-throughput droplet-based methods and high-sensitivity plate-based methods. There is a widespread belief that scRNA-seq will often fail to generate measurements for some genes from some cells owing to technical molecular inefficiencies. It is believed that this causes data to have an overabundance of zero values compared to what is expected from random sampling and that this effect is particularly pronounced in droplet-based methods. Here I present an investigation of published data for technical controls in droplet-based scRNA-seq experiments that demonstrates that the number of zero values in the data is consistent with common distributional models of molecule sampling counts. Thus, any additional zero values in biological data likely result from biological variation or may reflect variation in gene abundance among cell types or cell states.
 
In a study on simulation of scRNA-seq data, the investigators found droplet-based data to be sufficiently modeled using a negative binomial distribution, whereas plate-based data needed zero inflation to be accurately simulated.  
——MLASvensson, Valentine. Nature Biotechnology 38.2 (2020): 147-150

我这里偷偷懒， 直接放上原本第一段还有倒数第二段。 总的来说， 他用了一些只有同一种细胞来做droplet-bsed scRNA-seq的样本，然后研究了一下zero count的基因的数量。他的主要结论是在high droup-out rate其实是生物学上的差异导致的，而不是技术问题。在droplet-based的数据中用zero-infaltion并不会带来额外的提升，但是plate-based的还是需要的。

一年之后，Hoffmann实验室又在NBT上再次回应了这个issue（吵吵架就把NBT发了🥲)。

他们的结论也是简单明了，直接看标题就行了。实际上也不能说是吵架，因为他们的结果还是支持Svensson的。他们想要纠正的是不要直接把platform作为是否使用zero-inflation的依据。而是要依据建库的时候有没有用UMI来校正PCR扩增的偏差，没有的话就加上zero-inflation。

所以，到了2021年TotalVI发出来的时候， 他们也就因为CITE-seq数据是droplet-based而且RNA count也使用了UMI， 把模型由ZINB调整为了NB。OK，讲了这么多， 模型设计的部分终于讲完。下面我们就来看看和scVI相比，TotalVI 多了哪些功能。

填补未测量的蛋白

总的来看，TotalVI 和scVI功能总的来说是一脉相承的，基于隐变量的几个功能他都有， 基于后验分布的功能它也有，但是由于TotalVI 多了一个蛋白水平的建模，所以多了两个其特有的protein background decoupling和missing protein imputation。 前一个功能是我们已经介绍了原理，下面则着重介绍填补缺失蛋白。

首先需要明确问题，缺失的蛋白是怎么来的。CITE-seq测量的蛋白数量是有限， 不同的panel会测量不同的蛋白质集合。第一种情况，两个CITE-seq数据集使用了不同的蛋白panel，那么当他们整合在一起之后，我们就能用解耦了批次效应的隐变量来预测未测量的蛋白质水平（inference with ）。第二种情况，我们直接将一个CITE-seq数据和scRNA-seq数据进行整合， 第二个数据集直接没有任何蛋白数据。

这里的假设是，如果两个数据集能被很好地整合， 那么各自的细胞都能很好地在隐空间中混合。如果两个数据集的表征的分布很相似，那么从一个数据集中学习到的decoder就能给另一个数据集使用。我其实是用Domain Adaptation来理解这个问题的。

先验的小小改动

VAE模型的先验 ， 默认都是标准的高斯分布， scVI也继承了这个设定。而TotalVI却悄咪咪地改动了这里，它采用了probabilistic simplex作为隐变量的先验。

• scVI :
• TotalVI :

这个小小的改动，直接使得只有大部分细胞都被限制在几个极点界定的范围之内。类似于topics model， 一个细胞的隐变量可以理解为某几个主题的大杂烩。 比如说30%T细胞+50%细胞增殖+10%细胞毒性+10%淋巴细胞。而隐空间的每个维度也因为这个改进而获得了某种主题。

只要找出哪些细胞会让一个主题彻底激活到彻底失活（archetypal analysis极点分析）， 我们就可以知道隐变量的该维度和哪些基因相关，以此来给主题一些生物学意义。例如说， 文章的Fig 5i中就展示了隐空间第16维和Tranitional B cells的发育有关。第16维的激活值越高， B细胞就越不成熟，变现出更多的早起发育的marker Rag1的表达 。

总结

这一篇推送首先介绍了一下scVI的下游分析， 隐变量如何用来做数据整合，和计算细胞之间的关系。 然后是基于模型输出的后验估计。

然后我们将这个框架拓展到CITE-seq数据，并且详细介绍了TotalVI是如何建模由抗体共聚物捕捉的蛋白水平count，并且如果处理其非特异性结合带来的技术噪音。 之后又结合scVI介绍了其下游的分析。

最有趣的是，透过模型细节的改进（取消RNA部分的zero-inflation以及使用了simplex prior），我们仿佛看到了领域的进步和scVI整体的改进。说实话刚开始我感觉这篇文章有点水，但是仔细阅读之后反而察觉到这篇工作的扎实。 魔鬼往往藏在细节里啊～

说实话两篇推送分割的地方不是很好，其实完全把scVI的部分放到上一篇全部讲完会更好吧。而这篇推送磕磕绊绊地在各种拖拖拉拉中终于写完。希望下篇差异表达分析能写得更加简洁有趣。

参考文献

Lopez, Romain, et al. "Deep generative modeling for single-cell transcriptomics." Nature methods 15.12 (2018): 1053-1058.

Gayoso, Adam, et al. "Joint probabilistic modeling of single-cell multi-omic data with totalVI." Nature methods 18.3 (2021): 272-282.

Svensson, Valentine. "Droplet scRNA-seq is not zero-inflated." Nature Biotechnology  38.2 (2020): 147-150.

Cao, Yingying, et al. "UMI or not UMI, that is the question for scRNA-seq zero-inflation." Nature Biotechnology 39.2 (2021): 158-159.

Matuła, Kinga, Francesca Rivello, and Wilhelm TS Huck. "Single‐cell analysis using droplet microfluidics." : 1900188.