自学内容网 自学内容网
自学内容网 > 正文

Scanpy

🕗 发布于 01-11 00:15 机器学习  人工智能  

0.需求

一个可怜的医学生朋友突然开始学python，从hello world开始的那种，学了一周终于学会了print()的一百种应用。但是她导师说明天就要学会使用Scanpy做项目出图，于是我带着圣光来了。

官方教程：装置——简洁

前提：已经安装好anaconda。

1.创建环境，安装依赖

（1）创建新的conda环境，我这里起名scanpy

conda create -n scanpy

（2）激活该环境，也就是进入到这个环境

conda activate scanpy

（3）根据官网，安装依赖

conda install -c conda-forge scanpy python-igraph leidenalg

注意是在当前的 conda 环境里安装，这里涉及到三个生物信息学/单细胞分析常用的 Python 包

  • scanpy：单细胞转录组分析主力库

  • python-igraph：图论算法（Scanpy 做聚类、PAGA 时需要）

  • leidenalg：Leiden 社区检测算法（Scanpy 默认的聚类算法之一）

【非必要：拉取版本代码】

如果涉及scanpy库的实现原理这些，可以用该方法获取到代码包去了解。

git仓库在这：GitHub - scverse/scanpy: Single-cell analysis in Python. Scales to >100M cells.

【如果进不去】换个镜像会快一点，但是有可能不可用：GitHub - scverse/scanpy: Single-cell analysis in Python. Scales to >100M cells.

【如果没git】官网是用很炫酷的方法直接git clone下来，但是医学生可以用更直观的方式，直接下载，不至于再学git。直接点进去链接，手动下载。

这个下载的zip包里就是要运行的代码文件夹，解压到自己喜欢的地方，然后用编译器打开。

【如果有git】直接在目标存放代码的路径下执行git命令（使用对应的镜像）：

git clone https://github.com/scverse/scanpy.git

殊途同归，最后进入到主目录里：

cd scanpy

本地获取到整个目录结构：

2.预处理与聚类

官方教程使用notebook，那我们也创建一个新的ipynb文件。注意内核选择，即前面所创建的这个环境。接下来都按照官方教程走。

2.1 基础准备

（1）引入依赖包

# Core scverse libraries
from __future__ import annotations
import anndata as ad
# Data retrieval
import pooch
import scanpy as sc

【报错】No module named 'pooch'，说明环境里没有这个包，那就在环境里装

conda install -c conda-forge pooch

（2）配置画图默认样式

一个全局配置，把以后所有 scanpy.pl.* 出的图，一次性设成尺寸50、背景纯白。

sc.settings.set_figure_params(dpi=50, facecolor="white")

（3）下载数据

官方教程是联网下载了 figshare 上 DOI 对应的那个数据集。

EXAMPLE_DATA = pooch.create(
    path=pooch.os_cache("scverse_tutorials"),
    base_url="doi:10.6084/m9.figshare.22716739.v1/",
)
EXAMPLE_DATA.load_registry_from_doi()

【报错】国内网络连 DOI 解析问题。按照这个方法在本地下载：Scanpy教程-预处理和聚类数据集-CSDN博客
然后把下载后的2个h5文件放到缓存路径里。
【缓存路径的获取与创建】

import pooch, os
path = pooch.os_cache("scverse_tutorials")
print(path) 
os.makedirs(path, exist_ok=True)

放完两个文件之后运行以下代码，pooch会从缓存里取到数据，就不会再联网下载了。文件现在被存放在 EXAMPLE_DATA 这个变量里。

import hashlib
from pathlib import Path

files = [
    "s1d1_filtered_feature_bc_matrix.h5",
    "s1d3_filtered_feature_bc_matrix.h5",
]

registry = {}
for fname in files:
    fpath = Path(path) / fname
    if not fpath.exists():
        raise FileNotFoundError(f"文件不存在: {fpath}")
    
    # 计算 SHA256
    sha256_hash = hashlib.sha256()
    with open(fpath, "rb") as f:
        for chunk in iter(lambda: f.read(4096), b""):
            sha256_hash.update(chunk)
    hash_hex = sha256_hash.hexdigest()
    
    registry[fname] = f"sha256:{hash_hex}"
    print(f"{fname} -> {hash_hex}")

print("\n生成的 registry 字典：")
print(registry)

EXAMPLE_DATA = pooch.create(
    path=pooch.os_cache("scverse_tutorials"),
    base_url="",  # 不需要联网
    registry=registry,
)

（4）加载与查看数据

从 EXAMPLE_DATA 里逐个抓取文件读取数据。

samples = {
    "s1d1": "s1d1_filtered_feature_bc_matrix.h5",
    "s1d3": "s1d3_filtered_feature_bc_matrix.h5",
}
adatas = {}

for sample_id, filename in samples.items():
    path = EXAMPLE_DATA.fetch(filename)
    sample_adata = sc.read_10x_h5(path)
    sample_adata.var_names_make_unique()
    adatas[sample_id] = sample_adata

adata = ad.concat(adatas, label="sample")
adata.obs_names_make_unique()
print(adata.obs["sample"].value_counts())
adata

这段代码干了三件事：

  1. 把两个 10x 单细胞数据文件（.h5）读进来；

  2. 给每个细胞打上样本标签；

  3. 合并成一个大数据集，方便后续分析。

.h5 是 10x 的官方“过滤表达矩阵”，里面是一个稀疏的“细胞×基因”表。可以想象成一个细胞 × 基因的 Excel 表，每个格子是某个细胞中某个基因的UMI计数（表达量）。因为是稀疏数据，所以用 .h5 这种高效格式存。

输出说明第一个文件包含了8785个细胞数据，第二个文件包含8340个细胞数据。在合并之后，总共有17125个细胞（obs），然后每个细胞有36601个基因（var），所以形成一个17125 * 36601的矩阵。

2.2 质量控制

（1）单独标记mt、ribo、hb

查个资料，这三类基因常被用来判断细胞质量。

  • MT-：线粒体基因（人源）；

  • RPS/RPL：核糖体蛋白基因；

  • HB...：血红蛋白基因（红细胞污染指示）。

# mitochondrial genes, "MT-" for human, "Mt-" for mouse
adata.var["mt"] = adata.var_names.str.startswith("MT-")
# ribosomal genes
adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))
# hemoglobin genes
adata.var["hb"] = adata.var_names.str.contains("^HB[^(P)]")

相当于对全部的基因做字符匹配，标记符合匹配规则的基因。可以稍微看一下所存放的mt数据，看样子key为基因名，value为bool类型，即true/false，示意是否属于“mt”类型。因为总共有36601个基因，所以总长度为36601。

（2）质量控制计算

这里的输入就用到了前面得到的整个矩阵数据，和三类标记。

sc.pp.calculate_qc_metrics(adata, qc_vars=["mt", "ribo", "hb"], inplace=True, log1p=True)

涉及到知识盲区，问一问万能的大模型：

“写进obs”和"新增这些列"，说明这些qc指标是以细胞为单位的，即，对一个细胞的所有基因做统计，然后得到这个细胞的这些指标，存放到var里。

（3）小提琴图

官方教程里画出来的是这些指标中的

所以对应代码里用的这些指标默认的列名，画出小提琴图：

sc.pl.violin(
    adata,
    ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
    jitter=0.4,
    multi_panel=True,
)

一个图对应一个指标，在一个图里，一个点就是一个细胞。【所以一个图里会有17125个点】

（4）散点图

联合使用qc指标，即自定义几个指标，查看他们之间的关联。
如果是两个指标，在函数里自定义出x轴和y轴即可。
如果三个指标，就可以像官方教程里一样，用颜色来表示第三个维度。

sc.pl.scatter(adata, "total_counts", "n_genes_by_counts", color="pct_counts_mt")

好像可以看出来点结论：这个total和n呈正相关，total和pct负相关，n和pct负相关。【不太专业】

（5）过滤策略

像是一种异常值的过滤。官方教程里给的是：过滤表达少于100个基因的细胞和检测不到3个细胞的基因。这里主要就是阈值的控制。

sc.pp.filter_cells(adata, min_genes=100)
sc.pp.filter_genes(adata, min_cells=3)

运行这个代码之后，再去重新做前面的字符匹配，会发现数量减少了，因为被筛掉了一些基因：


原文地址：https://blog.csdn.net/zzzzz_12321/article/details/156766642

免责声明：本站文章内容转载自网络资源，如侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！

相关文章

自学内容网
Copyright © 2015-2026