# Insufficient batch correction for certain cell-types

**URL:** <https://discourse.scverse.org/t/insufficient-batch-correction-for-certain-cell-types/2237>\
**Category:** scvi-tools\
**Tags:** integration, scvi\
**Created:** [April 24, 2024, 9:08am UTC](https://discourse.scverse.org/t/insufficient-batch-correction-for-certain-cell-types/2237 "2024-04-24T09:08:59Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![dkioroglou](https://avatars.discourse-cdn.com/v4/letter/d/dfb087/32.png) [@dkioroglou](https://discourse.scverse.org/u/dkioroglou)\
**Post date:** [April 24, 2024, 9:08am UTC](https://discourse.scverse.org/t/insufficient-batch-correction-for-certain-cell-types/2237/1 "2024-04-24T09:08:59Z")

</div>

Hello,

## Dataset

I have a scRNASeq dataset that was based on the following experimental design:

| timepoint | Nsamples |
| --- | --- |
| t0 | 5 |
| t1 | 3 |
| t2 | 4 |

## Aim

Correct for interindividual variability across timepoints and study gene expression for each cell-type between timepoints.

## Batch correction

In order to remove interindividual variability, first I created two columns:

1. One column named **Replicate** containing sample IDs.
2. Second column denoting the timepoint **t1** , **t2** and **t3**.

then I performed batch correction as following:

1. Select 4000 highly variable genes

```auto
sc.pp.highly_variable_genes(adata, n_top_genes=4000, subset = True, layer = 'soupX_counts', flavor = "seurat_v3", batch_key="Replicate")

```

`soupX_counts` were provided by the `SoupX` package duing the QC steps to remove any possible mRNA contamination.

1. Setup SCVI model

```auto
scvi.model.SCVI.setup_anndata(adata, layer = "soupX_counts",
                             categorical_covariate_keys=["Replicate"],
                             continuous_covariate_keys=['pct_counts_mt', 'total_counts'])

```

I tried tweaking the parameters (as suggested in other posts on batch correction) and the following provided the best results for me:

```auto
model = scvi.model.SCVI(
        adata,
        n_layers=2,
        n_hidden=200,
        n_latent=20,
        gene_likelihood='zinb',
        dispersion="gene"
        )

```

1. Define training size

```auto
if 0.1 * adata.n_obs < 20000:
    train_size = 0.9
else:
    train_size = 1-(20000/adata.n_obs)

```

1. Train SCVI model

```auto
model.train(
        accelerator="gpu", 
        devices=1,
        early_stopping=True,
        train_size=train_size,
        early_stopping_patience=400,
        max_epochs=10000,
        batch_size=1024,
        limit_train_batches=20
        )

```

## Results

Despite the batch correction step, still there seems to be some batch effect for certain cell-types especially **classical monocytes**.

 ![all_PRE_UMAP_high_resolution](https://canada1.discourse-cdn.com/flex035/uploads/forum11/original/2X/0/05f909999d178d1252549ce8f284265274e724ea.jpeg)

In the figure above cells from timepoint **t0** (PRE) were selected and UMAP was performed with `scanpy`. Cell-type annotation was performed with `CellTypist`.

## Questions

1. Do you think the observed batch effect appears strong enough to affect differential expression?
2. Did the batch correction get affected negatively by the unbalanced dataset?
3. Do you have any suggestion to improve the batch correction process?

Thank you in advance.

---

<div class="post-metadata">

**Author:** ![martinkim0](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/martinkim0/32/881_2.png) [@martinkim0](https://discourse.scverse.org/u/martinkim0)\
**Post date:** [April 24, 2024, 3:42pm UTC](https://discourse.scverse.org/t/insufficient-batch-correction-for-certain-cell-types/2237/2 "2024-04-24T15:42:39Z")

</div>

Hi, I have a couple of clarifying questions:

- What is the size of your dataset? Depending on that, the number of maximum epochs might be too large. How long does it take for early stopping to kick in?
- What’s the reason behind setting `limit_train_batches=20`?

---

<div class="post-metadata">

**Author:** ![dkioroglou](https://avatars.discourse-cdn.com/v4/letter/d/dfb087/32.png) [@dkioroglou](https://discourse.scverse.org/u/dkioroglou)\
**Post date:** [April 24, 2024, 7:23pm UTC](https://discourse.scverse.org/t/insufficient-batch-correction-for-certain-cell-types/2237/3 "2024-04-24T19:23:33Z")

</div>

The size of the merged dataset (i.e t0+t1+t2) is 86500 cells. The activation of the early stopping depended on the `early_stopping_patience` value. When the latter was set at 45, early stopping happened within 7 minutes, at 150 within 20 minutes and at 400 close to 50 minutes.

Regarding the `limit_train_batches`, I’m not sure whether my reasoning is correct but the idea was to constrain RAM usage and instead of loading all cells during training, I wanted to load batches of 20 till all dataset was used. However, I’m not sure whether this is correct. Does `limit_train_batches` gets multiplied by `batch_size` and provides the total number of cells used in each epoch? If that’s true, then I’m using only 20k cells in each epoch.

---

<div class="post-metadata">

**Author:** ![martinkim0](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/martinkim0/32/881_2.png) [@martinkim0](https://discourse.scverse.org/u/martinkim0)\
**Post date:** [April 25, 2024, 4:02pm UTC](https://discourse.scverse.org/t/insufficient-batch-correction-for-certain-cell-types/2237/4 "2024-04-25T16:02:42Z")

</div>

Could you also specify the number of epochs for which your model is training? I have a feeling that with those values of `max_epochs` and `early_stopping_patience`, your model is most likely overfitting on the training set. This can be checked by plotting `"elbo_validation"` and `"elbo_train"` in `model.history` post-training. If the validation ELBO is significantly higher than its training counterpart, then it’s likely that the model overfit. Something in the ballpark of `max_epochs=400` and `early_stopping_patience=10` might be more appropriate for this dataset size.

My understanding of `limit_train_batches` is that it will run that number of minibatches during each training epoch, so yes the model would be seeing about 20k observations during training. Your main memory usage will mostly be set by your full dataset size, and your GPU RAM usage will be set by `batch_size` and the model size, so `limit_train_batches` will only affect how long it takes for each epoch to complete.

---

<div class="post-metadata">

**Author:** ![cane11](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/cane11/32/241_2.png) [@cane11](https://discourse.scverse.org/u/cane11)\
**Post date:** [April 25, 2024, 8:44pm UTC](https://discourse.scverse.org/t/insufficient-batch-correction-for-certain-cell-types/2237/5 "2024-04-25T20:44:42Z")

</div>

Hi, I think overall it looks fine integration wise. Yes, it will likely affect DE genes between replicates (between cell-types I guess it might be fine. If you want to get classical monocytes better integrated, I would check DE gene between those batches. It’s a very typical behaviour that monocytes get activated during the experiment (in humans higher in IL1B but also FOS and JUN). scVI usually does not integrate this (it’s helpful as it would otherwise do overintegration). If you think it’s critical to get better integration, I would reduce the number of genes and n\_latent (to maybe 1500 genes). I’m also not sure whether Replicate is a unique sample or whether it has the same value for different samples at different time points. I would in that case recommend to replace Replicate by a string concatenation of Replicate and Timepoint. I hope this is clear.  
ScANVI might also be helpful to get better integration (there you inform the model about the meaningful cellular variation).  
I would generally recommend using batch\_key instead of categorical covariate key if you use a single category as it offers more downstream capabilities (see transform\_batch in the codebase).

---

<div class="post-metadata">

**Author:** ![dkioroglou](https://avatars.discourse-cdn.com/v4/letter/d/dfb087/32.png) [@dkioroglou](https://discourse.scverse.org/u/dkioroglou)\
**Post date:** [April 26, 2024, 1:34pm UTC](https://discourse.scverse.org/t/insufficient-batch-correction-for-certain-cell-types/2237/6 "2024-04-26T13:34:52Z")

</div>

Thank you both for your suggestions. I have done a couple of optimizations by combining your recommendations.

## Previous implementation

The intuition of Martin was right and my previous implementation was overfitting.

 ![model_elbo_train_and_validation_score](https://canada1.discourse-cdn.com/flex035/uploads/forum11/original/2X/1/12819adad6b547357d9886b136000d147fd9bad3.jpeg)

## Optimization 1

I introduced the following changes in the code:

1. Reduced highly variable genes to 3000:

> **Code**
>
> ```auto
> sc.pp.highly_variable_genes(adata, n_top_genes=3000, subset = True, layer = 'soupX_counts', flavor = "seurat_v3", batch_key="Replicate")
> 
> ```

1. Used `SCANVI` and replaced `categorical_covariate_keys` with `batch_key`.

> **Code**
>
> ```auto
> scvi.model.SCANVI.setup_anndata(adata, layer = "soupX_counts", 
> batch_key="Replicate",
> labels_key="celltypist_cell_types_high_resolution", 
> unlabeled_category="Unknown",
> continuous_covariate_keys=['pct_counts_mt',
> 'total_counts'])
> 
> ```

1. Reduced `max_epochs` to `400`, `early_stopping_patience` to `10` and removed `limit_train_batches` parameter.

> **Code**
>
> ```auto
> model.train(
> accelerator="gpu", 
> devices=1,
> early_stopping=True,
> train_size=train_size,
> max_epochs=400,
> early_stopping_patience=10,
> batch_size=1024,
> )
> 
> ```

**Results** : less overfitting but batch effect remained.

 ![model_elbo_train_and_validation_score](https://canada1.discourse-cdn.com/flex035/uploads/forum11/original/2X/0/015b51568373675c27b94402a3412840c4846b80.jpeg)

 ![all_PRE_UMAP_high_resolution](https://canada1.discourse-cdn.com/flex035/uploads/forum11/original/2X/5/5469e5089eb73392e7b1ec49bf0d2b44aa84020f.jpeg)

## Optimization 2

I implemented the recommendation of Can and engineered a new Replicate feature for the batch correction called `Replicate2`.

> **Code**
>
> ```auto
> adata.obs['Replicate2'] = [f"{x}_{y}" for x,y in zip(adata.obs['Replicate'], adata.obs['Condition'])]
> 
> ```

**Results** : less overfitting (but seems to need a bit more time to converge), batch effect seems to have disappeared.

 ![model_elbo_train_and_validation_score](https://canada1.discourse-cdn.com/flex035/uploads/forum11/original/2X/0/07f40e9899494533b5142605c5275c38ab61b021.jpeg)

 ![all_PRE_UMAP_high_resolution](https://canada1.discourse-cdn.com/flex035/uploads/forum11/original/2X/6/67f29446139c7c3266c031f1d2b8dcaeae96ecb0.jpeg)

## Optimization 3

I reduced `n_latent` to `10` and increased the `early_stopping_patience` to `20`.

**Results** : better convergence and best integration so far.

 ![model_elbo_train_and_validation_score](https://canada1.discourse-cdn.com/flex035/uploads/forum11/original/2X/3/37c7a2ed8c4d549223fed82c91807bc7edd9b6b9.jpeg)

 ![all_PRE_UMAP_high_resolution](https://canada1.discourse-cdn.com/flex035/uploads/forum11/original/2X/a/a2200379078f52f862f36519bc654963ecf6a7bf.jpeg)

## New problem

Although optimization 3 yielded the best results, after performing differential gene expression I got no DEGs (which is very weird). I think the issue is with the way the feature `Replicate2` has been engineered which I find quite unintuitive. Let’s me give some clarifications on the experimental design:

I have the following columns in the dataset:  
**Replicate** : `{p1,p2,p4,p5,p6}`  
**Condition** : `{t0, t1, t2}`  
The timepoints are called Condition because they concern before treatment (t0), after treatment (t1) and after many treatments (t2). So DEGs should be identified between `t0 vs t1`, `t1 vs t2`.

Based on the suggestion of Can, the new “Replicate” feature should be like following:  
**Replicate2** : `{p1_t0, p1_t1, p1_t2, ..., p6_t0, p6_t1, p6_t2}`

I think by performing batch correction using `Replicate2` all the biological signal between the timepoints got disappeared. Therefore, I think that batch correction should be performed using the `Replicate` instead of the `Replicate2` feature. Could you please clarify why correcting using `Replicate2` is better?

---

<div class="post-metadata">

**Author:** ![cane11](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/cane11/32/241_2.png) [@cane11](https://discourse.scverse.org/u/cane11)\
**Post date:** [May 7, 2024, 10:19am UTC](https://discourse.scverse.org/t/insufficient-batch-correction-for-certain-cell-types/2237/7 "2024-05-07T10:19:40Z")

</div>

Can you also share your line of code for DE (I need to know how batch\_key is used there)?  
In general batch effects exist between sequencing samples/sample processings. I assume the time points of one replicate are independent (different patients or mice). If it’s the same patient sampled multiple times using replicate for integration makes sense (otherwise your design could contain e.g. male and female for the same replicate name and this will lead to low integration).  
I’m case, it’s a single replicate across time points. I still think your first integration looks fine and your main point was how to get more integration.

---

<div class="post-metadata">

**Author:** ![dkioroglou](https://avatars.discourse-cdn.com/v4/letter/d/dfb087/32.png) [@dkioroglou](https://discourse.scverse.org/u/dkioroglou)\
**Post date:** [May 7, 2024, 11:43am UTC](https://discourse.scverse.org/t/insufficient-batch-correction-for-certain-cell-types/2237/8 "2024-05-07T11:43:08Z")

</div>

Hi Can,

Regarding the experimental design, there are 5 patients (`p1,p2,p4,p5,p6`) each one provided samples at 3 timepoints (`t0, t1, t2`) as shown in the table below.

| Replicate | Condition | Replicate2 |
| --- | --- | --- |
| p1 | t0 | p1\_t0 |
| p1 | t1 | p1\_t1 |
| p1 | t2 | p1\_t2 |
| p2 | t0 | p2\_t0 |
| p2 | t1 | p2\_t1 |
| p2 | t2 | p2\_t2 |
| p4 | t0 | p4\_t0 |
| p4 | t1 | p4\_t1 |
| p4 | t2 | p4\_t2 |
| p5 | t0 | p5\_t0 |
| p5 | t1 | p5\_t1 |
| p5 | t2 | p5\_t2 |
| p6 | t0 | p6\_t0 |
| p6 | t1 | p6\_t1 |
| p6 | t2 | p6\_t2 |

The idea is to remove inter-individual differences between patients across timepoints (Condition) and get DEGs between timepoints (Condition). My initial implementation included the column “ **Replicate** ” as `batch_key` and my last implementation included “ **Replicate2** ” as `batch_key`.

Below is my code for DEGs:

```auto
pairs = [("t0", "t1"), ("t1", "t2"), ("t0", "t2")]

def getDEGs(cellType):
    global DEGs
    for pair in pairs:
        p1, p2 = pair
        cellIDX1 = ((adata.obs['celltypist_cell_types_low_resolution'] == cellType) & (adata.obs['Condition'] == p1))
        cellIDX2 = ((adata.obs['celltypist_cell_types_low_resolution'] == cellType) & (adata.obs['Condition'] == p2))
        tmp = model.differential_expression(idx1=cellIDX1, idx2=cellIDX2, mode='change', batch_correction=True)
        DEGs[f"{cellType}_{p1}_vs_{p2}"] = tmp

DEGs = {}
for cellType in celltypesLowResolution:
    getDEGs(cellType)

```

When “ **Replicate2** ” was used as `batch_key` no DEGs were identified. To my understanding if `batch_correction` is `True` then by default `bathcid1` and `batchid2` use all the categories specified in the `batch_key` with `setup_anndata`. So, I think batch correcting with “ **Replicate2** ” leads to betten integration but removes all biological signal.

---

<div class="post-metadata">

**Author:** ![cane11](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/cane11/32/241_2.png) [@cane11](https://discourse.scverse.org/u/cane11)\
**Post date:** [May 15, 2024, 8:25pm UTC](https://discourse.scverse.org/t/insufficient-batch-correction-for-certain-cell-types/2237/9 "2024-05-15T20:25:19Z")

</div>

I think in your setup it is fine to use replicate as batch\_key. The different donors will be the reason for most of the batch effect likely, so there’s little batch effect between time points. It is fine to use Replicate2 in this case though you want to set _batch\_correction=False_ in model.differential\_expression or better set batch\_id1 to all batches of condition1 and for batch\_id2 to condition2. This way the model won’t correct gene expression for the covariate of interest and will generate those counts for the correct batch.