Scanpy源码浅析之pp.filter_cells

版本

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

import scanpy as sc

sc.__version__
#'1.9.1'

功能

函数pp.filter_cells用于Filter cell outliers， 其源代码在scanpy/preprocessing/_simple.py
其主要功能为过滤异常的细胞。可以根据设置细胞的最大/最小表达量设置阈值，将大于/小于这个阈值的细胞过滤。也可根据细胞表达基因的最大/最小数量设定阈值，将大于/小于这个阈值的细胞给过滤掉。
参数设置：

• data: AnnData 对象
• min_counts细胞的最小表达counts
• min_genes细胞的基因最小表达数目
• max_counts细胞的最大表达counts
• max_genes细胞的基因最大表达数目
• inplace是否在原数据对象上修改

代码解析

主要代码

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

参数设置

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

def filter_cells(
    data: AnnData,
    min_counts: Optional[int] = None,
    min_genes: Optional[int] = None,
    max_counts: Optional[int] = None,
    max_genes: Optional[int] = None,
    inplace: bool = True,
    copy: bool = False,
) -> Optional[Tuple[np.ndarray, np.ndarray]]:

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

参数使用限制

n_given_options = sum(
    option is not None for option in [min_genes, min_counts, max_genes, max_counts]
)
if n_given_options != 1:
    raise ValueError(
        'Only provide one of the optional parameters `min_counts`, '
        '`min_genes`, `max_counts`, `max_genes` per call.'
    )

上面几行代码， 限定了min_genes, min_counts, max_genes, max_counts，四个参数不能同时使用，调用函数时，只能使用一个。

递归获取结果及处理

if isinstance(data, AnnData):
    adata = data.copy() if copy else data
    # materialize_as_ndarray是用于数据类型转换，并非过滤操作的逻辑代码，可忽略
    cell_subset, number = materialize_as_ndarray(
        filter_cells(adata.X, min_counts, min_genes, max_counts, max_genes)
    )
    if not inplace:
        return cell_subset, number
    if min_genes is None and max_genes is None:
        adata.obs['n_counts'] = number
    else:
        adata.obs['n_genes'] = number
    adata._inplace_subset_obs(cell_subset)
    return adata if copy else None

这部分代码 对输入data 类型做了一个判断， 如果为AnnData， 则将data.X 重新传入filter_cells函数本身来递归调用，得到过滤后的结果。
如果inplace=False, 则直接返回cell_subset, number
否则， 继续在adata操作，根据过滤条件参数， 添加n_counts或者n_genes
然后adata._inplace_subset_obs(cell_subset)是用过滤后的数据取代原来的数据

过滤操作

        
X = data  # proceed with processing the data matrix
min_number = min_counts if min_genes is None else min_genes
max_number = max_counts if max_genes is None else max_genes
number_per_cell = np.sum(
    X if min_genes is None and max_genes is None else X > 0, axis=1
)
# 稀疏矩阵的处理代码
if issparse(X):
    number_per_cell = number_per_cell.A1
if min_number is not None:
    cell_subset = number_per_cell >= min_number
if max_number is not None:
    cell_subset = number_per_cell <= max_number
    
# 返回细胞布尔索引，每个细胞的counts 总和/基因表达数量
return cell_subset, number_per_cell

这部分为实际进行过滤筛选的代码， 其中data不再为初始的AnnData对象的data,而是二次调用传入的adata.X。
下面两行对变量做个统一
min_genes，min_counts 统一为min_number； max_genes，max_genes统一为max_number

min_number = min_counts if min_genes is None else min_genes
max_number = max_counts if max_genes is None else max_genes

如果min_genes或 max_genes都没设置， 则计算每行的total count， 否则则判断基因是否>0，得到真值, 计算每行基因的表达数量。

number_per_cell = np.sum(
    X if min_genes is None and max_genes is None else X > 0, axis=1
)

根据min_number或max_number，进行大小判断，得到cell 行的布尔索引

if min_number is not None:
    cell_subset = number_per_cell >= min_number
if max_number is not None:
    cell_subset = number_per_cell <= max_number