Scanpy 进阶高级分析 — 单细胞数据深度挖掘工具包¶

# 创建专用 conda 环境（Scanpy 1.12 要求 Python >=3.12）
conda create -n scanpy_env python=3.12 -y

# 激活环境
conda activate scanpy_env

# 安装 Scanpy 及常用依赖（1.12.1 为最新稳定版）
pip install scanpy==1.12.1

# 安装高性能计算支持包
pip install leidenalg igraph numba

# 安装批次校正工具 harmonypy
pip install harmonypy

# 安装 GPU 加速支持（可选，需 CUDA 环境）
pip install rapids-singlecell

import scanpy as sc                    # 导入 Scanpy 主包
import numpy as np                     # 数值计算
import pandas as pd                    # 数据框操作

# 设置全局随机种子保证结果可重复
sc.settings.seed = 42
# 设置图片分辨率，适合发论文
sc.settings.set_figure_params(dpi=150, facecolor='white')

# ── 数据读入 ──────────────────────────────────────────
adata = sc.read_h5ad('processed.h5ad')          # 读取已处理的 h5ad 文件
# 或从 10X Cellranger 结果读入
adata = sc.read_10x_mtx('path/to/filtered_feature_bc_matrix/')

# ── 基础质控 ──────────────────────────────────────────
# 计算每个细胞的线粒体基因占比（mt- 开头）
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)

# 过滤低质量细胞（线粒体基因 <20%，基因数 200-6000）
sc.pp.filter_cells(adata, min_genes=200)        # 每细胞至少 200 个基因
sc.pp.filter_genes(adata, min_cells=3)          # 每基因至少 3 个细胞表达
adata = adata[adata.obs.pct_counts_mt < 20, :]  # 去除高线粒体比例细胞

# ── 标准化与降维 ──────────────────────────────────────
sc.pp.normalize_total(adata, target_sum=1e4)    # 每细胞标准化到总量 10000
sc.pp.log1p(adata)                              # log(x+1) 变换，稳定方差
sc.pp.highly_variable_genes(adata, n_top_genes=3000)  # 选高变基因
sc.pp.scale(adata, max_value=10)                # 标准化到均值0方差1，最大截断10
sc.tl.pca(adata, n_comps=50)                    # PCA 降维到 50 个主成分

import scanpy as sc
import matplotlib.pyplot as plt

# ── 完整进阶分析流程 ──────────────────────────────────

# 1. 读取并预处理数据
adata = sc.read_h5ad('raw_counts.h5ad')
sc.pp.filter_cells(adata, min_genes=500)        # 至少 500 个基因
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# 2. 保存原始计数矩阵（差异表达时需要）
adata.layers['counts'] = adata.X.copy()

# 3. 批次校正（使用 Harmony，适合多样本整合）
sc.pp.highly_variable_genes(adata, n_top_genes=3000, batch_key='batch')
sc.tl.pca(adata)
sc.external.pp.harmony_integrate(adata, 'batch')  # Harmony 批次校正
sc.pp.neighbors(adata, use_rep='X_pca_harmony')   # 用校正后的 PCA 建图

# 4. 聚类与可视化
sc.tl.umap(adata)                                  # UMAP 非线性降维
sc.tl.leiden(adata, resolution=0.6)                # Leiden 算法聚类
sc.pl.umap(adata, color=['leiden', 'batch'])        # 可视化聚类和批次

# 5. 细胞周期打分（过滤细胞周期效应）
s_genes = sc.queries.mitochondrial_genes('hsapiens')  # 占位，实际用 S 期基因列表
# 使用内置的细胞周期基因列表
cell_cycle_genes = [x.strip() for x in open('cell_cycle_genes.txt')]
s_genes = cell_cycle_genes[:43]                    # S 期基因
g2m_genes = cell_cycle_genes[43:]                  # G2M 期基因
sc.tl.score_genes_cell_cycle(adata, s_genes=s_genes, g2m_genes=g2m_genes)
sc.pl.umap(adata, color=['phase'])                  # 可视化细胞周期阶段

# 6. 差异表达分析（Wilcoxon 秩和检验）
sc.tl.rank_genes_groups(
    adata,
    groupby='leiden',                              # 按照聚类分组
    method='wilcoxon',                             # 推荐使用 Wilcoxon 检验
    layer='counts',                                # 用原始计数
    pts=True                                       # 计算每个基因的细胞表达比例
)
sc.pl.rank_genes_groups_dotplot(
    adata, n_genes=5, groupby='leiden'             # 每个 cluster 展示前5个 marker
)

# 7. 自动注释细胞类型（使用 CellTypist，见独立教程）
# 手动注释示例
new_cluster_names = {
    '0': 'CD4 T cells', '1': 'CD8 T cells',
    '2': 'B cells',     '3': 'NK cells',
    '4': 'Monocytes',   '5': 'DCs'
}
adata.obs['cell_type'] = adata.obs['leiden'].map(new_cluster_names)
sc.pl.umap(adata, color='cell_type', legend_loc='on data')

# 8. 保存结果
adata.write_h5ad('final_annotated.h5ad')
print("分析完成！")

# 安装
pip install scanpy==1.12.1 leidenalg igraph

# 常用一行命令
sc.pp.recipe_zheng17(adata)         # 标准化流程一键搞定
sc.tl.leiden(adata, resolution=0.8) # 聚类
sc.pl.umap(adata, color='leiden')   # UMAP 可视化

# 输出文件格式
adata.write_h5ad('result.h5ad')     # 保存为 h5ad 格式（推荐）
adata.write('result.loom')          # 保存为 loom 格式