作者：童蒙
编辑：angelica
scanpy代码解读来啦~
单细胞分析第一步是对数据进行标准化，标准化的方法有很多，下面给大家解读一下scanpy的一个：函数为：scanpy.pp.normalize_total
作用
对每个细胞的count进行标准化，使得每个细胞标准化后有相同的count
- 如果使用target_sum=1e6，那么相当于CPM标准化。
- 设置exclude_highly_expressed=True 后，非常高的表达基因会被排除到计算size factor中，这个会影响其他的基因。需要同max_fraction进行连用，只要有一个细胞里面超过了，就会被判别为高表达基因。默认为0.05。
常见参数
- adata：内置的AnnDta数据
- target_sum: 如果设置为none，那么会用所有样品的median值来替代；
- exclude_highly_expressed ：是否去除高表达基因
- max_fraction：高表达基因的阈值
- key_added ：是否再obs里面加一个属性
- layer：针对哪一个layer
案例
from anndata import AnnData
import scanpy as sc
sc.settings.verbosity = 2
np.set_printoptions(precision=2)
adata = AnnData(np.array([
[3, 3, 3, 6, 6],
[1, 1, 1, 2, 2],
[1, 22, 1, 2, 2],
]))
adata.X
array([[ 3., 3., 3., 6., 6.],
[ 1., 1., 1., 2., 2.],
[ 1., 22., 1., 2., 2.]], dtype=float32)
X_norm = sc.pp.normalize_total(adata, target_sum=1, inplace=False)['X']
X_norm
array([[0.14, 0.14, 0.14, 0.29, 0.29],
[0.14, 0.14, 0.14, 0.29, 0.29],
[0.04, 0.79, 0.04, 0.07, 0.07]], dtype=float32)
X_norm = sc.pp.normalize_total(
adata, target_sum=1, exclude_highly_expressed=True,
max_fraction=0.2, inplace=False
)['X']
The following highly-expressed genes are not considered during normalization factor computation:
['1', '3', '4']
X_norm
array([[ 0.5, 0.5, 0.5, 1. , 1. ],
[ 0.5, 0.5, 0.5, 1. , 1. ],
[ 0.5, 11. , 0.5, 1. , 1. ]], dtype=float32)
代码解读
对特定的代码进行重点介绍一下，有以下三个：
对于高表达基因的确定
if exclude_highly_expressed:
counts_per_cell = adata.X.sum(1) # original counts per cell
counts_per_cell = np.ravel(counts_per_cell)
# at least one cell as more than max_fraction of counts per cell
gene_subset = (adata.X > counts_per_cell[:, None] * max_fraction).sum(0)
gene_subset = np.ravel(gene_subset) == 0
msg += (
' The following highly-expressed genes are not considered during '
f'normalization factor computation:\n{adata.var_names[~gene_subset].tolist()}'
)
确定size factor
counts_per_cell = X.sum(1)
counts_per_cell = np.ravel(counts_per_cell).copy()
adata.X = _normalize_data(adata.X, counts_per_cell, target_sum)
标准化
def _normalize_data(X, counts, after=None, copy=False):
X = X.copy() if copy else X
if issubclass(X.dtype.type, (int, np.integer)):
X = X.astype(np.float32) # TODO: Check if float64 should be used
counts = np.asarray(counts) # dask doesn't do medians
after = np.median(counts[counts > 0], axis=0) if after is None else after
counts += counts == 0
counts = counts / after
if issparse(X):
sparsefuncs.inplace_row_scale(X, 1 / counts)
else:
np.divide(X, counts[:, None], out=X)
return X
相信大家看了代码，能够理解内部的运行方式，请继续关注我们吧。
参考资料
https://scanpy.readthedocs.io/en/stable/generated/scanpy.pp.normalize_total.html
