# What is the best way to extract a "full" batch effect corrected count matrix from scVI model?

**URL:** <https://discourse.scverse.org/t/what-is-the-best-way-to-extract-a-full-batch-effect-corrected-count-matrix-from-scvi-model/825>\
**Category:** scvi-tools\
**Tags:** scvi\
**Created:** [October 14, 2022, 3:31am UTC](https://discourse.scverse.org/t/what-is-the-best-way-to-extract-a-full-batch-effect-corrected-count-matrix-from-scvi-model/825 "2022-10-14T03:31:32Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![jjuhh](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/jjuhh/32/396_2.png) [@jjuhh](https://discourse.scverse.org/u/jjuhh)\
**Post date:** [October 14, 2022, 3:31am UTC](https://discourse.scverse.org/t/what-is-the-best-way-to-extract-a-full-batch-effect-corrected-count-matrix-from-scvi-model/825/1 "2022-10-14T03:31:32Z")

</div>

Hello!

I would like to get “full” batch corrected count matrix from trained scVI model for other scRNA analysis tools requiring count data. And as I understood, I could extract the batch corrected frequency of each transcript in each cell from `model.get_normalized_expression()` with `transform_batch` parameter. And to get the “corrected count”, I thought I should use the function with `library_size='latent'` as in the scVI paper(Expected count = 𝜌(frequency) X _I_(cell specific scaling factor))  
However, if the adata.X in the model had only hvgs which is used for training, the model return corrected count with only hvgs. But I could get full corrected matrix when I load model with adata which had full genes in adata.X.  
So, my questions are

1. Is using `get_normalized_expression()` function with `library_size='latent'` and `batch_size="my batches"` parameters the best way to get the batch corrected count matrix?
2. Is it okay to use the model trained only with hvgs(about 1000~10,000 genes) with full adata.X to get the full count matrix?
3. If all of my assumptions are valid, could I use the normalized data based on that batch corrected count for plotting expression based plots such as dot plot or feature plot?

Thanks a lot!!

---

<div class="post-metadata">

**Author:** ![Valentine\_Svensson](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/valentine_svensson/32/10_2.png) [@Valentine\_Svensson](https://discourse.scverse.org/u/Valentine_Svensson)\
**Post date:** [October 18, 2022, 3:02am UTC](https://discourse.scverse.org/t/what-is-the-best-way-to-extract-a-full-batch-effect-corrected-count-matrix-from-scvi-model/825/2 "2022-10-18T03:02:37Z")

</div>

Hi jjuhh,

1. 

First, `batch_size` is the parameter that determines how many cells are run through the model in each ‘chunk’ of data. Do you mean the parameter `transform_batch`?

Other than that, what you describe will get you scale parameters with a learned library size for each cell and gene in the model. If you want counts, you will need to take these scale parameter estimates, pass them to a count distribution, and sample from that.

The meaning of ‘correcting’ the expression estimates is that you wish to analyze the data as if it came from one single batch. Similarly, for counts, you would want to ‘correct’ the technical factor of count depth so that the count depth is the one same value for all cells.

What you can do is run `scales = get_normalized_expression(adata, transform_batch = "My favorite batch", library_size = 100000)`. Which batch you pick and which standardized library size you pick is arbitrary, the point is that all scales will be as if they came from the same batch so you can compare expression between all your samples. Again, this gives you scales that you need to pass to a function that samples from the given distribution. As an example for the sake of simplicity, if you trained the model with a parameter `gene_likelihood = 'Poisson'`, than to make counts you would sample from a Poisson distribution with your `scales` from above. Like this: `counts = torch.distributions.poisson.Poisson(scales).sample()`.

1. 

The model can only work on the genes you have trained it on. If you want to do this I would suggest training the model on all genes. If you are exclusively using it for this normalization task, the slightly more ‘blurry’ representation you get is probably not an issue.

1. 

I would suggest plotting the continuous values you get from `get_normalized_expression()` with a standardized library size directly instead. The issue with counts is that are discrete, which makes them hard to plot. The estimates of the expression levels you get from `get_normalized_expression()` will include all the normalization you would attempt to make the plot readable.

Hope this helps!  
/Valentine

---

<div class="post-metadata">

**Author:** ![adamgayoso](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/adamgayoso/32/100_2.png) [@adamgayoso](https://discourse.scverse.org/u/adamgayoso)\
**Post date:** [April 7, 2023, 11:30pm UTC](https://discourse.scverse.org/t/what-is-the-best-way-to-extract-a-full-batch-effect-corrected-count-matrix-from-scvi-model/825/3 "2023-04-07T23:30:14Z")

</div>

Valentine makes a lot of good points, I just want to standardize some of the terminology here.

# Single experimental batch case

**Denoised** - Denoised expression for scVI would represent a conditional expectation for cell n and gene g: \mathbb{E}\_{q(z\_n\mid x\_n)}[x\_{ng} | z\_n]. There are some caveats:

1. Denoised is not necessarily **normalized** , as the conditional expectation would use the observed or estimated library size of cell n
2. In the case of `gene_likelihood="zinb"`, my expression is inaccurate, it would technically be \mathbb{E}\_{q(z\_n\mid x\_n)}[x\_{ng} | z\_n, \text{not zero-inflated}], or the conditional expectation when also conditioning on the non-zero component of the ZINB distribution. This follows the original approach of the scVI manuscript where 0s captured by the ZI component were excess and undesirable. However, this hypothesis oversimplifies the biology and it’s currently recommended to train with `gene_likelihood="nb"` instead. @Valentine_Svensson can also likely explain very well why 🙂

**Normalized** - Denoised and scaled. We know that \mathbb{E}\_{q(z\_n\mid x\_n)}[x\_{ng} | z\_n] = \ell\_n [f(z\_n)]\_g , where \ell\_n is the scalar library size for cell n and f(z\_n) would be the denoised frequency of gene expression for gene g in cell n . The `library_size` param of `get_normalized_expression` allows manipulating \ell\_n.

# Multi-batch case

In the case of multiple batches f(z\_n), also called the decoder is actually represented as f(z\_n, s\_n), where s\_n is a one-hot experimental batch encoding. Therefore, while z\_n may be “batch-corrected”, the decoder always maps to the observed data, which clearly has batch specific effects.

To mitigate this, we “counterfactually decode”, or decode using f(z\_n, s\_n'), where s\_n' is the value of a different batch (different than the observed). This mechanism of decoding (`transform_batch`) allows us to produce **denoised** and **normalized** expression values.

If for some reason we wanted to create batch-corrected counts, we’d want to **denoise** and **not** normalize, and then sample from the likelihood distribution as Valentine described. This functionality is not currently in the codebase.

It’s also important to mention that `transform_batch` will not work perfectly when there is confounding in the experimental design. Making counterfactual decodings for cell types where s\_n' was not observed may have unintended consequences.

I can expand more on points of confusion, but I hope this helps clear up some of what this function does.

---

<div class="post-metadata">

**Author:** ![miyang](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/miyang/32/755_2.png) [@miyang](https://discourse.scverse.org/u/miyang)\
**Post date:** [July 21, 2023, 10:09pm UTC](https://discourse.scverse.org/t/what-is-the-best-way-to-extract-a-full-batch-effect-corrected-count-matrix-from-scvi-model/825/4 "2023-07-21T22:09:28Z")

</div>

> [@jjuhh](#):
>
> I

hello,  
Thanks for your reply. But what if we used a negative binomial instead of a Poisson ? How would you modify this line ?  
counts = torch.distributions.poisson.Poisson(scales).sample()

---

<div class="post-metadata">

**Author:** ![jjuhh](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/jjuhh/32/396_2.png) [@jjuhh](https://discourse.scverse.org/u/jjuhh)\
**Post date:** [August 16, 2023, 3:47am UTC](https://discourse.scverse.org/t/what-is-the-best-way-to-extract-a-full-batch-effect-corrected-count-matrix-from-scvi-model/825/5 "2023-08-16T03:47:23Z")

</div>

Thank you both for your detailed explanations! Standardizing the wording will be helpful for me. Thanks again!! 😊