bulk-rna-seq-differential-expression-with-omicverse

使用 omicverse 进行 Bulk RNA-seq 差异表达分析

概述

按照此技能运行 t_deg.ipynb 中展示的端到端差异表达（DEG）工作流。此技能假设用户提供了原始基因级计数矩阵（例如，来自 featureCounts），并希望在 omicverse 中分析 Bulk RNA-seq 队列。

说明

1. 设置会话
   • 导入 omicverse as ov、scanpy as sc 和 matplotlib.pyplot as plt。
   • 调用 ov.plot_set()，使后续绘图采用 omicverse 的样式。
2. 准备 ID 映射资源
   • 需要将基因 ID 转换为基因符号时，指导用户通过 ov.utils.download_geneid_annotation_pair() 下载映射对，并将其存放在 genesets/ 下。
   • 说明可用的预构建基因组包括 T2T-CHM13、GRCh38、GRCh37、GRCm39、danRer7 和 danRer11；如有需要，用户也可以根据 GTF 文件生成自己的映射。
3. 载入原始计数
   • 使用 ov.pd.read_csv(..., sep='\t', header=1, index_col=0) 读取以制表符分隔的 featureCounts 输出。
   • 使用列表推导式从列名中移除末尾的 .bam 部分，以清理样本 ID。
4. 映射基因标识符
   • 运行 ov.bulk.Matrix_ID_mapping(counts_df, 'genesets/pair_<GENOME>.tsv')，将 gene_id 条目替换为基因符号。
5. 初始化 DEG 对象
   • 创建 dds = ov.bulk.pyDEG(mapped_counts)。
   • 使用 dds.drop_duplicates_index() 处理重复的基因符号，以保留表达量最高的版本。
6. 标准化并估算 size factor
   • 执行 dds.normalize()，计算 DESeq2 size factor，以校正文库大小和批次差异。
7. 进行差异检验
   • 将处理组和对照组的重复标签分别收集到列表中。
   • 调用 dds.deg_analysis(treatment_groups, control_groups, method='ttest') 执行默认的 Welch t 检验。
   • 提供可选方法：method='edgepy' 用于类似 edgeR 的检验，method='limma' 用于 limma 风格的建模。
8. 筛选结果并设置阈值
   • 注意，默认会保留低表达基因；需要时，可使用 dds.result.loc[dds.result['log2(BaseMean)'] > 1] 进行筛选。
   • 使用 dds.foldchange_set(fc_threshold=-1, pval_threshold=0.05, logp_max=6) 设置动态倍数变化和显著性阈值（fc_threshold=-1 会根据 log2FC 分布自动选择阈值）。
9. 可视化差异表达
   • 使用 dds.plot_volcano(title=..., figsize=..., plot_genes=... or plot_genes_num=...) 绘制火山图，突出显示关键基因。
   • 使用 dds.plot_boxplot(genes=[...], treatment_groups=..., control_groups=..., figsize=..., legend_bbox=...) 为每个基因生成箱线图；必要时调整 y 轴刻度标签。
10. 进行通路富集分析（可选）
    • 通过 ov.utils.download_pathway_database() 下载经整理的通路库。
    • 使用 ov.utils.geneset_prepare(<path>, organism='Mouse'|'Human'|...) 加载基因集。
    • 从 dds.result.loc[dds.result['sig'] != 'normal'].index 构建 DEG 基因列表。
    • 使用 ov.bulk.geneset_enrichment(gene_list=deg_genes, pathways_dict=..., pvalue_type='auto', organism=...) 进行富集分析。建议无法访问互联网的用户提供 background 基因列表。
    • 使用 ov.bulk.geneset_plot(...) 可视化单个基因库的结果，并使用 ov.bulk.geneset_plot_multi(enr_dict, colors_dict, num=...) 合并多个本体的结果。
11. 记录输出
    • 建议将 dds.result 和富集分析表导出为 CSV，以便后续报告使用。
    • 在笔记本以外的环境中运行时，建议使用 plt.savefig(...) 保存 matplotlib 生成的图表。
12. 防御性验证

    # Before DEG: verify treatment/control groups exist as column names
    all_cols = set(dds.result.columns) if hasattr(dds, 'result') else set(counts_df.columns)
    for g in treatment_groups + control_groups:
        assert g in all_cols, f"Sample '{g}' not found in count matrix columns"
    # Verify groups don't overlap
    assert not set(treatment_groups) & set(control_groups), "Treatment and control groups must not overlap"
    

13. 故障排除提示
    • 确保 treatment_groups/control_groups 中的样本标签与清理后的列名完全一致。
    • 确认已安装富集可视化所需的软件包（omicverse、pyComplexHeatmap、gseapy）。
    • 提醒用户，首次下载基因映射或通路数据库时需要访问互联网。

示例

• “我有一个小鼠肿瘤样本的 featureCounts 矩阵——使用 DESeq2 进行标准化，运行 t 检验差异表达基因分析，并在火山图中突出显示排名前 8 的基因。”
• “使用 omicverse 计算处理组和对照组重复样本之间的 edgeR 风格差异表达，然后对显著基因进行 GO 富集分析。”
• “指导我将 Ensembl ID 转换为基因符号、进行 limma 差异表达基因分析，并绘制 Krtap9-5 和 Lef1 的箱线图。”

参考资料

• 详细操作指南笔记本：t_deg.ipynb
• 用于测试的示例计数矩阵：sample/counts.txt
• 可快速复制粘贴的命令：reference.md