# TOTALVI RNA/protein analysis for R users

**URL:** <https://discourse.scverse.org/t/totalvi-rna-protein-analysis-for-r-users/70>\
**Category:** scvi-tools\
**Created:** [April 6, 2021, 9:38pm UTC](https://discourse.scverse.org/t/totalvi-rna-protein-analysis-for-r-users/70 "2021-04-06T21:38:22Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![tkapello](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/tkapello/32/39_2.png) [@tkapello](https://discourse.scverse.org/u/tkapello)\
**Post date:** [April 6, 2021, 9:38pm UTC](https://discourse.scverse.org/t/totalvi-rna-protein-analysis-for-r-users/70/1 "2021-04-06T21:38:22Z")

</div>

Hi all,

I have a dataset with two modalities (RNA, protein) which I want to analyze using TOTALIV. I am more familiar with R so ideally I would like to run the model in Python and then transfer the latent repsresentations and denoised counts to continue my analysis there. I went through the user guide and I have a couple of questions that I needed some clarifications. I attach them below for anyone who may be able to assist:

1. how can one store the protein raw counts in the _adata_ file? It is obvious for the RNA counts

> adata.layers[“counts”] = adata.X.copy()

but not for the protein ones.

1. Is it possible to extract the denoised counts for RNA, proteins to add in a Seurat object?

2. Is it possible to calculate DE genes/proteins for each of the clusters AND one additional condition? For example, cluster 1 stim vs ctrl? I can see that this information can be filled in the

> groupby

argument, but I understand this can only apply for a one vs all comparison.

Thanks very much in advance,

Theo

---

<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, 2021, 4:16am UTC](https://discourse.scverse.org/t/totalvi-rna-protein-analysis-for-r-users/70/2 "2021-04-07T04:16:38Z")

</div>

> [@tkapello](#):
>
> how can one store the protein raw counts in the _adata_ file? It is obvious for the RNA counts

Something like `adata.obsm["protein_counts"] = protein_data_frame ` would work. But it all depends how you’re loading the protein data in python. Is it cell ranger output? If so we have [this](https://docs.scvi-tools.org/en/stable/api/reference/scvi.data.organize_cite_seq_10x.html#scvi.data.organize_cite_seq_10x) function to help organize.

> [@tkapello](#):
>
> Is it possible to extract the denoised counts for RNA, proteins to add in a Seurat object?

Yes, if you follow the tutorial command `rna, protein = vae.get_normalized_expression`, the `rna` and `protein` objects should be dataframes that you can write to csvs with Pandas.

> [@tkapello](#):
>
> Is it possible to calculate DE genes/proteins for each of the clusters AND one additional condition? For example, cluster 1 stim vs ctrl? I can see that this information can be filled in the

Yes, you can put the category of group 1 and group 2 (in a group 1 versus group 2 test) where group 1 and group 2 are string categories present in the `groupby` cell metadata vector.

See more about these options [here](https://docs.scvi-tools.org/en/stable/api/reference/scvi.model.TOTALVI.differential_expression.html#scvi.model.TOTALVI.differential_expression). As an example, `groupby="cell_types"` and `group1="cd4"` and `group2="cd8"`

---

<div class="post-metadata">

**Author:** ![tkapello](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/tkapello/32/39_2.png) [@tkapello](https://discourse.scverse.org/u/tkapello)\
**Post date:** [April 7, 2021, 9:08am UTC](https://discourse.scverse.org/t/totalvi-rna-protein-analysis-for-r-users/70/3 "2021-04-07T09:08:30Z")

</div>

Adam, thank you for the replies.

> Something like `adata.obsm["protein_counts"] = protein_data_frame ` would work. But it all depends how you’re loading the protein data in python. Is it cell ranger output? If so we have [this](https://docs.scvi-tools.org/en/stable/api/reference/scvi.data.organize_cite_seq_10x.html#scvi.data.organize_cite_seq_10x) function to help organize.

I was planning to load the protein data from the respective Seurat object slot. So, it would be:

> In R:
> 
> adata\<-sc$AnnData(X=t(as.matrix(pbmc@assays$RNA@counts)),  
> obs=pbmc[],  
> var=GetAssay(pbmc)[])
> 
> adata$obsm[‘protein\_counts’]=data.frame(pbmc@assays$ADT@counts)

Then save as ` .h5ad` and read in Python where I would set the up the anndata and run the model as shown in the user guide.

> Yes, you can put the category of group 1 and group 2 (in a group 1 versus group 2 test) where group 1 and group 2 are string categories present in the `groupby` cell metadata vector. See more about these options [here](https://docs.scvi-tools.org/en/stable/api/reference/scvi.model.TOTALVI.differential_expression.html#scvi.model.TOTALVI.differential_expression). As an example, `groupby="cell_types"` and `group1="cd4"` and `group2="cd8"`

Here I maybe did not explain myself well. Let’s assume I have a ‘disease’ factor with two levels: ‘control’ and ‘disease’. I also have three cell types, ‘cd4’, ‘cd8’, b cell’. I want to compare control vs disease for each of the three cell types. From the `TOTALVI.` `differential_expression` documentation, I understand that the `group1` and `group2` arguments must be levels of the same parameter (defined by `groupby`). Thus, I could only think of concatenating the two factors into one string for such a comparison. Would there be an alternative to this?

---

<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, 2021, 5:43pm UTC](https://discourse.scverse.org/t/totalvi-rna-protein-analysis-for-r-users/70/4 "2021-04-07T17:43:22Z")

</div>

> [@tkapello](#):
>
> adata$obsm[‘protein\_counts’]=data.frame(pbmc@assays$ADT@counts)

I’ve had trouble with this line in reticulate, you can make protein into a separate anndata and then combine in Python like

```python
adata.obsm["protein"] = protein_adata.to_df()

```

> [@tkapello](#):
>
> Would there be an alternative to this?

What you suggested would work, but you could also use the indices argument to manually select the cells of group1 and group2.

---

<div class="post-metadata">

**Author:** ![tkapello](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/tkapello/32/39_2.png) [@tkapello](https://discourse.scverse.org/u/tkapello)\
**Post date:** [April 8, 2021, 3:12pm UTC](https://discourse.scverse.org/t/totalvi-rna-protein-analysis-for-r-users/70/5 "2021-04-08T15:12:16Z")

</div>

Ok got it and all working fine now so far! 😀

Just a couple more (hopefully last) clarifications as I get more familiar with the tool.

1. In the DGE analysis of a dual modality experiment, the `model.differential_expression` compares features from both the e.g. RNA and protein and uses the corrected/denoised data for this comparison, am I correct? The reason why I am asking is that the results do not include any protein markers, but a few integers which I am not sure where they come from.

2. Reading other members’ posts, I wanted to ask what your recommendation for the `n_number` argument should be. Default is 25, but it feels to me like a small posterior sample for a DE comparison. On the other hand, higher values would be computationally more costly. Do you have a feeling for this?

---

<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 9, 2021, 5:27am UTC](https://discourse.scverse.org/t/totalvi-rna-protein-analysis-for-r-users/70/6 "2021-04-09T05:27:23Z")

</div>

> [@tkapello](#):
>
> The reason why I am asking is that the results do not include any protein markers, but a few integers which I am not sure where they come from

I think what is happening is that when you put the protein data inside the anndata, it’s a numpy array and not a pandas dataframe. If a dataframe is not provided scvi-tools will give numerical names to each protein. So “14” is the 15th column of the protein matrix (0 indexed in python)

> [@tkapello](#):
>
> Default is 25, but it feels to me like a small posterior sample for a DE comparison. On the other hand, higher values would be computationally more costly. Do you have a feeling for this?

You need to quadruple the number of samples to halve the standard deviation of the estimate. Too many samples becomes too computationally expensive and from our experience the variance is already fairly low. You should feel free to experiment with it though. At the end of the day the entire method is stochastic, so variance from run to run might even be higher than the monte carlo estimator variance.