bulk-rnaseq

バルク RNA-seq

概要

このスキルは、生のシーケンスリードからエンリッチメント解析済みのパスウェイと図まで、根拠を示せるバルク RNA-seq 差次的発現解析を一貫して実行します。これは再実装ではなくルーターです。このリポジトリには、ほとんどの工程に対応する専用スキルがすでにあり、このスキルはそれらを適切な順序でつなぎ、唯一の実質的な不足部分（生リード → 遺伝子レベルのカウント行列）を補い、最終結果の信頼性を左右する実験デザインと QC の判断を徹底します。

ここでいう「根拠を示せる」とは、全工程に適用される次の 3 つを意味します。

• 再現可能 — パイプライン／ツールのバージョンを固定し、可能な場合はコンテナを使用し、パラメーターを記録し、乱数シードを固定します。
• 品質ゲート付き — QC を省略せず、定量の前・実行中・後に確認して対応します。
• 統計的に妥当 — 十分な反復数、生物学的背景に合った実験デザイン、適切なカウント処理、FDR を制御した検定を行います。

パイプラインは次のとおりです。FastQC／トリミング → アラインメント／定量（STAR／Salmon）→ カウント → DE（pydeseq2）→ エンリッチメント（pathway-enrichment）→ 図

このスキルを使う場面

ユーザーが次のことを望んでいる場合に、このスキルを使います。

• FASTQ ファイル（またはシーケンスラン）から、差次的発現遺伝子とパスウェイを取得する。
• nf-core/rnaseq を実行または設定する、あるいは STAR、Salmon、featureCounts でアラインメント／定量する。
• Salmon／STAR／featureCounts の出力を、DESeq2／PyDESeq2 用のカウント行列に変換する。
• 計算リソースを確保する前に、バルク RNA-seq 実験のデザイン（反復、バッチ、ストランド特異性）を検討または妥当性確認する。
• RNA-seq 解析全体の範囲を定め、連携するツールとスキルを決める。

これはバルク RNA-seq（サンプル = 生物学的検体）です。シングルセル／核データには scanpy を使い、DE 統計のみには pydeseq2 を、エンリッチメントのみには pathway-enrichment を使います。

パイプラインの概要

flowchart TD
    fastq["Raw FASTQ + samplesheet"] --> qc["FastQC + MultiQC"]
    qc --> trim["Trim: fastp / Trim Galore"]
    trim --> align["Align + quant: STAR and/or Salmon"]
    align --> counts["Gene-level counts matrix"]
    counts --> de["Differential expression"]
    de --> enrich["Pathway / GSEA enrichment"]
    de --> fig["Figures"]
    enrich --> fig
    nfcore["nf-core/rnaseq via nextflow skill"] -.->|"path A"| align
    manual["Standalone recipes (this skill)"] -.->|"path B"| align
    bridge["build_counts_matrix.py (this skill)"] -.-> counts
    pydeseq2skill["pydeseq2 skill"] -.-> de
    pwskill["pathway-enrichment skill"] -.-> enrich
    vizskill["scientific-visualization skill"] -.-> fig

2 つの上流パス — どちらかを選択

リード → カウントの工程は、2 つの方法で実行できます。どちらも同等の遺伝子カウントが得られます。状況に応じて選び、選んだ方法を最後まで使ってください。

パス A — nf-core/rnaseq を使う場面…	パス B — スタンドアロンツール を使う場面…
業界標準で、監査可能かつ引用可能なパイプラインを 1 つのコマンドで実行したい	サンプル数が少なく、各工程を学びながら確認したい
サンプル数が多い、または HPC／クラウドに拡張する予定がある	Nextflow／コンテナを利用できない、または環境に制約がある
再現性と包括的な MultiQC レポートを最も重視する	パイプラインでは利用できない標準外の工程が必要
→ nextflow スキルを通して実行	→ references/upstream-manual.md に従う

どちらを選ぶか迷った場合は、パス A を推奨します。nf-core/rnaseq は FastQC → トリミング → STAR／Salmon → 定量 → tximport → MultiQC を、適切にレビューされた妥当なデフォルト設定で連携して実行するため、最も根拠を示しやすい選択肢です。パス B は、透明性を重視する場合や環境に制約がある場合のために用意されています。

どちらのパスも遺伝子レベルのカウント行列に到達し、その後のワークフローは共通です。

セットアップ

# This skill's glue (bridge + handoffs) — Python
uv pip install pytximport pandas

# Downstream skills install their own deps:
#   pydeseq2 skill           -> uv pip install pydeseq2
#   pathway-enrichment skill -> uv pip install gseapy gprofiler-official

# Path A (nf-core): only Nextflow + a container engine are needed — see the `nextflow` skill.

# Path B (standalone tools): install via bioconda. Pin versions for reproducibility.
conda create -n rnaseq -c bioconda -c conda-forge \
  fastqc fastp trim-galore "star=2.7.11b" "salmon=1.10.3" subread multiqc

使用した正確なバージョン（パイプラインのリビジョン、ツールのバージョン、参照ゲノムとアノテーションのリリース）を記録してください。これらは方法のセクションに記載するものであり、解析の再現性を確保します。

クイックスタート

パス A — nf-core/rnaseq（推奨）

# 0. Validate the samplesheet first (catches the most common failures early)
python scripts/validate_samplesheet.py --samplesheet samplesheet.csv

# 1. Smoke-test the environment with tiny bundled data
nextflow run nf-core/rnaseq -r 3.26.0 -profile test,docker --outdir test_results

# 2. Real run: pin the revision, pick an aligner, pass a samplesheet + reference
nextflow run nf-core/rnaseq -r 3.26.0 \
  -profile docker \
  --input samplesheet.csv \
  --genome GRCh38 \
  --aligner star_salmon \
  --outdir results \
  -resume

nf-core/rnaseq は tximport を内部で実行するため、遺伝子カウントはすでに統合された状態で出力され、ブリッジスクリプトは不要です。DE には results/star_salmon/salmon.merged.gene_counts_length_scaled.tsv を使います。サンプルシートの形式、アライナーの選択、出力については references/upstream-nfcore.md を参照してください。実行エンジン／HPC／クラウド／コンテナの詳細については、nextflow スキルを使います。

パス B — スタンドアロンの STAR／Salmon（簡略版）

fastqc -o qc/ reads/*.fastq.gz                      # 1. QC raw reads
fastp -i s1_R1.fq.gz -I s1_R2.fq.gz \
      -o s1_R1.trim.fq.gz -O s1_R2.trim.fq.gz \
      --thread 4 -j s1.fastp.json                   # 2. Trim adapters/low-quality
salmon quant -i salmon_index -l A \
      -1 s1_R1.trim.fq.gz -2 s1_R2.trim.fq.gz \
      --gcBias --seqBias -p 8 -o quant/s1            # 3. Quantify (per sample)

完全な手順（FastQC、fastp／Trim Galore、STAR インデックス作成＋アラインメント＋--quantMode GeneCounts、Salmon のデコイ対応インデックス、featureCounts、ストランド特異性）については、references/upstream-manual.md を参照してください。

カウント → DE → エンリッチメント（どちらのパスも共通）

# Path B only: assemble a gene x sample counts matrix + metadata template for PyDESeq2
python scripts/build_counts_matrix.py --from salmon \
  --quant-dir quant/ --tx2gene tx2gene.tsv --output-dir counts/

# Then hand off (see the dedicated skills):
#   pydeseq2:           counts.csv + metadata.csv -> DE table (log2FC, padj, stat)
#   pathway-enrichment: rank by `stat` (GSEA) or padj+|LFC| hit list (ORA)
#   scientific-visualization / matplotlib: volcano, MA, heatmap, PCA, enrichment dotplot

ステージごとのワークフロー

上から順に進めます。各ステージでは、詳細を担当するスキルまたはファイルを示しています。デザインやQCのステージは省略しないでください。バルクRNA-seqの研究で最も問題が起きやすいのは、これらのステージです。

1. デザインとサンプルシート。 各グループに生物学的反復が3つ以上あることを確認し、バッチや交絡因子を特定して、比較対象を選びます。サンプルシートを作成し、scripts/validate_samplesheet.py で検証します。根拠とルール：references/design-and-qc.md。
2. リードの生QC。 ファイルごとに FastQC を実行し、MultiQC で集約します。塩基ごとの品質、アダプター含有量、重複、過剰出現を確認します。しきい値：references/design-and-qc.md。
3. トリミング。 アダプターと品質の低い末端を除去します（fastp または Trim Galore を使用）。FastQC を再実行して確認します。手順：references/upstream-manual.md（Path A ではこの工程を自動で行います）。
4. アラインメント／定量。 STAR（ゲノムアラインメント＋--quantMode GeneCounts）および／または Salmon（転写産物の準マッピング、デコイ考慮）を使用します。ストランド性を判定します。誤りやすく、気づかないうちにカウントが半減することがあるためです。詳細：references/upstream-manual.md、パイプラインのパラメーター：references/upstream-nfcore.md。
5. カウント行列の作成。 定量結果を、遺伝子 × サンプルの整数行列とメタデータのテンプレートに変換します（scripts/build_counts_matrix.py）。推定カウントと遺伝子IDのマッピングに関する注意点は、references/counts-and-handoff.md に記載されています。
6. 差次的発現 → pydeseq2 スキル。 counts.csv と metadata.csv を読み込み、デザイン（例：~batch + condition）を設定してフィッティングし、FDRを制御して検定します。QCとしてPCAとp値のヒストグラムを確認します。
7. エンリッチメント → pathway-enrichment スキル。 GSEA では、DESeq2 の stat で全遺伝子リストを順位付けします。ORA では、しきい値で絞り込んだヒットリスト（padj < 0.05、必要に応じて |log2FC| > 1）を渡します。まず遺伝子IDをシンボルに変換します。
8. 図 → scientific-visualization スキル。 ボルケーノプロット、MAプロット、サンプル間距離のヒートマップ、PCA、エンリッチメントのドットプロットに加え、QCの説明には MultiQC レポートを使用します。

counts → DE の橋渡し（重要なつなぎ）

上流／下流に担当スキルがない唯一のステージであるため、このスキルが担当します。scripts/build_counts_matrix.py は定量結果を pydeseq2 が期待する形式に変換します。

• Salmon（--from salmon）：サンプルごとの quant.sf を、counts_from_abundance="length_scaled_tpm" を使って遺伝子レベルに集約します（遺伝子レベルのDEにはこちらが適しています）。tx2gene マップが必要です。
• STAR（--from star）：各 ReadsPerGene.out.tab を読み込み、--strandedness（unstranded／forward／reverse）に対応する列を選択します。
• featureCounts（--from featurecounts）：統合済みの featureCounts 行列を解析します。

counts.csv（遺伝子 × サンプル、整数）と metadata_template.csv（サンプルごとに1行）を出力します。これらに必要な情報を入力してください。Salmon／RSEM のカウントは推定値（整数ではありません）ですが、PyDESeq2 が整数カウントを要求するため、整数に丸められます。これが length_scaled_tpm で許容される理由と、オフセットベースの DESeq2+tximport の手法との違いについては、references/counts-and-handoff.md を参照してください。このリファレンスでは、Ensembl からシンボルへのマッピング（エンリッチメントの前に必要）と、PyDESeq2 が要求する正確な向きについても説明しています。

よくある落とし穴

次の点が、誤った結果や再現性のないバルク RNA-seq 結果の主な原因です。

1. 反復数が少なすぎる。 各グループの生物学的反復が 3 未満では、検出力がほとんどなく、分散推定も不安定になります。シーケンスを深くするより、反復数を増やすほうが効果的です。
2. バッチと条件が交絡している。 処理群のサンプルすべてが、対照群とは異なる日に、または異なるレーンで処理されている場合、影響を切り分けることはできません。ランダム化し、既知のバッチをモデル化してください（~batch + condition）。references/design-and-qc.md を参照してください。
3. ストランド特異性の指定が誤っている。 誤った STAR の列や、featureCounts の -s／Salmon のライブラリータイプを選ぶと、リードの約半分が気づかないうちに破棄されます。Salmon の -l A を使うか、ストランド特異性を推定し、割り当てられたリードの割合を確認してください。
4. TPM／FPKM を DESeq2 に渡す。 DESeq2 が必要とするのは生の（または長さでスケーリングした）カウントであり、TPM／FPKM／正規化済みの値ではありません。このブリッジが変換を処理します。
5. 非整数カウント。 PyDESeq2 には整数が必要です。Salmon の推定値は丸めてください（ブリッジがこれを処理します）。
6. エンリッチメントに渡す遺伝子 ID が一致していない。 DESeq2 の出力は Ensembl ID であることが多い一方、Enrichr／MSigDB ではシンボルが必要です。pathway-enrichment の実行前に ID をマッピングしないと、「有意な結果なし」となります。
7. 定量後の QC を省略する。 DE の結果を信頼する前に、必ず PCA とサンプル間距離のヒートマップを確認してください。ラベルの入れ替わり、外れ値、隠れたバッチを見つけられます。
8. サンプル間でアライナーを混在させる。 すべてのサンプルを同じツール、バージョン、リファレンス、パラメーターで定量してください。
9. バージョンを固定しない。 「latest」のパイプラインやゲノムを使うと、結果を再現できません。-r、ツールのバージョン、ゲノム／アノテーションのリリースを固定してください。

他のスキルとの連携

• 上流処理の実行： nextflow（nf-core/rnaseq を実行。Path A、HPC／クラウド／コンテナー）。
• リファレンスデータ／遺伝子 ID： gget（ゲノム＋GTF には gget ref、ID マッピングには gget info／gget search）、database-lookup（Ensembl／NCBI）、biopython／pysam（FASTA／BAM の処理）。
• 差次的発現： pydeseq2（このスキルがカウントを渡す DE エンジン）。
• エンリッチメント： pathway-enrichment（ORA ＋ GSEA。付属の scripts/run_enrichment.py は DESeq2 の結果 CSV を直接読み込みます）。
• 図表とレポート： scientific-visualization、matplotlib、seaborn。方法と結果の記述には scientific-writing。
• 関連するが別の用途： scanpy（シングルセル）、statistical-analysis（多重検定の詳細）。

リファレンスファイル

詳しい情報が必要な場合は、該当するファイルを参照してください。各ファイルは単独で完結しています。

• references/upstream-nfcore.md — Path A：サンプルシートの形式、--aligner／--pseudo_aligner の選択、主要なパラメーター、salmon.merged.gene_counts*.tsv の出力、MultiQC、および pydeseq2 に渡すデータ。
• references/upstream-manual.md — Path B：FastQC、fastp／Trim Galore、STAR のゲノムインデックス＋アライメント＋--quantMode GeneCounts、Salmon のデコイ対応インデックス＋quant、featureCounts、およびストランド特異性の判定方法。
• references/counts-and-handoff.md — 定量出力を PyDESeq2 用の counts.csv／metadata.csv に変換する方法（pytximport、STAR の列選択、featureCounts）、整数と推定カウントの違い、Ensembl からシンボルへのマッピング、DE からエンリッチメントへのランク／ヒットリストの作成手順。
• references/design-and-qc.md — 実験デザイン（反復、バッチ、交絡、デザイン式）と QC 指標の解釈（マッピング率、重複率、rRNA、複雑度、PCA／外れ値）。妥当性を説明できるパイプラインの基盤となる内容です。

リソース

• nf-core/rnaseq: https://nf-co.re/rnaseq · STAR: https://github.com/alexdobin/STAR · Salmon: https://salmon.readthedocs.io
• fastp: https://github.com/OpenGene/fastp · Trim Galore: https://github.com/FelixKrueger/TrimGalore · MultiQC: https://multiqc.info
• pytximport: https://pytximport.complextissue.com · featureCounts (Subread): https://subread.sourceforge.net
• 手法の背景：Love et al. 2014（DESeq2）DOI 10.1186/s13059-014-0550-8 · Soneson et al. 2015（tximport）DOI 10.12688/f1000research.7563.2

関連スキル

• pydeseq2：DESeq2のPython移植版であるPyDESeq2を使用して、バルクRNA-seqカウントデータの差次的発現解析を実行します。
• pathway-enrichment：遺伝子リストまたはランキング済み遺伝子データを対象に、パスウェイおよび遺伝子セットのエンリッチメント解析を実行し、結果を解釈します。
• nextflow：Nextflowデータパイプラインとnf-coreワークフローをエンドツーエンドで構築、実行、デバッグします。