Scanpy源码浅析之pp.filter_genes

版本

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

import scanpy as sc

sc.__version__
#'1.9.1'

功能

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

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

代码解析

主要代码

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

参数设置

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

def filter_genes(
    data: AnnData,
    min_counts: Optional[int] = None,
    min_cells: Optional[int] = None,
    max_counts: Optional[int] = None,
    max_cells: Optional[int] = None,
    inplace: bool = True,
    copy: bool = False,
) -> Union[AnnData, None, Tuple[np.ndarray, np.ndarray]]:

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

参数使用限制

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

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

递归获取结果及处理

if isinstance(data, AnnData):
    adata = data.copy() if copy else data
    gene_subset, number = materialize_as_ndarray(
        filter_genes(
            adata.X,
            min_cells=min_cells,
            min_counts=min_counts,
            max_cells=max_cells,
            max_counts=max_counts,
        )
    )
    if not inplace:
        return gene_subset, number
    if min_cells is None and max_cells is None:
        adata.var['n_counts'] = number
    else:
        adata.var['n_cells'] = number
    adata._inplace_subset_var(gene_subset)
    return adata if copy else None

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

过滤操作

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

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

min_number = min_counts if min_cells is None else min_cells
max_number = max_counts if max_cells is None else max_cells

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

number_per_gene = np.sum(
    X if min_cells is None and max_cells is None else X > 0, axis=0
)

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

if min_number is not None:
    gene_subset = number_per_gene >= min_number
if max_number is not None:
    gene_subset = number_per_gene <= max_number