mcpField
B
bulk-rna-seq-deseq2-analysis-with-omicverse

bulk-rna-seq-deseq2-analysis-with-omicverse

github.com/omicverse数据与数据库1.4k57 次安装2 天前更新
代码仓库

关于 bulk-rna-seq-deseq2-analysis-with-omicverse

面向 bulk RNA-seq 分析者的 omicverse 技能，复现基于 PyDESeq2 的差异表达分析，并提供基因标识符映射、可视化和可选的 GSEA 富集流程。

主要功能

  1. 01读取 featureCounts 矩阵
  2. 02基因 ID 映射为符号
  3. 03PyDESeq2 差异表达分析
  4. 04差异基因筛选与阈值设置
  5. 05火山图和基因箱线图
  6. 06GSEA 富集与绘图

适用场景

  1. 01比较处理组与对照组的表达差异
  2. 02筛选差异基因并绘制火山图
  3. 03查看指定基因在组间的表达箱线图
  4. 04对排名基因运行 GSEA 并绘制通路曲线

安装

复制下面的命令，在终端中运行。

$npx skills add https://github.com/omicverse/omicverse --skill bulk-rna-seq-deseq2-analysis-with-omicverse

安装后在 Claude Code 中输入/bulk-rna-seq-deseq2-analysis-with-omicverse确认已加载。

附带文件reference.md

详细介绍

使用 omicverse 进行 Bulk RNA-seq DESeq2 分析

概述

当用户希望复现 t_deseq2.ipynb 中展示的 DESeq2 工作流时，使用此技能。内容包括加载原始 featureCounts 矩阵、将 Ensembl ID 映射为基因符号、通过 ov.bulk.pyDEG 运行 PyDESeq2，以及探索下游富集图。

说明

  1. 导入并整理表达矩阵
    • 调用 import omicverse as ov 和 ov.style()，统一可视化样式。
    • 使用 ov.io.read(..., index_col=0, header=1) 读取 featureCounts 的制表符分隔计数数据。
    • 使用 [c.split('/')[-1].replace('.bam', '') for c in data.columns] 去除列名末尾的 .bam。
  2. 映射基因标识符
    • 运行 ov.utils.download_geneid_annotation_pair()，确保存在适当的映射对。
    • 使用 ov.bulk.Matrix_ID_mapping(data, 'genesets/pair_<GENOME>.tsv') 将 gene_id 替换为基因符号。
  3. 初始化 DEG 对象
    • 基于映射后的计数创建 dds = ov.bulk.pyDEG(data)。
    • 使用 dds.drop_duplicates_index() 处理重复的基因名称，并在日志中确认成功。
  4. 定义对比组并运行 DESeq2
    • 将样本标签整理到 treatment_groups 和 control_groups 列表中，确保与列名完全一致。
    • 执行 dds.deg_analysis(treatment_groups, control_groups, method='DEseq2') 以调用 PyDESeq2。
  5. 筛选并调整阈值
    • 检查结果形状（dds.result.shape），并可按需筛选低表达基因，例如 dds.result.loc[dds.result['log2(BaseMean)'] > 1]。
    • 使用 dds.foldchange_set(fc_threshold=-1, pval_threshold=0.05, logp_max=6) 设置阈值，自动选取倍数变化截断值。
  6. 可视化差异基因
    • 使用 dds.plot_volcano(...) 绘制火山图，并汇总关键基因。
    • 按基因生成箱线图：dds.plot_boxplot(genes=[...], treatment_groups=..., control_groups=..., figsize=(2, 3))。
  7. 运行富集分析（可选）
    • 使用 ov.utils.download_pathway_database() 下载富集分析库，并通过 ov.utils.geneset_prepare 加载。
    • 使用 rnk = dds.ranking2gsea() 为 GSEA 排序基因。
    • 实例化 gsea_obj = ov.bulk.pyGSEA(rnk, pathway_dict)，并调用 gsea_obj.enrichment() 计算术语。
    • 使用 gsea_obj.plot_enrichment(...) 绘制富集气泡图，并使用 gsea_obj.plot_gsea(term_num=..., ...) 绘制 GSEA 曲线。
  8. 防御性验证 
    # Before PyDESeq2: verify count matrix contains raw integers (not log-transformed)
import numpy as np
if hasattr(data, 'values'):
    sample = data.values.flatten()[:1000]
else:
    sample = np.array(data).flatten()[:1000]
if np.any(sample != sample.astype(int)):
    print("WARNING: Data may not be raw counts. PyDESeq2 requires integer counts, not log-transformed.")
# Verify treatment/control groups match column names
for g in treatment_groups + control_groups:
    assert g in data.columns, f"Sample '{g}' not in count matrix columns: {list(data.columns)}"
  9. 故障排除
    • 如果 PyDESeq2 因 size factor 报错，请提醒用户提供原始计数（而非经过 log 转换的数据）。
    • gene_id 映射取决于物种；如果结果稀疏，请指导用户下载正确的基因组映射对。
    • 大型通路库可能需要提高递归限制，或在绘图前筛选前 N 个术语。

示例

  • “对处理组与对照组的重复样本运行 PyDESeq2，并突出显示富集程度最高的 WikiPathways 术语。”
  • “将 DEG 筛选为 log2(BaseMean) > 1 的基因，自动选取倍数变化截断值，并生成火山图和箱线图。”
  • “生成用于 GSEA 的排序基因列表，并绘制排名靠前通路的富集曲线。”

参考

评论

支持 Markdown