作者：童蒙
編輯：angelica
scanpy代碼解讀來啦~
單細(xì)胞分析第一步是對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化缓升，標(biāo)準(zhǔn)化的方法有很多帅矗，下面給大家解讀一下scanpy的一個：函數(shù)為：scanpy.pp.normalize_total
作用
對每個細(xì)胞的count進(jìn)行標(biāo)準(zhǔn)化，使得每個細(xì)胞標(biāo)準(zhǔn)化后有相同的count
- 如果使用target_sum=1e6，那么相當(dāng)于CPM標(biāo)準(zhǔn)化搂捧。
- 設(shè)置exclude_highly_expressed=True 后，非常高的表達(dá)基因會被排除到計算size factor中懂缕，這個會影響其他的基因允跑。需要同max_fraction進(jìn)行連用，只要有一個細(xì)胞里面超過了，就會被判別為高表達(dá)基因聋丝。默認(rèn)為0.05索烹。
常見參數(shù)
- adata：內(nèi)置的AnnDta數(shù)據(jù)
- target_sum: 如果設(shè)置為none，那么會用所有樣品的median值來替代弱睦；
- exclude_highly_expressed ：是否去除高表達(dá)基因
- max_fraction：高表達(dá)基因的閾值
- 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)
代碼解讀
對特定的代碼進(jìn)行重點介紹一下百姓，有以下三個：
對于高表達(dá)基因的確定
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)
標(biāo)準(zhǔn)化
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
相信大家看了代碼，能夠理解內(nèi)部的運行方式况木，請繼續(xù)關(guān)注我們吧垒拢。
參考資料
https://scanpy.readthedocs.io/en/stable/generated/scanpy.pp.normalize_total.html
