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

omicverse を使用したバルク RNA-seq の DESeq2 解析

概要

ユーザーが t_deseq2.ipynb で紹介されている DESeq2 のワークフローを再現したい場合に、このスキルを使用します。raw 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. 遺伝子 ID をマッピングする
   • 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 がサイズファクターに関するエラーを返す場合は、対数変換済みデータではなく raw カウントを指定するようユーザーに案内します。
   • gene_id のマッピングは種に依存します。結果がまばらな場合は、正しいゲノムのペアをダウンロードするよう案内します。
   • 大規模なパスウェイライブラリを扱う場合、再帰制限を引き上げるか、プロットの前に上位 N 件の用語に絞り込む必要があります。

使用例

• 「処理群と対照群の反復サンプルに対して PyDESeq2 を実行し、エンリッチメント上位の WikiPathways 用語を強調表示してください。」
• 「log2(BaseMean) > 1 の遺伝子に DEG を絞り込み、フォールドチェンジのカットオフを自動選択して、ボルケーノプロットと箱ひげ図を作成してください。」
• 「GSEA 用の順位付き遺伝子リストを生成し、上位パスウェイのエンリッチメント曲線をプロットしてください。」

参考資料

• チュートリアルノートブック：t_deseq2.ipynb
• featureCounts のサンプルマトリックス：sample/counts.txt
• コピー＆ペースト用のクイックコマンド：reference.md