标签：pp log1p expr obs metrics 源码 var type 浅析

版本

导入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， 否则直接返回指标数据

标签：pp,log1p,expr,obs,metrics,源码,var,type,浅析
来源： https://www.cnblogs.com/huanping/p/16684870.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。