微精选

收到一个粉丝的提问，如下：

这里直接看 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.ndarray with   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)

今天分享到这~