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

omicverse によるバルク RNA-seq の差次的発現解析

概要

このスキルに従って、t_deg.ipynb で紹介されている差次的発現解析（DEG）のエンドツーエンドのワークフローを実行します。ユーザーが生の遺伝子レベルのカウント行列（例：featureCounts の出力）を提供し、omicverse 内でバルク 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. 正規化とサイズファクターの推定を行う
   • dds.normalize()を実行してDESeq2のサイズファクターを計算し、ライブラリサイズとバッチの差を補正します。
7. 差次的検定を実行する
   • 処理群と対照群の反復サンプルラベルをリストにまとめます。
   • デフォルトのWelchのt検定には、dds.deg_analysis(treatment_groups, control_groups, method='ttest')を呼び出します。
   • 代替方法として、edgeR形式の検定にはmethod='edgepy'、limma形式のモデリングにはmethod='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にエクスポートすることを提案します。
    • ノートブック以外で実行する場合は、matplotlibで生成した図をplt.savefig(...)で保存するよう勧めます。
12. 防御的な検証

    # DEGの前に、処理群／対照群が列名として存在することを確認する
    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"
    # 群が重複していないことを確認する
    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 検定による DEG を実行して、ボルケーノプロットで上位 8 個の遺伝子を強調表示してください。」
• 「omicverse を使用して、処理群と対照群の反復間で edgeR 形式の差次的発現を計算し、有意な遺伝子に対して GO エンリッチメントを実行してください。」
• 「Ensembl ID をシンボルに変換し、limma による DEG を実行して、Krtap9-5 と Lef1 の箱ひげ図を作成する手順を教えてください。」

参考資料

• 詳細な手順を解説したノートブック：t_deg.ipynb
• テスト用サンプルカウントマトリックス：sample/counts.txt
• コピー＆ペースト用のクイックコマンド：reference.md