scanpy源码浅析：或谈python面向对象结构

我们已经知道scanpy是一个由Python组织的功能强大的单细胞数据分析工具，用的顺手的同时，我们不禁要问：她是如何组织的？

于是，我们尝试阅读他的源代码。在此之前我们要知道scanpy的源码在哪里放着：

python 
Python 3.7.4 (default, Aug 13 2019, 20:35:49) 
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import scanpy as sc 
>>> sc
<module 'scanpy' from 'pathto/lib/python3.7/site-packages/scanpy/__init__.py'>

映入我们眼帘的是这样的结构，第一步就是要理解这里的文件都是什么作用及其调用关系：

.
├── api【调用接口】
├── cli.py
├── _compat.py
├── datasets 【数据集】
├── external【外部接口】
├── get.py
├── __init__.py
├── logging.py
├── __main__.py
├── neighbors【计算图结构】
├── plotting【绘图模块】
├── preprocessing【预处理模块】
├── __pycache__
├── queries【查询，富集等】
├── readwrite.py
├── _settings.py
├── sim_models
├── tools
└── _utils.py

首先我们要理解一下init.py这个文件。在scanpy中我们看到主要的面向用户的函数都在这里的，这里的 __init__.py是这样的：

# some technical stuff
import sys
from ._utils import pkg_version, check_versions, annotate_doc_types

__author__ = ', '.join([
    'Alex Wolf',
    'Philipp Angerer',
    'Fidel Ramirez',
    'Isaac Virshup',
    'Sergei Rybakov',
    'Gokcen Eraslan',
    'Tom White',
    'Malte Luecken',
    'Davide Cittaro',
    'Tobias Callies',
    'Marius Lange',
    'Andrés R. Muñoz-Rojas',
])
__email__ = ', '.join([
    'f.alex.wolf@gmx.de',
    'philipp.angerer@helmholtz-muenchen.de',
    # We don’t need all, the main authors are sufficient.
])
try:
    from setuptools_scm import get_version
    __version__ = get_version(root='..', relative_to=__file__)
    del get_version
except (LookupError, ImportError):
    __version__ = str(pkg_version(__name__))

check_versions()
del pkg_version, check_versions

# the actual API
from ._settings import settings, Verbosity  # start with settings as several tools are using it
from . import tools as tl
from . import preprocessing as pp
from . import plotting as pl
from . import datasets, logging, queries, external, get

from anndata import AnnData
from anndata import read_h5ad, read_csv, read_excel, read_hdf, read_loom, read_mtx, read_text, read_umi_tools
from .readwrite import read, read_10x_h5, read_10x_mtx, write, read_visium
from .neighbors import Neighbors

set_figure_params = settings.set_figure_params

# has to be done at the end, after everything has been imported
annotate_doc_types(sys.modules[__name__], 'scanpy')
del sys, annotate_doc_types

__init__.py该文件的作用就是相当于把自身整个文件夹当作一个包来管理，每当有外部import的时候，就会自动执行里面的函数。通常init.py 文件为空，但是我们还可以为它增加其他的功能。我们在导入一个包时，实际上是导入了它的init.py文件。这样我们可以在init.py文件中批量导入我们所需要的模块，不再需要一一导入。

api

整个scanpy中，api是其他模块的接口，也就是api中包含了其他模块的调用方式。

tree api/
api/
├── datasets.py
├── export_to.py
├── __init__.py
├── logging.py
├── pl.py
├── pp.py
├── __pycache__
│   ├── datasets.cpython-37.pyc
│   ├── export_to.cpython-37.pyc
│   ├── __init__.cpython-37.pyc
│   ├── logging.cpython-37.pyc
│   ├── pl.cpython-37.pyc
│   ├── pp.cpython-37.pyc
│   ├── queries.cpython-37.pyc
│   └── tl.cpython-37.pyc
├── queries.py
└── tl.py

1 directory, 16 files

api就像一个管理器，没有明显地定义具体的功能。具体的方法在下面的模块中定义，主要有三个核心功能：

pp 【预处理】
tl 【计算】
pl【绘图】

以及datasets、logging、queries、export_to四个不常用的功能。

datasets/
├── _ebi_expression_atlas.py
├── __init__.py
└── __pycache__
    ├── _ebi_expression_atlas.cpython-37.pyc
    └── __init__.cpython-37.pyc

1 directory, 4 files

datasets中__init__.py是这样的，调用不同的数据，这些数据不在本地第一次使用函数时会下载。

"""Builtin Datasets.
"""
from ._datasets import (
    blobs,
    burczynski06,
    krumsiek11,
    moignard15,
    paul15,
    toggleswitch,
    pbmc68k_reduced,
    pbmc3k,
    pbmc3k_processed,
    visium_sge,
)
from ._ebi_expression_atlas import ebi_expression_atlas

preprocessing中定义了预处理的相关模块：highly_variable_genes、normalization、qc等。

tree   preprocessing/ 
preprocessing/
├── _combat.py
├── _deprecated
│   ├── highly_variable_genes.py
│   ├── __init__.py
│   └── __pycache__
│       ├── highly_variable_genes.cpython-37.pyc
│       └── __init__.cpython-37.pyc
├── _distributed.py
├── _docs.py
├── _highly_variable_genes.py
├── __init__.py
├── _normalization.py
├── __pycache__
│   ├── _combat.cpython-37.pyc
│   ├── _distributed.cpython-37.pyc
│   ├── _docs.cpython-37.pyc
│   ├── _highly_variable_genes.cpython-37.pyc
│   ├── __init__.cpython-37.pyc
│   ├── _normalization.cpython-37.pyc
│   ├── _qc.cpython-37.pyc
│   ├── _recipes.cpython-37.pyc
│   ├── _simple.cpython-37.pyc
│   ├── _utils.cpython-37.pyc
│   ├── _utils.sparse_mean_var_minor_axis-49.py37m.1.nbc
│   └── _utils.sparse_mean_var_minor_axis-49.py37m.nbi
├── _qc.py
├── _recipes.py
├── _simple.py
└── _utils.py

如均一化脚本：

from typing import Optional, Union, Iterable, Dict

import numpy as np
from anndata import AnnData
from scipy.sparse import issparse
from sklearn.utils import sparsefuncs

from .. import logging as logg
from .._compat import Literal


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

在scanpy函数中，我们经常看到，这样的定义结构：

def read(
    filename: Union[Path, str],
    backed: Optional[Literal['r', 'r+']] = None,
    sheet: Optional[str] = None,
    ext: Optional[str] = None,
    delimiter: Optional[str] = None,
    first_column_names: bool = False,
    backup_url: Optional[str] = None,
    cache: bool = False,
    cache_compression: Union[Literal['gzip', 'lzf'], None, Empty] = _empty,
    **kwargs,
) -> AnnData:
    """\
    Read file and return :class:`~anndata.AnnData` object.

    To speed up reading, consider passing ``cache=True``, which creates an hdf5
    cache file.

    Parameters
    ----------
    filename
        If the filename has no file extension, it is interpreted as a key for
        generating a filename via ``sc.settings.writedir / (filename +
        sc.settings.file_format_data)``.  This is the same behavior as in
        ``sc.read(filename, ...)``.
    backed

def A() ->B:"""C"""。这是什么意思呢？

在一些Python的工程项目中，我们会看到函数参数中会有冒号，有的函数后面会跟着一个箭头，你可能会疑惑，这些都是什么东西？

其实函数参数中的冒号是参数的类型建议符，告诉程序员希望传入的实参的类型。函数后面跟着的箭头是函数返回值的类型建议符，用来说明该函数返回的值是什么类型。

更官方的解释：此为type hints，是Python 3.5新加的功能，作用如上所述，官方文档为 https://www.python.org/dev/peps/pep-0484/

值得注意的是，类型建议符并非强制规定和检查，也就是说即使传入的实际参数与建议参数不符，也不会报错。我认为类型建议符的作用更多的体现在软件工程方面：在多人合作的时候，我们对他人开发的代码并不熟悉，没有对类型的解释说明的话，往往需要花费更多的时间才能看出函数的参数和返回值是什么类型，有了说明符，可以方便程序员理解函数的输入与输出（具体涉及到的工作，比如静态分析与代码重构）。

Python函数参数中的冒号与箭头

在tools中写了scanpy的核心计算函数，比如paga、louvain、score_genes、tsne、umap等。

tree  tools/  
tools/
├── _dendrogram.py
├── _diffmap.py
├── _dpt.py
├── _draw_graph.py
├── _embedding_density.py
├── _ingest.py
├── __init__.py
├── _leiden.py
├── _louvain.py
├── _marker_gene_overlap.py
├── _paga.py
├── _pca.py
├── __pycache__
│   ├── _dendrogram.cpython-37.pyc
│   ├── _diffmap.cpython-37.pyc
│   ├── _dpt.cpython-37.pyc
│   ├── _draw_graph.cpython-37.pyc
│   ├── _embedding_density.cpython-37.pyc
│   ├── _ingest.cpython-37.pyc
│   ├── __init__.cpython-37.pyc
│   ├── _leiden.cpython-37.pyc
│   ├── _louvain.cpython-37.pyc
│   ├── _marker_gene_overlap.cpython-37.pyc
│   ├── _paga.cpython-37.pyc
│   ├── _pca.cpython-37.pyc
│   ├── _rank_genes_groups.cpython-37.pyc
│   ├── _score_genes.cpython-37.pyc
│   ├── _sim.cpython-37.pyc
│   ├── _top_genes.cpython-37.pyc
│   ├── _tsne.cpython-37.pyc
│   ├── _tsne_fix.cpython-37.pyc
│   ├── _umap.cpython-37.pyc
│   ├── _utils_clustering.cpython-37.pyc
│   └── _utils.cpython-37.pyc
├── _rank_genes_groups.py
├── _score_genes.py
├── _sim.py
├── _top_genes.py
├── _tsne_fix.py
├── _tsne.py
├── _umap.py
├── _utils_clustering.py
└── _utils.py

1 directory, 42 files

如umap的计算：

# 这里省略210行代码。

。。。代码去谈恋爱了。

趁着代码谈恋爱的阶段我们来介绍一下python中的_(单下划线)以及__（双下划线）。

单下划线:单下划线的变量是一种程序员之间美丽的约定——只要是这种变量就不要随便在类外部去访问它！！！
但是如果我们在导入模块时来看这个单下划线开头的变量，那就不一样了，在这里这种特殊名字的变量就变成了类似一种某个模块的“私有”变量，因为我们在使用from 模块名 import *语句导入模块时，这些单下划线开头的变量默认是不会被导入的，所以实际上这个单下划线对python的解释器有了影响。
双下划线：解析器自动转换为：_类名_成员名，代替原有成员，访问需要在原有成员名字前加上类名。如:Python自动将__name 解释成\ _student__name,我们可以用\ _student__name访问.
python中的单下划线,双下划线以及两端双下划线

`_ 的含义`

在python的类中没有真正的私有属性或方法，没有真正的私有化。

但为了编程的需要，我们常常需要区分私有方法和共有方法以方便管理和调用。那么在Python中如何做呢？

一般Python约定加了下划线 _ 的属性和方法为私有方法或属性，以提示该属性和方法不应在外部调用，也不会被from ModuleA import * 导入。如果真的调用了也不会出错,但不符合规范。

`__ 的含义`

Python中的__和一项称为name mangling的技术有关，name mangling (又叫name decoration命名修饰).在很多现代编程语言中,这一技术用来解决需要唯一名称而引起的问题,比如命名冲突/重载等.

Python中双下划线开头,是为了不让子类重写该属性方法.通过类的实例化时自动转换,在类中的双下划线开头的属性方法前加上”_类名”实现.

构建图结构

tree neighbors/
neighbors/
├── __init__.py
└── __pycache__
    └── __init__.cpython-37.pyc

1 directory, 2 files

绘图函数

tree  plotting/  
plotting/
├── _anndata.py
├── _docs.py
├── __init__.py
├── palettes.py
├── _preprocessing.py
├── __pycache__
│   ├── _anndata.cpython-37.pyc
│   ├── _docs.cpython-37.pyc
│   ├── __init__.cpython-37.pyc
│   ├── palettes.cpython-37.pyc
│   ├── _preprocessing.cpython-37.pyc
│   ├── _qc.cpython-37.pyc
│   ├── _rcmod.cpython-37.pyc
│   └── _utils.cpython-37.pyc
├── _qc.py
├── _rcmod.py
├── _tools
│   ├── __init__.py
│   ├── paga.py
│   ├── __pycache__
│   │   ├── __init__.cpython-37.pyc
│   │   ├── paga.cpython-37.pyc
│   │   └── scatterplots.cpython-37.pyc
│   └── scatterplots.py
└── _utils.py

3 directories, 22 files

我们看看著名的stacked_violin函数：

@_doc_params(show_save_ax=doc_show_save_ax, common_plot_args=doc_common_plot_args)
def stacked_violin(
    adata: AnnData,
    var_names: Union[_VarNames, Mapping[str, _VarNames]],
    groupby: Optional[str] = None,
    log: bool = False,
    use_raw: Optional[bool] = None,
    num_categories: int = 7,
    figsize: Optional[Tuple[float, float]] = None,
    dendrogram: Union[bool, str] = False,
    gene_symbols: Optional[str] = None,
    var_group_positions: Optional[Sequence[Tuple[int, int]]] = None,
    var_group_labels: Optional[Sequence[str]] = None,
    standard_scale: Optional[Literal['var', 'obs']] = None,
    var_group_rotation: Optional[float] = None,
    layer: Optional[str] = None,
    stripplot: bool = False,
    jitter: Union[float, bool] = False,
    size: int = 1,
    scale: Literal['area', 'count', 'width'] = 'width',
    order: Optional[Sequence[str]] = None,
    swap_axes: bool = False,
    show: Optional[bool] = None,
    save: Union[bool, str, None] = None,
    row_palette: str = 'muted',
    ax: Optional[_AxesSubplot] = None,
    **kwds,
):

这里我们插播一条python的知识点：装饰器。

装饰器就是一个函数对另一个函数的管理和持有。

装饰器背后的主要动机源自python面向对象编程，装饰器是在函数调用之上的修饰，这些修饰仅是当声明一个函数或者方法的时候，才会应用的额外调用。装饰器的语法以@开头，接着是装饰器显式的名字和可选的参数。紧跟着装饰器声明的是被修饰的函数，和修饰函数的可选参数。

装饰器实际上是一个函数，他们接受函数对象，但他们是怎么处理那些函数？

当包装一个函数的时候，你最终会调用它，最棒的是我们能够在包装的环境下在合适的时机调用它，我们在执行函数之前，可以运行那些预备代码，也可以在执行代码之后做个清理工作。所以一个装饰器函数，很可能在里面是这样一些代码：它定义了某个函数并在定义内的某处嵌入了对目标函数的调用或者至少一些引用。从本质上看，这些特征引入了java开发者称之为AOP的概念，可以考虑在装饰器中置入通用功能的代码来降低程序复杂度。

可以用装饰器来：

引入日志
增加计时逻辑来检测性能
给函数加入事物能力。

查询功能模块

tree   queries/  
queries/
├── __init__.py
├── __pycache__
│   ├── __init__.cpython-37.pyc
│   └── _queries.cpython-37.pyc
└── _queries.py

1 directory, 4 files

from ._queries import (
    biomart_annotations,
    gene_coordinates,
    mitochondrial_genes,
)  # Biomart queries
from ._queries import enrich  # gprofiler queries

scanpy也是可以做富集分析的啊。我去，python能做的都可以在python中做。

sc.queries.enrich(['Klf4', 'Pax5', 'Sox2', 'Nanog'], org="hsapiens")
Out[166]: 
   source  ...                                            parents
0    REAC  ...                               [REAC:R-HSA-1266738]
1   GO:BP  ...                           [GO:0001704, GO:0045165]
2    REAC  ...                                [REAC:R-HSA-452723]
3   GO:MF  ...                           [GO:0000981, GO:0001216]
4   GO:MF  ...                                       [GO:0003700]
5      WP  ...                                        [WP:000000]
6   GO:BP  ...                                       [GO:0019827]
7      WP  ...                                        [WP:000000]
8   GO:BP  ...                           [GO:0001708, GO:0001711]
9   GO:BP  ...                           [GO:0009653, GO:0009790]

external 外部功能模块

在这里，我们可以发挥想象力了。这里可以成为我们自己的接口，接入我们自己的分析模块。

tree  external/ 
external/
├── exporting.py
├── __init__.py
├── pl.py
├── pp
│   ├── _bbknn.py
│   ├── _dca.py
│   ├── __init__.py
│   ├── _magic.py
│   ├── _mnn_correct.py
│   └── __pycache__
│       ├── _bbknn.cpython-37.pyc
│       ├── _dca.cpython-37.pyc
│       ├── __init__.cpython-37.pyc
│       ├── _magic.cpython-37.pyc
│       └── _mnn_correct.cpython-37.pyc
├── __pycache__
│   ├── exporting.cpython-37.pyc
│   ├── __init__.cpython-37.pyc
│   └── pl.cpython-37.pyc
└── tl
    ├── __init__.py
    ├── _palantir.py
    ├── _phate.py
    ├── _phenograph.py
    ├── __pycache__
    │   ├── __init__.cpython-37.pyc
    │   ├── _palantir.cpython-37.pyc
    │   ├── _phate.cpython-37.pyc
    │   ├── _phenograph.cpython-37.pyc
    │   ├── _pypairs.cpython-37.pyc
    │   └── _trimap.cpython-37.pyc
    ├── _pypairs.py
    └── _trimap.py

5 directories, 28 files

哇，许多可以分析单细胞的工具啊。

from ._pypairs import cyclone, sandbag
from ._phate import phate
from ._phenograph import phenograph
from ._palantir import palantir
from ._trimap import trimap
from ._harmony_timeseries import harmony_timeseries
from ._sam import sam
from ._wishbone import wishbone

sim_models

tree  sim_models/  
sim_models/
├── __init__.py
└── __pycache__
    └── __init__.cpython-37.pyc

1 directory, 2 files

pycache

tree __pycache__/
__pycache__/
├── cli.cpython-37.pyc
├── _compat.cpython-37.pyc
├── get.cpython-37.pyc
├── __init__.cpython-37.pyc
├── logging.cpython-37.pyc
├── __main__.cpython-37.pyc
├── readwrite.cpython-37.pyc
├── _settings.cpython-37.pyc
└── _utils.cpython-37.pyc

0 directories, 9 files

这里再插播一条python的知识点：__pycache__的作用。

Python程序运行时不需要编译成二进制代码，而直接从源码运行程序，简单来说是，Python解释器将源码转换为字节码，然后再由解释器来执行这些字节码。

解释器的具体工作：

1、完成模块的加载和链接；
2、将源代码编译为PyCodeObject对象(即字节码)，写入内存中，供CPU读取；
3、从内存中读取并执行，结束后将PyCodeObject写回硬盘当中，也就是复制到.pyc或.pyo文件中，以保存当前目录下所有脚本的字节码文件。

之后若再次执行该脚本，它先检查【本地是否有上述字节码文件】和【该字节码文件的修改时间是否在其源文件之后】，是就直接执行，否则重复上述步骤。

那么，__pycache__文件夹的意义何在呢？
因为第一次执行代码的时候，Python解释器已经把编译的字节码放在__pycache__文件夹中，这样以后再次运行的话，如果被调用的模块未发生改变，那就直接跳过编译这一步，直接去_pycache_文件夹中去运行相关的 *.pyc 文件，大大缩短了项目运行前的准备时间。

基本事实：

阅读代码的时间远大于写代码的时间
学习编程这类实践类的知识就是基本的就是试错
不谋全局着不足以谋一隅

最后编辑于：2020.04.02 07:52:38

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 213,014评论 6赞 492
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 90,796评论 3赞 386
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 158,484评论 0赞 348
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 56,830评论 1赞 285
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 65,946评论 6赞 386
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 50,114评论 1赞 292
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 39,182评论 3赞 412
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 37,927评论 0赞 268
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 44,369评论 1赞 303
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 36,678评论 2赞 327
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 38,832评论 1赞 341
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 34,533评论 4赞 335
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 40,166评论 3赞 317
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 30,885评论 0赞 21
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,128评论 1赞 267
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 46,659评论 2赞 362
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 43,738评论 2赞 351

scanpy源码浅析：或谈python面向对象结构

_ 的含义

__ 的含义

推荐阅读更多精彩内容

`_ 的含义`

`__ 的含义`