收到一个粉丝的提问,如下:
这里直接看 https://github.com/brianhie/scanorama,下载源码看着两个函数都干了什么。顺便通过一个示例教程看具体的结果区别。「官网有一个解释,但是我感觉这个描述跟我运行的结果不相符。两个都没有返回scanorama-transformed 的矩阵。」
直接比较结果
以这个教程为准:https://colab.research.google.com/drive/12hNry9nlMgZRu-bGUiXbz0Veqeh2WhcU?usp=sharing
数据 Small dataset demo: Jurkat and 293T cells integration with Scanorama
1.输入数据准备
示例数据下载地址:https:///record/7968485/files/small_dataset.zip
这个网址也不好下载,有需要的call me:Biotree123 或者邮箱 [email protected]
# bash
# 解压:
unzip small_dataset.zip -d scanorama_demo
# 运行环境,我这里使用前面准备好的 conda小环境 sc:
conda activate sc
下面开始都是python代码啦,如果你还不知道在哪里运行python既简单又方便,可以喊我。如果遇到报错,缺啥就安装啥。
import scanorama
import scanpy as sc
import anndata as ad
import scanpy.external as sce
import matplotlib.pyplot as plt
import numpy as np
## 读取数据,使用 scanpy包的 read_h5ad 函数:
%%time
adata = sc.read_h5ad('scanorama_demo/small_dataset/small_data_293T_Jurkat.h5ad')
adata
2.按照 batch 分割对象并生成一个list的AnnData
原始的 adata 对象被拆分成一个子集列表(adatas),其中每个子集包含 batch 列中唯一值对应的数据。
开头的 %%time 命令用于测量此操作的执行时间。
%%time
batch_key = 'batch'
adatas = [adata[adata.obs[batch_key] == batch_value].copy() for batch_value in adata.obs[batch_key].unique()]
adatas
3.batch校正
“adatas”列表中的数据集使用“scanorama.integrate_scanpy”函数以及默认参数进行整合。
在完成整合后,调用“adatas”变量以在输出中显示其内容。
整合后的结果以“X_scanorama”的形式重新存储在“adatas”变量中。
dimred:dimred: int, optional (default: 100), Dimensionality of integrated embedding.
%%time
scanorama.integrate_scanpy(adatas, dimred = 100)
adatas
这一步运行完后adatas list里面每个对象都会生成一个 “X_scanorama”的结果:a list of
numpy.ndarraywith integrated low dimensional embeddings
直接取出来看看:adatas[0].obsm['X_scanorama'],类似于pca的降维坐标,维度为100。
adatas[0]
adatas[0].obsm['X_scanorama']
adatas[0].obs_names
# 将 numpy.ndarray 转换为 pandas DataFrame
import pandas as pd
X_scanorama_df = pd.DataFrame(adatas[0].obsm['X_scanorama'], index=adatas[0].obs_names)
X_scanorama_df
4.添加 scanorama corrected matrix
将list中每个adata的低维坐标 adatas[0].obsm['X_scanorama'] 整合成一个整的:
%%time
adata_sc = adata.copy()
scanorama_int = [ad.obsm['X_scanorama'] for ad in adatas]
scanorama_int
all_s = np.concatenate(scanorama_int)
adata_sc.obsm['X_scanorama'] = all_s
adata_sc
import pandas as pd
X_scanorama_df = pd.DataFrame(adata_sc.obsm['X_scanorama'], index=adata_sc.obs_names)
X_scanorama_df
5.使用 'X_scanorama' 降维
第一行使用存储在 X_scanorama 中的整合数据计算了 adata_sc 对象中细胞的邻域图。
第二行使用 X_scanorama 的整合数据计算了 adata_sc 对象中细胞的 t 分布随机邻域嵌入(t-SNE)。
第三行计算了 adata_sc 对象中细胞的均匀流形近似和投影(UMAP)表示,使用 X_scanorama 表示中的整合数据。
%%time
sc.pp.neighbors(adata_sc, use_rep='X_scanorama')
sc.tl.tsne(adata_sc, use_rep='X_scanorama')
sc.tl.umap(adata_sc)
adata_sc
6.与非矫正数据比对
首先复制原始的 adata 对象,并将其赋值给新变量 adata_raw。
在存储于 adata_raw 对象中的数据上计算主成分分析(PCA)。使用 PCA 降维后的数据,为 adata_raw 对象中的细胞计算邻域图、t 分布随机邻域嵌入(t-SNE)和均匀流形近似和投影(UMAP)表示。
%%time
adata_raw = adata.copy()
sc.tl.pca(adata_raw)
sc.pp.neighbors(adata_raw)
sc.tl.tsne(adata_raw)
sc.tl.umap(adata_raw)
adata_raw
7.整合与非整合结果可视化
fig, axs = plt.subplots(2, 2, figsize=(10,8),constrained_layout=True)
sc.pl.embedding(adata_raw, basis='tsne', color='batch', title='Uncorrected batch t-SNE', ax=axs[0,0], show=False)
sc.pl.embedding(adata_raw, basis='tsne', color='cell_types', title='Uncorrected cell types t-sNE', ax=axs[1,0], show=False)
sc.pl.embedding(adata_sc, basis='tsne', color='batch', title='Scanorama corrected batch t-SNE', ax=axs[0,1], show=False)
sc.pl.embedding(adata_sc, basis='tsne', color='cell_types', title='Scanorama corrected cell types t-SNE', ax=axs[1,1], show=False)
如果是 correct
这个参数控制返回结果的不同
return_dimred=False
%%time
# Batch correction.
#corrected = scanorama.correct_scanpy(adatas, return_dimred=True)
corrected = scanorama.correct_scanpy(adatas, return_dimred=False)
integrate_scanpy vs correct_scanpy :
integrate_scanpy 没有返回值,直接在 adata 对象里面添加了 adata.obsm['X_scanorama']
# Integration with scanpy's AnnData object.
def integrate_scanpy(adatas, **kwargs):
"""Integrate a list of `scanpy.api.AnnData`.
Parameters
----------
adatas : `list` of `scanpy.api.AnnData`
Data sets to integrate.
kwargs : `dict`
See documentation for the `integrate()` method for a full list of
parameters to use for batch correction.
Returns
-------
None
"""
datasets_dimred, genes = integrate(
[adata.X for adata in adatas],
[adata.var_names.values for adata in adatas],
**kwargs
)
for adata, X_dimred in zip(adatas, datasets_dimred):
adata.obsm['X_scanorama'] = X_dimred
correct_scanpy 有返回值,返回的是 新整理的 new_adatas 对象
# Batch correction with scanpy's AnnData object.
def correct_scanpy(adatas, **kwargs):
"""Batch correct a list of `scanpy.api.AnnData`.
Parameters
----------
adatas : `list` of `scanpy.api.AnnData`
Data sets to integrate and/or correct.
`adata.var_names` must be set to the list of genes.
return_dimred : `bool`, optional (default=`False`)
When `True`, the returned `adatas` are each modified to
also have the integrated low-dimensional embeddings in
`adata.obsm['X_scanorama']`.
kwargs : `dict`
See documentation for the `correct()` method for a full list of
parameters to use for batch correction.
Returns
-------
corrected
By default (`return_dimred=False`), returns a list of new
`scanpy.api.AnnData`.
When `return_dimred=True`, `corrected` also includes the
integrated low-dimensional embeddings in
`adata.obsm['X_scanorama']`.
"""
if'return_dimred'in kwargs and kwargs['return_dimred']:
datasets_dimred, datasets, genes = correct(
[adata.X for adata in adatas],
[adata.var_names.values for adata in adatas],
**kwargs
)
else:
datasets, genes = correct(
[adata.X for adata in adatas],
[adata.var_names.values for adata in adatas],
**kwargs
)
from anndata import AnnData
new_adatas = []
for i in range(len((adatas))):
adata = AnnData(datasets[i])
adata.obs = adatas[i].obs
adata.obsm = adatas[i].obsm
# Ensure that variables are in the right order,
# as Scanorama rearranges genes to be in alphabetical
# order and as the intersection (or union) of the
# original gene sets.
adata.var_names = genes
gene2idx = { gene: idx for idx, gene in
zip(adatas[i].var.index,
adatas[i].var_names.values) }
var_idx = [ gene2idx[gene] for gene in genes ]
adata.var = adatas[i].var.loc[var_idx]
adata.uns = adatas[i].uns
new_adatas.append(adata)
if'return_dimred'in kwargs and kwargs['return_dimred']:
for adata, X_dimred in zip(new_adatas, datasets_dimred):
adata.obsm['X_scanorama'] = X_dimred
return new_adatas
总结
但是从上面的运行结果来看 「scanorama.integrate_scanpy(adatas, dimred = 100) 与 corrected = scanorama.correct_scanpy(adatas, return_dimred=True)的结果一样,前面的降维坐标直接添加在 adatas中,后者是新生成的一个对象 corrected。」
Integration.
scanorama.integrate_scanpy(adatas) # 返回类似pca的降维坐标,存储在 adata_sc.obsm['X_scanorama'] 里面,可以用于下游的KNN图构建,降维聚类分群
Batch correction.
corrected = scanorama.correct_scanpy(adatas) # 返回的一个新对象corrected,但是没有obsm['X_scanorama']
Integration and batch correction.
corrected = scanorama.correct_scanpy(adatas, return_dimred=True) # 返回的一个新对象corrected,但是有obsm['X_scanorama']
需要注意的是,这两个函数与 scanorama 自带的 scanorama.integrate 和 scanorama.correct 返回结果不一样:
# List of datasets (matrices of cells-by-genes):
datasets = [ list of scipy.sparse.csr_matrix or numpy.ndarray ]
# List of gene lists:
genes_list = [ list of list of string ]
import scanorama
# Integration.
integrated, genes = scanorama.integrate(datasets, genes_list)
# Batch correction.
corrected, genes = scanorama.correct(datasets, genes_list)
# Integration and batch correction.
integrated, corrected, genes = scanorama.correct(datasets, genes_list, return_dimred=True)
今天分享到这~