570 Scanpy 单细胞分析从零入门¶

filtered_feature_bc_matrix/
├── barcodes.tsv.gz    # 每个细胞的条码（barcode），一行一个细胞
├── features.tsv.gz    # 基因列表（以前叫 genes.tsv），一行一个基因
└── matrix.mtx.gz      # 稀疏矩阵，记录每个细胞中每个基因的 UMI 计数

# 创建一个专门的单细胞分析环境
conda create -n scrnaseq python=3.12 -y  # 创建名为 scrnaseq 的环境，Python 3.12

# 激活环境
conda activate scrnaseq  # 切换到新环境

# 安装 scanpy（当前最新版 1.12.1）
conda install -c conda-forge scanpy=1.12.1 -y  # 从 conda-forge 频道安装

# 安装可视化和辅助依赖
conda install -c conda-forge leidenalg python-igraph -y  # Leiden 聚类算法依赖
pip install scrublet  # 双细胞检测工具（可选）

import scanpy as sc  # 导入 scanpy
print(sc.__version__)  # 打印版本号，应该显示 1.12.1

import scanpy as sc  # 导入 scanpy，习惯简写为 sc
import numpy as np   # 导入 numpy，用于数值计算
import pandas as pd  # 导入 pandas，用于表格操作

sc.settings.verbosity = 3  # 设置日志级别：3 = 输出提示信息
sc.settings.set_figure_params(dpi=80, facecolor='white')  # 设置画图参数：分辨率和背景色
sc.logging.print_header()  # 打印 scanpy 版本和依赖信息

# 方法一：直接用 scanpy 内置数据集（推荐新手用这个，自动下载）
adata = sc.datasets.pbmc3k()  # 下载并加载 PBMC 3k 数据
# adata 就是一个 AnnData 对象，包含了所有信息

# 方法二：从 Cell Ranger 输出目录读取
# adata = sc.read_10x_mtx(
#     'data/filtered_feature_bc_matrix/',  # Cell Ranger 输出目录的路径
#     var_names='gene_symbols',            # 用基因名（而不是 Ensembl ID）作为变量名
#     cache=True                           # 缓存为 h5ad 格式，下次读取更快
# )

# 方法三：读取已保存的 h5ad 文件
# adata = sc.read_h5ad('pbmc3k.h5ad')

print(adata)  # 查看数据基本信息
# 输出类似：AnnData object with n_obs × n_vars = 2700 × 32738
# n_obs = 细胞数（2700），n_vars = 基因数（32738）

# 确保基因名是唯一的（有些基因有重名）
adata.var_names_make_unique()  # 给重名基因加后缀，确保唯一

# 标记线粒体基因（线粒体基因名字以 "MT-" 开头）
adata.var['mt'] = adata.var_names.str.startswith('MT-')  # 布尔标记

# 计算 QC 指标
sc.pp.calculate_qc_metrics(
    adata,                     # 数据对象
    qc_vars=['mt'],            # 计算线粒体相关指标
    percent_top=None,          # 不计算 top N 基因的占比
    log1p=False,               # 不取对数
    inplace=True               # 直接写入 adata
)

# 画 QC 小提琴图，直观看数据质量
sc.pl.violin(
    adata,
    ['n_genes_by_counts',   # 每个细胞检测到的基因数
     'total_counts',        # 每个细胞的总 UMI 数
     'pct_counts_mt'],      # 线粒体基因占比（%）
    jitter=0.4,             # 加点抖动，看分布
    multi_panel=True        # 每个指标画一个面板
)

# 过滤低质量细胞
sc.pp.filter_cells(adata, min_genes=200)    # 去掉基因数 < 200 的细胞（可能是空液滴）
sc.pp.filter_cells(adata, max_genes=2500)   # 去掉基因数 > 2500 的细胞（可能是双细胞）
adata = adata[adata.obs.pct_counts_mt < 5, :]  # 去掉线粒体占比 > 5% 的细胞（可能是死细胞）

# 过滤低表达基因
sc.pp.filter_genes(adata, min_cells=3)  # 去掉在少于 3 个细胞中表达的基因

print(f"质控后：{adata.n_obs} 个细胞, {adata.n_vars} 个基因")
# 典型输出：质控后：2638 个细胞, 13714 个基因

# 归一化：让不同细胞之间的 UMI 总数可比
sc.pp.normalize_total(
    adata,
    target_sum=1e4  # 把每个细胞的总 UMI 数缩放到 10000
)

# 对数化：压缩数据范围，让高表达基因不会主导分析
sc.pp.log1p(adata)  # log1p = log(x + 1)，加 1 是为了避免 log(0)

# 保存归一化后的数据（后面差异表达要用原始数据）
adata.raw = adata  # .raw 保存一份完整数据的副本

# 找出高变基因（Highly Variable Genes, HVGs）
# 这些基因在不同细胞间表达差异大，最有信息量
sc.pp.highly_variable_genes(
    adata,
    min_mean=0.0125,   # 平均表达量下限
    max_mean=3,         # 平均表达量上限
    min_disp=0.5        # 离散度下限（变异程度）
)

# 画高变基因图：红色点是被选中的高变基因
sc.pl.highly_variable_genes(adata)

# 只保留高变基因用于后续分析
adata = adata[:, adata.var.highly_variable]  # 筛选列（基因维度）

print(f"高变基因数：{adata.n_vars}")
# 典型输出：高变基因数：约 1800 个

# 缩放数据（让每个基因的均值为 0，方差为 1）
sc.pp.scale(
    adata,
    max_value=10  # 把极端值截断到 10，避免异常值干扰
)

# PCA 降维：从几千个基因压缩到 50 个主成分
sc.tl.pca(
    adata,
    svd_solver='arpack',  # 使用 arpack 求解器，适合大规模稀疏矩阵
    n_comps=50             # 计算 50 个主成分
)

# 画碎石图（Elbow Plot）：帮你决定用多少个 PC
sc.pl.pca_variance_ratio(adata, log=True, n_pcs=50)
# 看图找"拐点"，一般在第 10-15 个 PC 附近曲线变平

# 计算邻域图（KNN图）：基于 PCA 结果找每个细胞的邻居
sc.pp.neighbors(
    adata,
    n_neighbors=10,  # 每个细胞找 10 个最近邻
    n_pcs=10         # 使用前 10 个主成分
)

# UMAP 降维：把高维数据压到 2D 用于可视化
sc.tl.umap(adata)

# Leiden 聚类：自动把细胞分成若干群
sc.tl.leiden(
    adata,
    resolution=0.5,  # 分辨率参数：越大分的群越多，0.5 是常用起点
    flavor='igraph',  # 使用 igraph 后端
    n_iterations=2    # 迭代次数
)

# 可视化聚类结果
sc.pl.umap(
    adata,
    color=['leiden'],  # 用聚类标签着色
    title='Leiden Clustering'  # 图标题
)

# PBMC 中常见细胞类型的 marker 基因
marker_genes = {
    'T cells CD4': ['IL7R', 'CD4'],          # CD4+ T 细胞
    'T cells CD8': ['CD8A', 'CD8B'],          # CD8+ T 细胞
    'NK cells': ['GNLY', 'NKG7'],             # 自然杀伤细胞
    'B cells': ['MS4A1', 'CD79A'],            # B 细胞（MS4A1 = CD20）
    'Monocytes CD14+': ['CD14', 'LYZ'],       # CD14+ 单核细胞
    'Monocytes FCGR3A+': ['FCGR3A', 'MS4A7'], # CD16+ 单核细胞
    'Dendritic cells': ['FCER1A', 'CST3'],    # 树突状细胞
    'Platelets': ['PPBP'],                     # 血小板
}

# 用 dotplot 展示每一群中 marker 基因的表达情况
sc.pl.dotplot(
    adata,
    marker_genes,          # marker 基因字典
    groupby='leiden',      # 按 leiden 聚类分组
    dendrogram=True        # 加上树状图
)

# 根据 dotplot 结果手动注释
# 假设 cluster 0 高表达 IL7R → 是 CD4+ T 细胞
new_cluster_names = [
    'CD4 T',          # cluster 0
    'CD14 Monocytes', # cluster 1
    'B cells',        # cluster 2
    'CD8 T',          # cluster 3
    'NK cells',       # cluster 4
    'FCGR3A Mono',    # cluster 5
    'Dendritic',      # cluster 6
    'Platelets',      # cluster 7
]

# 把聚类编号替换成细胞类型名字
adata.rename_categories('leiden', new_cluster_names)  # 重命名聚类

# 再画一次 UMAP，现在标签是细胞类型
sc.pl.umap(
    adata,
    color='leiden',
    title='Cell Types (PBMC 3k)',
    frameon=False   # 去掉边框，更好看
)

# 找每一群的 marker 基因（差异表达基因）
sc.tl.rank_genes_groups(
    adata,
    groupby='leiden',       # 按聚类/细胞类型分组
    method='wilcoxon',      # 使用 Wilcoxon 秩和检验（非参数检验，最常用）
    use_raw=True            # 使用归一化后的原始数据（.raw 里存的）
)

# 查看每群 top 5 marker 基因
sc.pl.rank_genes_groups(adata, n_genes=5, sharey=False)

# 以表格形式查看结果
result = adata.uns['rank_genes_groups']  # 结果存在 .uns 里
groups = result['names'].dtype.names     # 获取所有组名
df = pd.DataFrame({
    group + '_' + key: result[key][group]
    for group in groups
    for key in ['names', 'pvals_adj', 'logfoldchanges']
}).head(10)  # 查看前 10 个基因
print(df)

# 1. UMAP 图：展示特定基因的表达分布
sc.pl.umap(
    adata,
    color=['CD4', 'CD8A', 'MS4A1', 'CD14'],  # 想看的基因
    ncols=2,     # 每行 2 列
    cmap='Reds'  # 颜色方案：红色梯度
)

# 2. 小提琴图（Violin Plot）：看某个基因在不同群中的表达分布
sc.pl.violin(
    adata,
    ['CD4', 'CD8A', 'MS4A1'],  # 基因列表
    groupby='leiden',           # 按细胞类型分组
    rotation=45                 # 横轴标签旋转 45 度，防止重叠
)

# 3. 点图（Dot Plot）：同时展示表达量和表达比例
sc.pl.dotplot(
    adata,
    marker_genes,      # 前面定义的 marker 基因字典
    groupby='leiden',  # 按细胞类型分组
    standard_scale='var'  # 按基因标准化颜色
)

# 4. 热图（Heatmap）：展示 top marker 基因的表达模式
sc.pl.rank_genes_groups_heatmap(
    adata,
    n_genes=5,        # 每群展示 top 5 基因
    groupby='leiden',
    show_gene_labels=True  # 显示基因名
)

# 5. 保存结果
adata.write('pbmc3k_analyzed.h5ad')  # 保存为 h5ad 文件，下次直接读取

AnnData 对象结构：
┌─────────────────────────────────────────────┐
│  .X          表达矩阵（细胞 × 基因）          │  ← 就像 Excel 里的成绩表
│              shape: (n_cells, n_genes)       │
├─────────────────────────────────────────────┤
│  .obs        细胞元数据（DataFrame）           │  ← 学生信息表（姓名、班级、成绩）
│              每行一个细胞                      │
│              例如：聚类标签、QC指标             │
├─────────────────────────────────────────────┤
│  .var        基因元数据（DataFrame）           │  ← 科目信息表（科目名、是否必修）
│              每行一个基因                      │
│              例如：是否高变基因、线粒体标记      │
├─────────────────────────────────────────────┤
│  .obsm       细胞的多维嵌入（dict）             │  ← 降维结果
│              'X_pca': PCA 坐标                │
│              'X_umap': UMAP 坐标              │
├─────────────────────────────────────────────┤
│  .uns        非结构化数据（dict）               │  ← 杂物柜
│              聚类颜色、差异表达结果等            │
├─────────────────────────────────────────────┤
│  .obsp       细胞间的关系矩阵                  │  ← 邻域图
│              'connectivities': KNN 连接        │
│              'distances': 距离矩阵             │
├─────────────────────────────────────────────┤
│  .raw        归一化后的原始数据备份              │  ← 用于差异表达
└─────────────────────────────────────────────┘

adata.shape          # 查看维度：(细胞数, 基因数)
adata.obs.head()     # 查看前 5 个细胞的元数据
adata.var.head()     # 查看前 5 个基因的元数据
adata.obs_names[:5]  # 查看前 5 个细胞的 barcode
adata.var_names[:5]  # 查看前 5 个基因名
adata.obsm.keys()    # 查看有哪些降维结果
adata.uns.keys()     # 查看有哪些非结构化数据

MemoryError: Unable to allocate XX GiB
# 或者进程直接被系统杀掉

# 方法一：读取时使用稀疏矩阵（默认就是，确认一下）
import scipy.sparse as sp  # 导入 scipy 稀疏矩阵模块
print(type(adata.X))  # 应该是 csr_matrix 或 csc_matrix

# 方法二：尽早过滤基因，减少数据量
sc.pp.filter_genes(adata, min_cells=10)  # 把阈值提高

# 方法三：使用 backed 模式（不把全部数据加载到内存）
adata = sc.read_h5ad('large_data.h5ad', backed='r')  # 只读模式

KeyError: "Key 'CD4' not found in axis"

# 查看数据里基因名长什么样
print(adata.var_names[:5])  # 看看是 gene symbol 还是 Ensembl ID

# 如果是 Ensembl ID，读取时指定使用 gene_symbols
adata = sc.read_10x_mtx(
    'filtered_feature_bc_matrix/',
    var_names='gene_symbols'  # 关键参数：用基因名而不是 ID
)

ImportError: Please install the leiden algorithm: `conda install -c conda-forge leidenalg`

conda install -c conda-forge leidenalg python-igraph -y  # 安装 Leiden 算法和 igraph

_tkinter.TclError: no display name and no $DISPLAY environment variable

# 在脚本最开头（import 之前）加上这两行
import matplotlib  # 导入 matplotlib
matplotlib.use('Agg')  # 使用非交互式后端，不需要显示器
import matplotlib.pyplot as plt  # 然后再导入 pyplot

# 或者保存图片而不是显示
sc.settings.set_figure_params(dpi=150)  # 设置分辨率
sc.pl.umap(adata, color='leiden', save='_clusters.png')  # save 参数保存图片
# 图片会保存到 ./figures/ 目录下

TypeError: AnnData expects .X to be np.ndarray or sparse matrix, not <class 'xxx'>

import scipy.sparse as sp  # 导入稀疏矩阵模块

# 如果 .X 是 numpy array，转成稀疏矩阵节省内存
adata.X = sp.csr_matrix(adata.X)

# 如果 .X 是其他格式，先转成 numpy array
import numpy as np  # 导入 numpy
adata.X = np.array(adata.X)