版本
导入Scanpy, 其版本为'1.9.1'，如果你看到的源码和下文有差异，其可能是由于版本差异。
import scanpy as sc
sc.__version__
#'1.9.1'
功能
函数pp.calculate_qc_metrics其源代码在scanpy/preprocessing/_qc.py
其主要功能为计算一些质控指标。详细指标见下文的小标题
代码解析
主要代码
以下为calculate_qc_metrics 主要逻辑代码，为方便理解主要逻辑，其中删除了一些即将废弃的，异常处理，日志打印，稀疏矩阵处理等代码。
[图片上传失败...(image-defb81-1663405843670)]
其中质控指标计算由另外两个函数完成，我们在下文另外展示它们的代码。
代码说明
代码前几行是函数的参数设置：
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， 否则直接返回指标数据
