Data integration and batch correction with SIMBA¶
Here we will use three scRNA-seq human pancreas datasets of different studies as an example to illustrate how SIMBA performs scRNA-seq batch correction for multiple batches
We follow the corresponding tutorial at SIMBA. We do not provide much explanation, and instead refer to the original tutorial.
In [1]:
Copied!
import omicverse as ov
from omicverse.utils import mde
workdir = 'result_human_pancreas'
ov.utils.ov_plot_set()
import omicverse as ov
from omicverse.utils import mde
workdir = 'result_human_pancreas'
ov.utils.ov_plot_set()
/Users/fernandozeng/miniforge3/envs/scbasset/lib/python3.8/site-packages/phate/__init__.py
We need to install simba at first
conda install -c bioconda simba
or
pip install git+https://github.com/huidongchen/simba
pip install git+https://github.com/pinellolab/simba_pbg
Read data¶
The anndata object was concat from three anndata in simba: simba.datasets.rna_baron2016(), simba.datasets.rna_segerstolpe2016(), and simba.datasets.rna_muraro2016()
It can be downloaded from figshare: https://figshare.com/ndownloader/files/41418600
In [2]:
Copied!
adata=ov.utils.read('simba_adata_raw.h5ad')
adata=ov.utils.read('simba_adata_raw.h5ad')
We need to set workdir to initiate the pySIMBA object
In [3]:
Copied!
simba_object=ov.single.pySIMBA(adata,workdir)
simba_object=ov.single.pySIMBA(adata,workdir)
simba have been install version: 1.2
Preprocess¶
Follow the raw tutorial, we set the paragument as default.
In [4]:
Copied!
simba_object.preprocess(batch_key='batch',min_n_cells=3,
method='lib_size',n_top_genes=3000,n_bins=5)
simba_object.preprocess(batch_key='batch',min_n_cells=3,
method='lib_size',n_top_genes=3000,n_bins=5)
Before filtering: 8569 cells, 15558 genes Filter genes based on min_n_cells After filtering out low-expressed genes: 8569 cells, 14689 genes 3000 variable genes are selected. Before filtering: 2122 cells, 15558 genes Filter genes based on min_n_cells After filtering out low-expressed genes: 2122 cells, 14766 genes 3000 variable genes are selected. Before filtering: 2127 cells, 15558 genes Filter genes based on min_n_cells After filtering out low-expressed genes: 2127 cells, 15208 genes 3000 variable genes are selected.
Generate a graph for training¶
Observations and variables within each Anndata object are both represented as nodes (entities).
the data store in simba_object.uns['simba_batch_edge_dict']
In [5]:
Copied!
simba_object.gen_graph()
simba_object.gen_graph()
#shared features: 2894 Performing randomized SVD ... Searching for mutual nearest neighbors ... 26156 edges are selected #shared features: 2966 Performing randomized SVD ... Searching for mutual nearest neighbors ... 25924 edges are selected relation0: source: C, destination: G #edges: 1032345 relation1: source: C, destination: G #edges: 786551 relation2: source: C, destination: G #edges: 390188 relation3: source: C, destination: G #edges: 154188 relation4: source: C, destination: G #edges: 34417 relation5: source: C2, destination: G #edges: 687963 relation6: source: C2, destination: G #edges: 404623 relation7: source: C2, destination: G #edges: 197409 relation8: source: C2, destination: G #edges: 73699 relation9: source: C2, destination: G #edges: 15752 relation10: source: C3, destination: G #edges: 752037 relation11: source: C3, destination: G #edges: 377614 relation12: source: C3, destination: G #edges: 180169 relation13: source: C3, destination: G #edges: 77739 relation14: source: C3, destination: G #edges: 13948 relation15: source: C, destination: C2 #edges: 26156 relation16: source: C, destination: C3 #edges: 25924 Total number of edges: 5230722 Writing graph file "pbg_graph.txt" to "./result_simba/pbg/graph0" ... Finished.
PBG training¶
Before training, let’s take a look at the current parameters:
- dict_config['workers'] = 12 #The number of CPUs.
In [10]:
Copied!
simba_object.train(num_workers=6)
simba_object.train(num_workers=6)
Auto-estimating weight decay ... `.settings.pbg_params['wd']` has been updated to 0.006774 Weight decay being used for training is 0.006774 Converting input data ... [2023-06-30 22:44:21.584124] Found some files that indicate that the input data has already been preprocessed, not doing it again. [2023-06-30 22:44:21.584434] These files are in: result_human_pancreas/pbg/graph0/input/entity, result_human_pancreas/pbg/graph0/input/edge Starting training ... Finished
In [6]:
Copied!
simba_object.load('result_human_pancreas/pbg/graph0')
simba_object.load('result_human_pancreas/pbg/graph0')
Batch correction¶
Here, we use simba_object.batch_correction() to perform the batch correction
Note
If the batch is greater than 10, then the batch correction is less effective
In [7]:
Copied!
adata=simba_object.batch_correction()
adata
adata=simba_object.batch_correction()
adata
Performing softmax transformation for query data 0; Performing softmax transformation for query data 1;
Out[7]:
AnnData object with n_obs × n_vars = 12818 × 15558
obs: 'cell_type1', 'batch'
var: 'feature_symbol'
uns: 'simba_batch_dict', 'simba_batch_edge_dict', 'simba_Gen'
obsm: 'X_simba'
Visualize¶
We also use mde instead umap to visualize the result
In [8]:
Copied!
adata.obsm["X_mde"] = mde(adata.obsm["X_simba"])
adata.obsm["X_mde"] = mde(adata.obsm["X_simba"])
In [11]:
Copied!
sc.pl.embedding(adata,basis='X_mde',color=['cell_type1','batch'])
sc.pl.embedding(adata,basis='X_mde',color=['cell_type1','batch'])
Certainly, umap can also be used to visualize
In [10]:
Copied!
import scanpy as sc
sc.pp.neighbors(adata, use_rep="X_simba")
sc.tl.umap(adata)
sc.pl.umap(adata,color=['cell_type1','batch'])
import scanpy as sc
sc.pp.neighbors(adata, use_rep="X_simba")
sc.tl.umap(adata)
sc.pl.umap(adata,color=['cell_type1','batch'])
computing neighbors
finished: added to `.uns['neighbors']`
`.obsp['distances']`, distances for each pair of neighbors
`.obsp['connectivities']`, weighted adjacency matrix (0:00:01)
computing UMAP
finished: added
'X_umap', UMAP coordinates (adata.obsm) (0:00:05)
In [ ]:
Copied!