tooluniverse-single-cell

单细胞基因组学与表达矩阵分析

规则零——先检查是否已有预计算结果

在执行以下任何指令之前，先扫描数据文件夹，查找：

• *_executed.ipynb → 使用 tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}' 读取，并将其单元格输出作为权威答案引用
• 预计算结果文件（名称类似 *results*、*deseq*、*enrich*、*stats*、*_simplified.csv 的 CSV/TSV 文件）→ 直接读取并报告所需值
• 规范分析脚本（analysis.R、run_*.py、find_*.R、*.Rmd）→ 按原样执行并读取输出

仅当以上内容均不存在时，才按照本技能下方的重新分析流程操作。从原始数据重新运行会产生不同于已发布答案的数值，而且耗时更长（通常需要多出 5–10 倍的轮次）。

---

使用 scanpy、anndata、scipy 和 ToolUniverse 进行全面的单细胞 RNA-seq 分析与表达矩阵处理。

---

查证，不要猜测

不确定科学事实时，应先搜索数据库（PubMed、UniProt、ChEMBL、ClinVar 等），而不是凭记忆推断。经数据库核实的答案总是比猜测更可靠。

---

何时使用本技能

适用于用户：

• 拥有 scRNA-seq 数据（h5ad、10X、CSV 计数矩阵）并希望进行分析
• 需要 scRNA-seq 质量控制 / QC 门控：根据线粒体比例（pct_counts_mt）、基因/UMI 计数、双细胞、环境 RNA、空液滴确定细胞筛选条件
• 询问细胞类型鉴定、聚类或注释
• 需要按细胞类型或条件进行差异表达分析
• 希望分析基因表达相关性（例如，按细胞类型分析基因长度与表达量的关系）
• 询问表达数据的 PCA、UMAP、t-SNE
• 需要对表达矩阵进行 Leiden/Louvain 聚类
• 希望比较细胞类型之间的统计差异（t 检验、ANOVA、倍数变化）
• 询问标记基因、批次校正、轨迹或细胞间通信

不适用于（请改用其他技能）：

• 仅进行 Bulk RNA-seq DESeq2 分析 -> tooluniverse-rnaseq-deseq2
• 仅进行基因富集分析 -> tooluniverse-gene-enrichment
• VCF/变异分析 -> tooluniverse-variant-analysis
• 深度学习整合/批次校正（scVI、scANVI、totalVI、PeakVI、MultiVI、DestVI、veloVI）-> tooluniverse-scvi-tools

---

核心原则

1. 数据优先——分析前加载、检查并验证数据
2. 以 AnnData 为中心——所有数据都通过 anndata 对象流转
3. 关注细胞类型——需要时按细胞类型进行子集筛选
4. 统计严谨——进行归一化、多重检验校正并计算效应量
5. 以问题为导向——解析用户实际提出的问题

---

必需软件包

import scanpy as sc, anndata as ad, pandas as pd, numpy as np
from scipy import stats
from scipy.cluster.hierarchy import linkage, fcluster
from sklearn.decomposition import PCA
from statsmodels.stats.multitest import multipletests
import gseapy as gp  # enrichment
import harmonypy     # batch correction (optional)

安装：pip install scanpy anndata leidenalg umap-learn harmonypy gseapy pandas numpy scipy scikit-learn statsmodels

---

工作流决策树

START: User question about scRNA-seq data
|
+-- FULL PIPELINE (raw counts -> annotated clusters)
|   Workflow: QC -> Normalize -> HVG -> PCA -> Cluster -> Annotate -> DE
|   See: references/scanpy_workflow.md
|
+-- DIFFERENTIAL EXPRESSION (per-cell-type comparison)
|   Most common pattern: per-cell-type DE
|   See: analysis_patterns.md "Pattern 1"
|
+-- CORRELATION ANALYSIS (gene property vs expression)
|   Pattern: Gene length vs expression correlation
|   See: analysis_patterns.md "Pattern 2"
|
+-- CLUSTERING & PCA (expression matrix analysis)
|   See: references/clustering_guide.md
|
+-- CELL COMMUNICATION (ligand-receptor interactions)
|   See: references/cell_communication.md
|
+-- TRAJECTORY ANALYSIS (pseudotime)
    See: references/trajectory_analysis.md

数据格式处理：

• h5ad -> sc.read_h5ad()
• 10X -> sc.read_10x_mtx() 或 sc.read_10x_h5()
• CSV/TSV -> pd.read_csv() -> 转换为 AnnData（检查方向！）

---

数据加载

AnnData 要求：细胞/样本为行（obs），基因为列（var）

adata = sc.read_h5ad("data.h5ad")  # h5ad already oriented

# CSV/TSV: check orientation
df = pd.read_csv("counts.csv", index_col=0)
if df.shape[0] > df.shape[1] * 5:  # genes > samples by 5x => transpose
    df = df.T
adata = ad.AnnData(df)

# Load metadata
meta = pd.read_csv("metadata.csv", index_col=0)
common = adata.obs_names.intersection(meta.index)
adata = adata[common].copy()
for col in meta.columns:
    adata.obs[col] = meta.loc[common, col]

---

质量控制与 QC 门控（在下游分析前执行）

QC 门控用于在归一化、聚类或 DE 之前确定哪些细胞和基因是真实的。跳过或草率处理会造成隐性影响：双细胞会变成虚假的“中间”状态，环境 RNA 会将标记基因扩散到各个簇中，空液滴则会抬高细胞计数。应用门控前，绝不要报告过滤后的细胞数；也不要报告实际未运行过的截断值。

如实执行：QC 通过 Bash/Python 运行 scanpy/AnnData。如果未安装 scanpy，不要编造指标——打印安装计划并停止： python scripts/scrna_qc.py --install-plan（以状态码 0 退出，列出缺失的软件包，并建议运行 pip install scanpy anndata scrublet）。

每个细胞的 QC 指标（要计算的内容）

vn = adata.var_names.str.upper()
adata.var['mt']   = vn.str.startswith('MT-')            # mitochondrial
adata.var['ribo'] = vn.str.startswith(('RPS', 'RPL'))   # ribosomal protein
adata.var['hb']   = vn.str.contains(r'^HB[^P]', regex=True)  # hemoglobin (RBC)
sc.pp.calculate_qc_metrics(
    adata, qc_vars=['mt', 'ribo', 'hb'],
    percent_top=None,   # REQUIRED for small gene panels (<500) — else IndexError
    log1p=True, inplace=True)

adata.obs 中的关键指标：n_genes_by_counts、total_counts、 pct_counts_mt、pct_counts_ribo、pct_counts_hb。

每项指标为何能标示问题（从生物学原理出发，而非依赖神奇数字）

• pct_counts_mt 高 -> 细胞正在死亡或受到应激。 细胞膜破裂会让胞质 mRNA 泄漏，而线粒体转录本仍被困在细胞内，导致捕获的 RNA 中线粒体 RNA 占比升高。阈值取决于组织（心肌细胞／肝细胞在基线状态下线粒体含量就很高——一刀切地设为 10% 会误删健康细胞）。
• n_genes_by_counts／total_counts 低 -> 空液滴或碎片 （捕获到的只有环境 RNA；基因数少、测序深度低）。
• 计数／基因数极高 -> 双细胞（两个转录组约为两倍测序深度和多样性）——但仅凭计数高判断的可靠性较低；应使用双细胞检测工具（见下文）。
• pct_counts_hb 高 -> 实体组织受到红细胞／血液污染。

选择阈值——根据分布调整（MAD），而非硬编码

硬编码阈值（mt<5%、n_genes<2500）只能作为起点；它们不适用于富含线粒体的组织，也不适用于测序深度不同的文库。优先使用基于 MAD 的规则（对要去除的极端离群值具有稳健性）：标记偏离中位数超过 nmads 个中位数绝对偏差的细胞。在经过 log1p 转换的计数／基因数上使用 nmads=5（上下两端），在 pct_counts_mt 上仅对上端使用 nmads=3，并为线粒体指标设置一个符合生物学特征的上限，避免整体退化的样本通过筛选。 务必先可视化分布（小提琴图 + total_counts 与 pct_counts_mt 的散点图——濒死细胞位于低计数／高线粒体占比区域）。

运行辅助工具（计算指标 + 基于 MAD 的筛选，并报告每一步移除的细胞数）：

python scripts/scrna_qc.py data.h5ad --doublets        # or --install-plan first

双细胞、环境 RNA、空液滴（单细胞指标无法检测这些问题）

• 双细胞：使用 Scrublet（sc.pp.scrublet，scanpy >=1.10）或 scDblFinder（R）。 合并样本前，按样本分别运行；标记后按簇剔除（双细胞会形成桥接簇）。expected_doublet_rate 约为每回收 1,000 个细胞中 0.8%（10x）。
• 环境 RNA：每个液滴中都有无细胞的“汤”状 mRNA——应进行计数校正 （SoupX／DecontX，R），而非过滤细胞。单细胞 QC 无法检测环境 RNA；如果标记基因呈现不合常理的普遍表达，应怀疑有环境 RNA。
• 空液滴：应在单细胞 QC 之前处理。CellRanger 的过滤矩阵已应用类似 EmptyDrops 的判定；如果只有原始矩阵，应在单细胞 QC 前运行 EmptyDrops（DropletUtils，R）或通过条形码排序曲线的拐点进行判定。

QC 解读表（指标 -> 关注问题 -> 处理方式）

QC 指标	方向	典型问题	建议操作
pct_counts_mt	高	细胞正在死亡或受到应激（细胞膜破裂）	过滤（MAD 上限 + 考虑组织特征的上限；富含线粒体的组织应提高上限）
n_genes_by_counts	极低	空液滴／碎片	过滤（min_genes 约为 200 + 低端 MAD）
total_counts	极低	测序深度不足／捕获失败	过滤（低端 MAD；检查条形码排序曲线的拐点）
n_genes_by_counts／total_counts	极高	双细胞（一个液滴中有两个细胞）	标记，运行 Scrublet/scDblFinder，然后剔除——不要仅凭计数设置硬上限
predicted_doublet（Scrublet）	True	多细胞液滴	按样本过滤；不确定时先分簇再剔除
pct_counts_hb	高	红细胞／血液污染	非血液组织中过滤；血液样本中调查
pct_counts_ribo	极高／极低	低复杂度／应激，或细胞类型信号	标记／调查（核糖体含量因细胞类型而异；很少作为硬性过滤条件）
各簇中标记基因普遍表达	—	环境 RNA 污染	调查——运行 SoupX/DecontX，不要在未处理的情况下继续
细胞很多，但计数中位数低	—	空液滴未被去除	调查——应用 EmptyDrops／拐点法，然后重新过滤

完整原理、MAD 代码、Scrublet/SoupX/EmptyDrops 示例及操作顺序：references/scrna_qc.md。内联流程 QC：references/scanpy_workflow.md 第 2 阶段。

完整流程（快速参考）

import scanpy as sc

adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")

# QC
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)
adata = adata[adata.obs['pct_counts_mt'] < 20].copy()
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)

# Normalize + HVG + PCA
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata.copy()
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
sc.tl.pca(adata, n_comps=50)

# Cluster + UMAP
sc.pp.neighbors(adata, n_pcs=30)
sc.tl.leiden(adata, resolution=0.5)
sc.tl.umap(adata)

# Find markers + Annotate + Per-cell-type DE
sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon')

---

差异表达决策树

Single-Cell DE (many cells per condition):
  Use: sc.tl.rank_genes_groups(), methods: wilcoxon, t-test, logreg
  Best for: Per-cell-type DE, marker gene finding

Pseudo-Bulk DE (aggregate counts by sample):
  Use: R DESeq2 via `tu run run_deseq2_analysis` or Rscript (NOT pydeseq2 — gives different DEG counts)
  Best for: Sample-level comparisons with replicates

Statistical Tests Only:
  Use: scipy.stats (ttest_ind, f_oneway, pearsonr)
  Best for: Correlation, ANOVA, t-tests on summaries

---

统计检验（快速参考）

from scipy import stats
from statsmodels.stats.multitest import multipletests

# Pearson/Spearman correlation
r, p = stats.pearsonr(gene_lengths, mean_expression)

# Welch's t-test
t_stat, p_val = stats.ttest_ind(group1, group2, equal_var=False)

# ANOVA
f_stat, p_val = stats.f_oneway(group1, group2, group3)

# Multiple testing correction (BH)
reject, pvals_adj, _, _ = multipletests(pvals, method='fdr_bh')

---

批次校正（Harmony）

import harmonypy
sc.tl.pca(adata, n_comps=50)
ho = harmonypy.run_harmony(adata.obsm['X_pca'][:, :30], adata.obs, 'batch', random_state=0)
adata.obsm['X_pca_harmony'] = ho.Z_corr.T
sc.pp.neighbors(adata, use_rep='X_pca_harmony')
sc.tl.leiden(adata, resolution=0.5)
sc.tl.umap(adata)

---

ToolUniverse 集成

数据发现（分析前）

• CxGDisc_search_datasets：按疾病、组织、物种搜索 CELLxGENE Discover 中的 scRNA-seq 数据集。使用宽泛的疾病术语（例如，“breast cancer”而非“triple-negative”）。
• GEO_search_rnaseq_datasets / geo_search_datasets：搜索 GEO 中的 scRNA-seq 研究
• NCBI_SRA_search_runs：搜索 SRA 中的测序运行（query="single cell RNA-seq [disease]"）
• OmicsDI_search_datasets：跨存储库搜索数据集
• SCXA_list_experiments / SCXA_search_gene：EBI Single Cell Expression Atlas——380 多个经过整理的 scRNA-seq 数据集，可按物种/关键词搜索，也可搜索哪些实验中表达了某个基因。SCXA_search_gene 可能较慢（实测每次调用需 10-40 秒以上）——不要仅仅因为它比其他查询耗时更长，就认为调用失败。
• SCXA_get_cluster_marker_genes：获取特定 SCXA 实验中预先计算的、按簇/细胞类型划分的标记基因。marker_type 仅接受 "clusters" 或 "cell_types"（不接受 "cell_type"）；此外，marker_type="cell_types" 还必须提供 organism_part（例如 "pancreas"、"lung"）——实测确认，这两种情况都会触发参数验证错误，而非一般警告。
• SCP_search_studies / SCP_list_studies / SCP_get_study：Broad Institute Single Cell Portal——1000 多项公开研究，包含 8000 多万个细胞。SCP_list_studies 按细胞数量降序排序，并可通过 min_cells 筛选。
• UCSCCellBrowser_search_datasets / UCSCCellBrowser_get_dataset / UCSCCellBrowser_list_facets：300 多个经过整理的单细胞数据集，可按物种/身体部位/疾病筛选。筛选 search_datasets 前，先调用 UCSCCellBrowser_list_facets 获取有效的筛选值（facets 为 "organisms"、"body_parts"、"diseases"）。

细胞类型标记

• CellMarker_search_by_cell_type：组织特异性细胞标记（先使用 CellMarker_list_cell_types——必须使用准确名称，例如 "Regulatory T(Treg) cell"，而不是 "Regulatory T cell"）
• CellMarker_search_cancer_markers：具有实验依据的癌症背景标记
• CellMarker_search_by_gene：反向查询——哪些细胞类型表达某个基因？
• HPA_search_genes_by_query：细胞类型标记基因搜索
• PanglaoDB_list_cell_types / PanglaoDB_markers_for_cell_type / PanglaoDB_cell_types_for_gene：PanglaoDB 整理的标记基因，每个都带有 canonical_marker 标志，以及按物种划分的 sensitivity_*/specificity_* 分数。先调用 PanglaoDB_list_cell_types（可按 organ 筛选）获取准确有效的 cell_type 字符串。PanglaoDB_cell_types_for_gene 的结果不按特异性排序——经实时验证，查询 CD3E（一种经典 T 细胞标记）时，返回的第一条结果是 “Macrophages”，其 canonical_marker: true，而 specificity_human: 0.1 较低；始终逐行检查 canonical_marker 和 specificity_*/sensitivity_* 字段，不要依赖结果顺序；对注释至关重要的查询，还要与上面的 CellMarker/HPA 交叉核对。
• CellTypist_search_models / CellTypist_get_model：60 多个预训练 CellTypist 分类器的目录，用于自动细胞类型注释（与上面的标记基因查询不同）——按 keyword（组织/状况，例如 "lung"、"immune"、"fetal"）和/或 min_celltypes 筛选，或通过准确的 filename 获取某个模型的完整元数据（例如 "Immune_All_Low.pkl"，这是 CellTypist 的默认模型，包含 98 种细胞类型）。每个条目都提供 n_celltypes、source_doi 和序列化模型的 download_url——这些仅是目录元数据，并非注释结果。若要实际运行所选模型来处理表达矩阵，请使用远程工具 run_celltypist_annotate。当问题是“组织 X 有哪些预训练分类器？”而非“哪些基因标记细胞类型 Y？”时，请使用此工具（后者请使用上面的 PanglaoDB/CellMarker）。

基因注释

• MyGene_query_genes / MyGene_batch_query：基因 ID 转换
• ensembl_lookup_gene：Ensembl 基因详情
• UniProt_get_function_by_accession：蛋白质功能

细胞间通信

• OmniPath_get_ligand_receptor_interactions：配体-受体（L-R）对（CellPhoneDB、CellChatDB）
• OmniPath_get_signaling_interactions：下游信号传导
• OmniPath_get_complexes：多亚基受体

富集分析（差异表达后）

• PANTHER_enrichment：GO 富集（BP、MF、CC）
• STRING_functional_enrichment：基于网络的富集分析
• ReactomeAnalysis_pathway_enrichment：Reactome 通路

临床背景（用于肿瘤免疫学）

• DGIdb_get_drug_gene_interactions：免疫检查点靶点的药物相互作用（PD-L1 对应 genes=["CD274"]）
• civic_search_evidence_items：突变/生物标志物的临床证据
• TIMER2_immune_estimation：TCGA 免疫浸润相关性
• search_clinical_trials：临床试验匹配
• GTEx_get_expression_summary：正常组织基线表达
• PubMed_search_articles：文献背景

---

Scanpy 与 Seurat 的对应关系

操作	Seurat (R)	Scanpy (Python)
加载数据	Read10X()	sc.read_10x_mtx()
归一化	NormalizeData()	sc.pp.normalize_total() + sc.pp.log1p()
查找高变基因	FindVariableFeatures()	sc.pp.highly_variable_genes()
PCA	RunPCA()	sc.tl.pca()
聚类	FindClusters()	sc.tl.leiden()
UMAP	RunUMAP()	sc.tl.umap()
查找标记	FindMarkers()	sc.tl.rank_genes_groups()
批次校正	RunHarmony()	harmonypy.run_harmony()

---

结果解读推理框架

证据分级

等级	标准	示例
高置信度	标记基因 padj < 0.01，log2FC > 1，且在 > 25% 的簇内细胞中表达	CD3D 作为 T 细胞标记基因，padj = 1e-50，log2FC = 3.2，pct = 0.85
中等置信度	padj < 0.05，log2FC > 0.5，或在 10-25% 的簇内细胞中表达	FOXP3 在 Treg 簇中，padj = 0.001，pct = 0.18
低置信度	padj < 0.05，但 log2FC < 0.5，或簇间 pct_diff 较低	普遍表达的基因，仅有边缘性富集
不可靠	簇内细胞少于 20 个，或 QC 指标表明存在双细胞	平均 nGenes > 6000 且双细胞评分较高的簇

解读指南

• QC 指标阈值：标准筛选条件为 nGenes > 200（去除空液滴）、nGenes < 5000-6000（去除双细胞）、pct_counts_mt < 20%（去除濒死细胞）。这些阈值取决于组织类型：免疫细胞适用更严格的 nGene 筛选条件；神经元的线粒体含量可能自然较高。设置阈值前，务必先可视化分布。
• 簇分辨率指南：Leiden 分辨率 0.3-0.5 可得到较宽泛的细胞类型（T 细胞、B 细胞、髓系细胞）。分辨率 0.8-1.2 可解析亚型（CD4 初始 T 细胞、CD4 记忆 T 细胞、Treg）。分辨率 > 2.0 可能导致过度聚类（将生物学上同质的群体拆分）。通过检查每个簇是否有独特的标记基因来验证。
• 标记基因置信度等级：强标记基因具有高度特异性（簇与其余细胞之间的 pct_diff 较高）且表达量高（log2FC 较高）。在许多簇中表达、倍数变化较小的基因不是好的标记基因。请与 CellMarker 或 HPA 数据库中的已知标记基因交叉比对。
• 伪 bulk 与单细胞差异表达：比较不同条件（处理组与对照组）时，伪 bulk 差异表达分析（按样本聚合后使用 DESeq2）在统计学上比单细胞差异表达分析更可靠；后者会因同一样本中的细胞并非相互独立而夸大显著性。
• 批次效应：如果样本在 UMAP 上按批次而非生物学特征聚类，请在进行生物学解读前应用 Harmony 或其他校正方法。

综合问题

1. 基于经典标记基因，识别出的簇是否对应已知细胞类型？是否有些簇缺乏明确的生物学身份（可能是双细胞或低质量细胞）？
2. 在所选的聚类分辨率下，降低分辨率时是否有簇会合并，表明它们可能是被技术噪声拆分的同一种细胞类型？
3. 对不同条件进行差异表达分析时，单细胞与伪 bulk 方法的结果是否一致？排名靠前的差异表达基因是否具有已知的生物学意义？
4. 被 QC 标记的细胞（线粒体含量高、基因计数极端）是否集中在特定簇中？移除这些细胞后，聚类结构是否发生变化？
5. 如果应用了批次校正，校正后的簇是否仍保持预期的细胞类型特异性标记基因表达？

---

故障排除

问题	解决方案
ModuleNotFoundError: leidenalg	pip install leidenalg
稀疏矩阵错误	.toarray()：X = adata.X.toarray() if issparse(adata.X) else adata.X
矩阵方向错误	基因数多于样本数？转置
相关性中出现 NaN	筛选：valid = ~np.isnan(x) & ~np.isnan(y)
用于 DE 的细胞太少	每种细胞类型、每个条件至少需要 3 个细胞
内存错误	使用 sc.pp.highly_variable_genes() 减少特征数

---

参考文档

详细分析模式：analysis_patterns.md（按细胞类型进行差异表达分析、相关性分析、PCA、ANOVA、细胞通讯）

核心工作流：

• references/scrna_qc.md - scRNA-seq QC 门控（线粒体比例、双细胞、环境 RNA、空液滴、MAD 阈值）
• references/scanpy_workflow.md - 完整的 scanpy 流程
• references/seurat_workflow.md - Seurat 到 Scanpy 的转换
• references/clustering_guide.md - 聚类方法
• references/marker_identification.md - 标记基因、注释
• references/trajectory_analysis.md - 伪时间
• references/cell_communication.md - OmniPath/CellPhoneDB 工作流
• references/troubleshooting.md - 详细的错误解决方案

---

分析规范

DESeq2 库的选择：与权威脚本保持一致

如果数据文件夹中包含权威脚本（run_*.py、analysis.R），请使用该脚本所用的 DESeq2 库（pydeseq2 或 R DESeq2）。这两个库在相同阈值下得出的差异表达基因数量略有不同（约 2-10%），因此保持一致很重要。如果没有脚本，优先通过 run_deseq2_analysis 工具或 Rscript 使用 R DESeq2：

tu run run_deseq2_analysis '{"operation":"deseq2","counts_file":"pseudo_bulk_counts.csv","metadata_file":"sample_meta.csv","design":"~ sex","contrast":"sex, M, F","lfc_shrinkage":true}'