tooluniverse-epigenomics

基因组学和表观基因组学数据处理

⚠️ 首要规则：长格式甲基化 CSV——统计行数，而非唯一位置数

当输入是长格式甲基化 CSV（每个 (sample, CpG_position) 占一行）， 例如包含 Pos, Chromosome, MethylationPercentage 等列时，“过滤时删除了多少个位点” 几乎总是指删除的行数，而非删除的唯一位置数。两种答案相差约 n_samples 倍。

问题表述	含义
“过滤时删除了多少个位点……”	删除的行数（= 样本数 × 未通过筛选的位置数）
“有多少个唯一 CpG 位点通过筛选”	唯一位置数（先按 Pos 去重，再筛选）

❌ 错误：先执行 df.drop_duplicates(["Pos"]).query("MethylationPercentage<10 or >90")，再用 len(filtered) → 统计唯一位置数（通常为 100–1500）

✅ 正确：先执行 df.query("MethylationPercentage<10 or MethylationPercentage>90")，再用 len(df) - len(filtered) → 统计行数（通常为 10k–30k）

如果数据有 1000+ 个位置 × 20+ 个样本，而你的答案小于 2000，说明你过早进行了去重。报告之前，请重新确认问题中使用的名词。

---

规则零——首先检查预计算结果

遵循以下任何说明之前，请先扫描数据文件夹，查找：

• *_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 倍的轮次）。

适用于生产环境的技能，结合 Python 计算（pandas、scipy、numpy、pysam、statsmodels）与 ToolUniverse 注释工具，用于表观基因组学分析。

查询，不要猜测

不确定科学事实时，请先搜索数据库。

适用场景

甲基化数据、ChIP-seq 峰、ATAC-seq、多组学整合、全基因组表观基因组统计。关键词：甲基化、CpG、ChIP-seq、ATAC-seq、组蛋白、染色质、表观遗传。

不适用于：RNA-seq DEG、变异检测、基因富集、蛋白质结构。

---

核心原则

1. 数据优先——先加载并检查数据，再进行分析
2. 问题驱动——提取具体数值答案
3. 坐标系统意识——跟踪基因组版本（hg19/hg38/mm10）、chr 前缀
4. 统计严谨性——进行 FDR 校正，筛选效应量
5. CpG 识别——解析 Illumina 探针 ID、基因组坐标

主脚本 — methylation_density.py（回答 CpG 密度问题时优先使用）

对于长格式甲基化 CSV（Pos, Chromosome, MethylationPercentage） 以及与之配对的染色体长度 CSV，始终先运行捆绑的脚本，再手动编写 pandas 代码。它可通过单次运行确定性地计算所有常见指标，避免因行数与位点数混淆而得出看似合理、实际错误的答案。

python skills/tooluniverse-epigenomics/scripts/methylation_density.py \
  --cpg <CpG csv> --chr-lengths <chr lengths csv> \
  --filter-meth-extremes 90 10

完整的 JSON 输出包含所有指标。选择与问题措辞相符的字段（不要选择名称相似的字段）：

问题措辞	脚本字段
“筛选时移除了多少个位点”	rows_removed
“有多少个唯一 CpG 位点通过筛选”	unique_pos_after_filter
“全基因组染色体平均密度”	density_avg_per_chr
“染色体 X 上的密度”	density_chromosome（传入 --chromosome X）
“全基因组总密度”	density_total_over_genome

两个密度数值（density_avg_per_chr 与 density_total_over_genome）通常相差约 2 倍，因为 CpG 在染色体上的分布并不均匀；问题询问其中一个时却报告另一个，是这里最常见的错误。

对于“移除了多少个位点”这类问题，长格式 CSV 中每个 CpG 位置有多行（每个样本一行），因此 rows_removed 的数值为数万，而 unique_pos_removed 的数值为数百。请根据问题选择相应的粒度。

区分“行”与“唯一位点”——甲基化 CSV 通常为长格式

CpG 甲基化 CSV 通常是每个（样本 × CpG 位点）一行，因此 len(df) >> n_unique_sites。在进行任何计算之前，先判断问题询问的是哪个维度：

问题措辞	维度	操作
“筛选时移除了多少个位点”	样本行	筛选后统计行数；不要按 Pos 去重。CSV 为长格式；此处的“位点”按行计数。用 len(df) 减去 len(df_filtered)。
“有多少个唯一 CpG 位点通过筛选”	唯一位置	按位置（或 Pos 列）去重，然后筛选
“全基因组染色体平均密度”	每条染色体的密度	计算每条染色体密度的平均值：(n_unique_per_chr / chr_length).mean()。不要使用 total_unique / total_genome——对于分布不均匀的 CpG，这会得到不同的答案（通常约为每条染色体平均值的一半）。
“染色体 X 上的密度”	单条染色体	X 上的唯一位置数 / X 的长度。注意物种是否正确——请检查问题中写的是 “Zebra Finch” 还是 “Jackdaw”。
“染色体间均匀分布的卡方检验”	每条染色体的唯一位置数	先筛选行，再按 (Chromosome, Pos) 去重，然后统计每条染色体的唯一位置数，以此进行卡方检验；期望值 = chr_length / total_length × n_unique_filtered

合理性检查：如果筛选后的计数比 GT 范围小两个数量级，你很可能在问题要求按行计数时进行了去重（或反之）。请用另一个维度重新计算并比较。

对于卡方均匀性检验，期望计数 = 染色体长度 / 基因组总长度 × 唯一位点数。卡方统计量取决于计数粒度（行数或唯一位点数）——由于 n 较大，按行进行的卡方检验会比按唯一位置进行的检验得到高得多的卡方值。

优先级：如果存在 *_executed.ipynb，请逐字读取其中的筛选代码——df[(df.MethylationPercentage > 90) | (df.MethylationPercentage < 10)]（不去重）和 df.drop_duplicates('Pos')（去重）在同一数据集上会得出截然不同的计数。

---

工作流

第 0 阶段：问题解析

识别数据文件、具体统计量、阈值和基因组版本。按关键词分类。 决策树请参阅 ANALYSIS_PROCEDURES.md。

第 1 阶段：甲基化处理

• 加载 beta/M 值矩阵（CSV/TSV/parquet/HDF5）
• 按方差、缺失率、探针类型、染色体、CpG 岛关系进行筛选
• 差异甲基化：对组间数据进行 T 检验/Wilcoxon 检验 + FDR
• 年龄相关 CpG：Pearson/Spearman 相关性 + FDR
• 染色体密度：CpG 数量 / 染色体长度

第 2 阶段：ChIP-seq 峰分析

• 加载 BED/narrowPeak/broadPeak，统一染色体格式
• 峰统计、基因注释、重叠分析（Jaccard）

第 3 阶段：ATAC-seq

• 检测 NFR（峰 <150bp），对区域进行分类

第 4 阶段：多组学整合

• 按探针-基因计算甲基化与表达的相关性（Pearson/Spearman + FDR）
• ChIP-seq + 表达：启动子峰与表达水平

第 5 阶段：临床数据

• 分析不同组学数据中的缺失数据，识别完整病例

第 6 阶段：ToolUniverse 注释

ENCODE 工具：

• ENCODE_search_rnaseq_experiments：assay_type（默认为 "total RNA-seq"；否则回退到 "polyA plus RNA-seq"）、biosample、limit
• ENCODE_search_histone_experiments：target（例如 "H3K27ac"）、cell_type/tissue/biosample、limit

GEO 工具：GEO_search_rnaseq_datasets、GEO_search_atacseq_datasets——均接受 limit 或 max_results

GTEx 工具：

• GTEx_get_median_gene_expression：gene_symbol（不能使用 Ensembl ID）
• GTEx_query_eqtl：gene_symbol、tissue_id（区分大小写且必须完全匹配，例如 "Whole_Blood"）

EWAS Catalog（已发表的表观基因组关联研究结果，而非原始甲基化数据）：EWASCatalog_search_by_cpg（cpg_id，例如 "cg00000029"）或 EWASCatalog_search_by_gene（gene_symbol）——当问题是“该 CpG/基因是否与已发表 EWAS 中的某种性状相关联”时使用，而不是用于处理已有的甲基化数据。

EWAS Atlas（独立于 EWAS Catalog 的精选 EWAS 知识库；当 EWASCatalog_* 返回 HTTP 503 时使用，该接口自 2026-10-05 起一直如此；也可用于交叉核验）：EWASAtlas_get_probe（probe_id）、EWASAtlas_search_by_gene（gene_symbol；探针按关联数量排序，并提供 metadata.top_traits）、EWASAtlas_search_by_region（GRCh37/hg19 上的 chromosome、start、end，窗口最大为 1 Mb）、EWASAtlas_get_study（study_id，例如 "ES00033"：队列、组织、平台、每个探针的 p 值）、EWASAtlas_get_publication（pmid）。

hits = tu.tools.EWASAtlas_search_by_gene(gene_symbol="AHRR", limit=5)
top = hits["metadata"]["top_traits"][:3]          # e.g. smoking, maternal smoking
probe = hits["data"][0]                            # cg05575921, 131 associations
study = tu.tools.EWASAtlas_get_study(study_id=probe["associations"][0]["study_id"])

其他：ensembl_lookup_gene（需要 species='homo_sapiens'）、ensembl_get_regulatory_features（不能带 "chr" 前缀）、SCREEN_get_regulatory_elements、ChIPAtlas_*（需要 operation 参数）、SRA_search_experiments（library_strategy："ChIP-Seq"/"Bisulfite-Seq"/"ATAC-seq"）

第 7 阶段：全基因组统计

全局 beta 均值/中位数、探针方差、染色体密度、DMP 数量。

完整实现请参阅 CODE_REFERENCE.md。

---

常见模式

模式	关键步骤
差异甲基化	筛选探针 → 分组 → t 检验 → FDR → 阈值
年龄相关 CpG 密度	与年龄进行相关性分析 → FDR → 映射到 chr → 密度比
多组学缺失数据	提取 ID → 求交集 → 检查 NaN → 完整病例计数
ChIP-seq 注释	加载峰 → 注释基因 → 区域分类
甲基化-表达	对齐样本 → 计算相关性 → FDR → 负相关

---

GTEx 组织 ID

Whole_Blood、Liver、Lung、Breast_Mammary_Tissue、Brain_Cortex、Heart_Left_Ventricle、Kidney_Cortex、Thyroid、Adipose_Subcutaneous、Muscle_Skeletal

---

证据评级

等级	标准
Strong	padj < 0.01 且 abs(delta-beta) >= 0.2，并经重复验证
Moderate	padj < 0.05 且 abs(delta-beta) >= 0.1
Weak	padj < 0.05，但 delta-beta < 0.1
Insufficient	padj >= 0.05 或未进行重复验证

Delta-beta >= 0.2 表示效应显著。ChIP-seq：q < 0.01、FE >= 2 才有可信度。ATAC-seq 中 NFR < 150bp 表示调控活跃。始终应用 BH FDR。确认基因组版本一致。

---

局限性

• 不支持 pybedtools/pyBigWig：使用纯 Python 区间
• 以 Illumina 为主（450K/EPIC）；使用 t 检验/Wilcoxon（不使用 limma）
• 不进行峰值识别（假定已完成峰值识别）
• API 速率限制：每批约 20 个基因

参考文件

CODE_REFERENCE.md、TOOLS_REFERENCE.md、ANALYSIS_PROCEDURES.md、QUICK_START.md