Scanpy源码浅析之pp.calculate_qc_metrics

版本

导入Scanpy, 其版本为'1.9.1'，如果你看到的源码和下文有差异，其可能是由于版本差异。

import scanpy as sc

sc.__version__
#'1.9.1'

功能

函数pp.calculate_qc_metrics其源代码在scanpy/preprocessing/_qc.py
其主要功能为计算一些质控指标。详细指标见下文的小标题

代码解析

主要代码

以下为calculate_qc_metrics 主要逻辑代码，为方便理解主要逻辑，其中删除了一些即将废弃的，异常处理，日志打印，稀疏矩阵处理等代码。

其中质控指标计算由另外两个函数完成，我们在下文另外展示它们的代码。

代码说明

代码前几行是函数的参数设置：

def calculate_qc_metrics(
    adata: AnnData,
    *,
    expr_type: str = "counts",
    var_type: str = "genes",
    qc_vars: Collection[str] = (),
    percent_top: Optional[Collection[int]] = (50, 100, 200, 500),
    layer: Optional[str] = None,
    use_raw: bool = False,
    inplace: bool = False,
    log1p: bool = True,
) -> Optional[Tuple[pd.DataFrame, pd.DataFrame]]:

adata, expr_type, ..., log1p是函数参数， 冒号后面跟的是参数类型注解，表明这个参数应该传递什么类型的值给函数。

# def _choose_mtx_rep(adata, use_raw=False, layer=None):
#     is_layer = layer is not None
#     if use_raw and is_layer:
#         raise ValueError(
#             "Cannot use expression from both layer and raw. You provided:"
#             f"'use_raw={use_raw}' and 'layer={layer}'"
#         )
#     if is_layer:
#         return adata.layers[layer]
#     elif use_raw:
#         return adata.raw.X
#     else:
#         return adata.X

X = _choose_mtx_rep(adata, use_raw, layer)

该行代码用到参数adata, use_raw, layer, 根据参数设置来选择对应的数据。其中注释部分是调用函数源代码。

obs_metrics = describe_obs(
    adata,
    expr_type=expr_type,
    var_type=var_type,
    qc_vars=qc_vars,
    percent_top=percent_top,
    inplace=inplace,
    X=X,
    log1p=log1p,
)
var_metrics = describe_var(
    adata,
    expr_type=expr_type,
    var_type=var_type,
    inplace=inplace,
    X=X,
    log1p=log1p,
)

if not inplace:
    return obs_metrics, var_metrics

在上面的代码中分别调用这两个函数describe_obs, describe_var，来进行对基因，细胞进行质控指标的计算。
最后inplace=False 则直接返回两个质控指标。

describe_obs

函数源码

def describe_obs(
    adata: AnnData,
    *,
    expr_type: str = "counts",
    var_type: str = "genes",
    qc_vars: Collection[str] = (),
    percent_top: Optional[Collection[int]] = (50, 100, 200, 500),
    layer: Optional[str] = None,
    use_raw: bool = False,
    log1p: Optional[str] = True,
    inplace: bool = False,
    X=None,
) -> Optional[pd.DataFrame]:

    # Handle whether X is passed
    if X is None:
        X = _choose_mtx_rep(adata, use_raw, layer)
        if isspmatrix_coo(X):
            X = csr_matrix(X)  # COO not subscriptable
        if issparse(X):
            X.eliminate_zeros()
    obs_metrics = pd.DataFrame(index=adata.obs_names)
    if issparse(X):
        obs_metrics[f"n_{var_type}_by_{expr_type}"] = X.getnnz(axis=1)
    else:
        obs_metrics[f"n_{var_type}_by_{expr_type}"] = np.count_nonzero(X, axis=1)
    if log1p:
        obs_metrics[f"log1p_n_{var_type}_by_{expr_type}"] = np.log1p(
            obs_metrics[f"n_{var_type}_by_{expr_type}"]
        )
    obs_metrics[f"total_{expr_type}"] = np.ravel(X.sum(axis=1))
    if log1p:
        obs_metrics[f"log1p_total_{expr_type}"] = np.log1p(
            obs_metrics[f"total_{expr_type}"]
        )
    if percent_top:
        percent_top = sorted(percent_top)
        proportions = top_segment_proportions(X, percent_top)
        for i, n in enumerate(percent_top):
            obs_metrics[f"pct_{expr_type}_in_top_{n}_{var_type}"] = (
                proportions[:, i] * 100
            )
    for qc_var in qc_vars:
        obs_metrics[f"total_{expr_type}_{qc_var}"] = np.ravel(
            X[:, adata.var[qc_var].values].sum(axis=1)
        )
        if log1p:
            obs_metrics[f"log1p_total_{expr_type}_{qc_var}"] = np.log1p(
                obs_metrics[f"total_{expr_type}_{qc_var}"]
            )
        obs_metrics[f"pct_{expr_type}_{qc_var}"] = (
            obs_metrics[f"total_{expr_type}_{qc_var}"]
            / obs_metrics[f"total_{expr_type}"]
            * 100
        )
    if inplace:
        adata.obs[obs_metrics.columns] = obs_metrics
    else:
        return obs_metrics

处理X参数

    # Handle whether X is passed
    if X is None:
        X = _choose_mtx_rep(adata, use_raw, layer)
        if isspmatrix_coo(X):
            X = csr_matrix(X)  # COO not subscriptable
        if issparse(X):
            X.eliminate_zeros()

如果没传递X参数， 重新在adata里根据use_raw, layer获取数据。

生成obs指标DataFrame

obs_metrics = pd.DataFrame(index=adata.obs_names)

该行代码生成一个DataFrame, 其中行名 为细胞名（adata.obs_names）

n_genes_by_counts

n_genes_by_counts为每个细胞的基因表达量>0的基因数目

    if issparse(X):
        obs_metrics[f"n_{var_type}_by_{expr_type}"] = X.getnnz(axis=1)
    else:
        obs_metrics[f"n_{var_type}_by_{expr_type}"] = np.count_nonzero(X, axis=1)

该部分代码if else两个分支所作用目的是一样的，只是为了支持不同的数据类似，形成了两个分支，
该部分前面两行为支持稀疏矩阵处理，暂且不管，当前的源码解析主要关注numpy.ndarray类型。
我们可以从源码中发现n_genes_by_counts由f"n_{var_type}_by_{expr_type}"生成，其中

• var_type 为可传递改变的参数，默认为"genes"
• expr_type为可传递改变的参数，默认为"counts"

np.count_nonzero(X, axis=1)计算了每行细胞中表达量非0的基因的数量

log1p_n_genes_by_counts

    if log1p:
        obs_metrics[f"log1p_n_{var_type}_by_{expr_type}"] = np.log1p(
            obs_metrics[f"n_{var_type}_by_{expr_type}"]
        )

如果log1p为True, 则对**n_genes_by_counts **进行log1p转换处理，得到log1p_n_genes_by_counts
log1p表示 log(X+1), 为防止为0值出现（log(0)）

total_counts

obs_metrics[f"total_{expr_type}"] = np.ravel(X.sum(axis=1))

计算每个细胞的total counts

log1p_total_counts

    if log1p:
        obs_metrics[f"log1p_total_{expr_type}"] = np.log1p(
            obs_metrics[f"total_{expr_type}"]
        )

如果log1p为True, 则对total_counts 进行log1p转换处理，得到log1p_total_counts

pct_counts_in_top_{n}_genes

    if percent_top:
        percent_top = sorted(percent_top)
        proportions = top_segment_proportions(X, percent_top)
        for i, n in enumerate(percent_top):
            obs_metrics[f"pct_{expr_type}_in_top_{n}_{var_type}"] = (
                proportions[:, i] * 100
            )

percent_top默认值为 (50, 100, 200, 500)该参数用于设定寻找每个细胞中前n个基因的表达量和占总的基因中表达量和的比例。
函数top_segment_proportions用于计算这个比例。for循环将percent_top中每个n值，所计算的比例转换成百分比，并保存在obs_metrics 这个DataFrame中。

def top_segment_proportions(
    mtx: Union[np.array, spmatrix], ns: Collection[int]
) -> np.ndarray:

    # Pretty much just does dispatch
    if not (max(ns) <= mtx.shape[1] and min(ns) > 0):
        raise IndexError("Positions outside range of features.")
    if issparse(mtx):
        if not isspmatrix_csr(mtx):
            mtx = csr_matrix(mtx)
        return top_segment_proportions_sparse_csr(mtx.data, mtx.indptr, np.array(ns))
    else:
        return top_segment_proportions_dense(mtx, ns)

def top_segment_proportions_dense(
    mtx: Union[np.array, spmatrix], ns: Collection[int]
) -> np.ndarray:
    # Currently ns is considered to be 1 indexed
    ns = np.sort(ns)
    sums = mtx.sum(axis=1)
    partitioned = np.apply_along_axis(np.partition, 1, mtx, mtx.shape[1] - ns)[:, ::-1][
        :, : ns[-1]
    ]
    values = np.zeros((mtx.shape[0], len(ns)))
    acc = np.zeros(mtx.shape[0])
    prev = 0
    for j, n in enumerate(ns):
        acc += partitioned[:, prev:n].sum(axis=1)
        values[:, j] = acc
        prev = n
    return values / sums[:, None]

top_segment_proportions源码见上面， 其中根据传入矩阵类型分别调用了两个函数进行处理:

• top_segment_proportions_sparse_csr处理sparse 矩阵
• top_segment_proportions_dense处理dense矩阵

我们关注下dense矩阵处理方式，理解top_segment_proportions_dense源码，有几个要点：

• np.apply_along_axis函数作用
• np.partition函数作用
• [:, ::-1]取反操作
• 其他代码

qc_vars 相关指标计算

    for qc_var in qc_vars:
        obs_metrics[f"total_{expr_type}_{qc_var}"] = np.ravel(
            X[:, adata.var[qc_var].values].sum(axis=1)
        )
        if log1p:
            obs_metrics[f"log1p_total_{expr_type}_{qc_var}"] = np.log1p(
                obs_metrics[f"total_{expr_type}_{qc_var}"]
            )
        obs_metrics[f"pct_{expr_type}_{qc_var}"] = (
            obs_metrics[f"total_{expr_type}_{qc_var}"]
            / obs_metrics[f"total_{expr_type}"]
            * 100
        )

qc_vars 用于指定adata.var里的特定字段，该字段需要为布尔值，来进行相关指标计算。例如，假设adata.var有个字段为"mt" 用于判断基因是否为线粒体基因。就会得到三个相关指标：

• total_counts_mt 细胞中，线粒体基因表达量总和
• log1p_total_counts_mt log1p(细胞中线粒体基因表达量总和)
• pct_counts_mt 细胞中，线粒体基因表达量总和 占总的基因表达和的百分比

返回

    if inplace:
        adata.obs[obs_metrics.columns] = obs_metrics
    else:
        return obs_metrics

若是inplace为真，则将计算的这些指标添加到adata.obs， 否则直接返回指标数据

describe_var

函数源码

def describe_var(
    adata: AnnData,
    *,
    expr_type: str = "counts",
    var_type: str = "genes",
    layer: Optional[str] = None,
    use_raw: bool = False,
    inplace=False,
    log1p=True,
    X=None,
) -> Optional[pd.DataFrame]:

    # Handle whether X is passed
    if X is None:
        X = _choose_mtx_rep(adata, use_raw, layer)
        if isspmatrix_coo(X):
            X = csr_matrix(X)  # COO not subscriptable
        if issparse(X):
            X.eliminate_zeros()
    var_metrics = pd.DataFrame(index=adata.var_names)
    if issparse(X):
        # Current memory bottleneck for csr matrices:
        var_metrics["n_cells_by_{expr_type}"] = X.getnnz(axis=0)
        var_metrics["mean_{expr_type}"] = mean_variance_axis(X, axis=0)[0]
    else:
        var_metrics["n_cells_by_{expr_type}"] = np.count_nonzero(X, axis=0)
        var_metrics["mean_{expr_type}"] = X.mean(axis=0)
    if log1p:
        var_metrics["log1p_mean_{expr_type}"] = np.log1p(
            var_metrics["mean_{expr_type}"]
        )
    var_metrics["pct_dropout_by_{expr_type}"] = (
        1 - var_metrics["n_cells_by_{expr_type}"] / X.shape[0]
    ) * 100
    var_metrics["total_{expr_type}"] = np.ravel(X.sum(axis=0))
    if log1p:
        var_metrics["log1p_total_{expr_type}"] = np.log1p(
            var_metrics["total_{expr_type}"]
        )
    # Relabel
    new_colnames = []
    for col in var_metrics.columns:
        new_colnames.append(col.format(**locals()))
    var_metrics.columns = new_colnames
    if inplace:
        adata.var[var_metrics.columns] = var_metrics
    else:
        return var_metrics

处理X参数

    # Handle whether X is passed
    if X is None:
        X = _choose_mtx_rep(adata, use_raw, layer)
        if isspmatrix_coo(X):
            X = csr_matrix(X)  # COO not subscriptable
        if issparse(X):
            X.eliminate_zeros()

如果没传递X参数， 重新在adata里根据use_raw, layer获取数据。

生成var指标DataFrame

var_metrics = pd.DataFrame(index=adata.var_names)

该行代码生成一个DataFrame, 其中行名 为基因名（adata.var_names）

n_cells_by_counts和mean_counts

    if issparse(X):
        # Current memory bottleneck for csr matrices:
        var_metrics["n_cells_by_{expr_type}"] = X.getnnz(axis=0)
        var_metrics["mean_{expr_type}"] = mean_variance_axis(X, axis=0)[0]
    else:
        var_metrics["n_cells_by_{expr_type}"] = np.count_nonzero(X, axis=0)
        var_metrics["mean_{expr_type}"] = X.mean(axis=0)

• n_cells_by_counts 为计算 所有细胞中表达该基因的的细胞数目
• mean_counts 为所有细胞中的该基因表达量的平均值

log1p_mean_counts

    if log1p:
        var_metrics["log1p_mean_{expr_type}"] = np.log1p(
            var_metrics["mean_{expr_type}"]
        )

log1p(mean_counts)

pct_dropout_by_counts

    var_metrics["pct_dropout_by_{expr_type}"] = (
        1 - var_metrics["n_cells_by_{expr_type}"] / X.shape[0]
    ) * 100

n_cells_by_counts 为 所有细胞中表达该基因的的细胞数目， pct_dropout_by_counts为细胞中未表达基因占总的细胞总数的百分比

total_counts

var_metrics["total_{expr_type}"] = np.ravel(X.sum(axis=0))

所有细胞中，基因的表达量总和

log1p_total_counts

    if log1p:
        var_metrics["log1p_total_{expr_type}"] = np.log1p(
            var_metrics["total_{expr_type}"]
        )

log1p(所有细胞中基因的表达量总和)

收尾

    # Relabel
    new_colnames = []
    for col in var_metrics.columns:
        new_colnames.append(col.format(**locals()))
    var_metrics.columns = new_colnames

    if inplace:
        adata.var[var_metrics.columns] = var_metrics
    else:
        return var_metrics

若是inplace为真，则将计算的这些指标添加到adata.var， 否则直接返回指标数据