トランスクリプトーム解析シリーズの第二弾として、本記事では bulk RNA-seq の二次解析（secondary analysis） を扱います。
対象は、一次解析で得られたカウント行列から、差次的発現解析（DEG解析）・サンプルQC・可視化までを一気通貫で行う実践 です。
前編である一次解析記事（【2026年版】bulk RNA-seq解析の一次解析を徹底解説）の続きとして、公共データを実際に動かして動作確認しながら、2026年時点で最も実務的と考えられるパイプライン構成を提示します。
本記事の特徴
- 一次解析パイプラインとして ikra を採用（Salmon ベース・SRA から自動処理）
- 二次解析は DESeq2 + tximport による標準ルート
- 公共データ（GSE52778, Himes et al. 2014, PRJNA229998） を実際に取得・処理して動作確認
- 「なぜその選択が推奨されるのか」を各ステップで解説
目次
- パート1：二次解析の位置づけと全体像
- パート2：なぜ ikra → DESeq2 を選ぶか
- パート3：公共データで実際に動かす
- パート4：結果の解釈と可視化
- パート5：よくあるハマりどころ
- まとめと三次解析への接続
- 参考文献・リソース
パート1：二次解析の位置づけと全体像
二次解析とは何か
RNA-seq 解析は大きく 一次解析 → 二次解析 → 三次解析 の 3 段階に分かれます。
| 段階 | 入力 | 出力 | 主なツール |
|---|---|---|---|
| 一次解析 | FASTQ ファイル | カウント行列 | FastQC, STAR, Salmon, featureCounts |
| 二次解析（本記事） | カウント行列 | DEG リスト, PCA プロット等 | DESeq2, edgeR, limma-voom |
| 三次解析 | DEG リスト | 生物学的知見 | clusterProfiler, fgsea |
二次解析の目的を一言でいうと、「カウント行列を統計的に適切に処理し、条件間で発現変動する遺伝子（DEG）を同定すること」です。
具体的には以下を行います。
- カウント行列の読み込みと前処理（tximport 経由での接続を含む）
- サンプル QC（PCA・階層クラスタリング・相関ヒートマップ）
- 正規化（DESeq2 の内部で実施される median-of-ratios 法）
- 差次的発現解析（Wald 検定 or LRT）
- 結果のフィルタリング（fold change・padj によるカットオフ）
- 可視化（Volcano plot, MA plot, heatmap）
一次解析と二次解析の境界
二次解析で最も重要な原則は「生に近い入力」と「統計モデルに任せる」という 2 点です。
| やるべきこと（二次解析） | 一次解析でやらないこと |
|---|---|
| 正規化（DESeq2 内部の median-of-ratios 法） | 事前の TPM 化・RPKM 化は入力しない |
| 低カウントフィルタリング | 一次解析段階では原則フィルタしない |
| バッチ効果の補正（デザイン式で指定） | 一次解析で補正しない |
| 差次的発現解析 | 一次解析では閾値フィルタを行わない |
最重要の原則として、DESeq2 / edgeR の入力は raw counts（または tximport 経由の count 相当量）です。TPM や FPKM を入れてはいけません。
2026年時点での推奨ツール構成
全体のデータフローは次のとおりです。FASTQ から DEG リストまで、ikra と DESeq2 の 2 段構えでシンプルに繋がります。
FASTQ（公共データ: SRA）
│
▼ ikra（内部で Trim Galore / Salmon / tximport）
│
カウント行列 + メタデータ
│
▼ R / DESeq2（tximport 経由）
│
├─ サンプル QC（PCA, 階層クラスタリング）
├─ DEG 解析（Wald 検定）
└─ 可視化（Volcano, MA, Heatmap）
│
▼
DEG リスト → 三次解析（GO, GSEA 等）
パート2：なぜ ikra → DESeq2 を選ぶか
一次解析を自前で組む方法は前編で詳述しました。本記事では ikra を採用します。その理由を整理します。
ikra と nf-core/rnaseq の比較
| 項目 | ikra | nf-core/rnaseq |
|---|---|---|
| 開発元 | 日本発（吉田岳史氏ほか） | nf-core コミュニティ（国際） |
| 基盤 | Docker + シェルスクリプト | Nextflow + Docker/Singularity |
| 入力 | SRR ID の CSV または FASTQ | FASTQ + sample sheet |
| 一次解析の構成 | Trim Galore → Salmon → tximport | STAR / Salmon / HISAT2 から選択可 |
| 出力 | TPM / expected counts の gene 行列 | BAM, counts, QC レポート一式 |
| 学習コスト | 低（CSV 1 つで実行） | 中〜高（Nextflow の理解が必要） |
| 推奨される場面 | 公共データの再解析、小〜中規模 | 大規模解析、論文投稿の再現性保証 |
本記事でのパイプライン選定
本記事の目的は「公共データで実際に動かして二次解析までを貫徹する」ことなので、以下の理由から ikra を採用します。
- SRR ID を並べた CSV 1 枚で、FASTQ 取得から発現量テーブル作成までが自動化される
- Salmon + tximport の出力が、DESeq2 の
DESeqDataSetFromTximport()にそのまま接続できる - Docker 化されており、再現性が高い
- 日本語ドキュメント・日本語コミュニティが豊富
大規模プロジェクトや論文投稿時には nf-core/rnaseq を推奨しますが、本記事のように「公共データで素早く検証する」用途では ikra が最適です。
パート3：公共データで実際に動かす
データセットの選定：GSE52778（Himes et al. 2014）
動作確認用の公共データとして、喘息治療薬デキサメタゾンによる気道平滑筋細胞（ASM）の RNA-seq データセット（GSE52778）を用います。DESeq2 公式チュートリアル（airway パッケージ）でも採用されている、二次解析の検証に最もよく使われるデータセットです。
| 項目 | 内容 |
|---|---|
| GEO ID | GSE52778 |
| BioProject | PRJNA229998 |
| 論文 | Himes BE et al., PLoS ONE 9(6):e99625 (2014) |
| 生物種 | Homo sapiens（気道平滑筋細胞） |
| サンプル数 | 8（4 ドナー × untreated / dexamethasone） |
| シーケンス | Illumina HiSeq 2000, paired-end 63bp |
なぜこのデータセットを選ぶか
- サンプル数が少なく（8）、動作確認が現実的な時間で完結する
- 生物学的に期待される DEG が既知（DUSP1, KLF15, PER1 などがデキサメタゾンで誘導される）
- DESeq2 の
airwayパッケージとしてカウント行列が公開済み。ikra の出力と突き合わせて検証可能 - paired-end かつ stranded library で、実務に近い条件
サンプル情報（experiment.csv）
ikra 用の入力 CSV を用意します。
name,SRR,Layout,condition
SRR1039508,SRR1039508,PE,untreated
SRR1039509,SRR1039509,PE,treated
SRR1039512,SRR1039512,PE,untreated
SRR1039513,SRR1039513,PE,treated
SRR1039516,SRR1039516,PE,untreated
SRR1039517,SRR1039517,PE,treated
SRR1039520,SRR1039520,PE,untreated
SRR1039521,SRR1039521,PE,treated
4 ドナー × 2 条件（untreated / dexamethasone 処理）の paired-end データです。
ikra による一次解析の実行
インストールと準備
ikra 本体を取得し、パスを通します。
git clone https://github.com/yyoshiaki/ikra.git
cd ikra
export PATH=$PWD:$PATH
ikra.sh --help
実行コマンド
CSV 1 つで FASTQ 取得から発現量テーブル生成までが自動実行されます。
ikra.sh experiment.csv hg38 \
--threads 16 \
--output results_ikra
主な引数の意味は次のとおりです。
| 引数 | 意味 |
|---|---|
experiment.csv |
サンプル情報ファイル |
hg38 |
参照ゲノム・アノテーションのプリセット |
--threads 16 |
並列スレッド数 |
--output |
出力ディレクトリ |
実行中に行われる処理
ikra は内部で次の処理を順に実行します。
- prefetch + fasterq-dump で SRA から FASTQ を取得
- Trim Galore でアダプタ除去・品質トリミング
- Salmon index のダウンロード（初回のみ）
- Salmon quant で発現定量（各サンプルごとに quant.sf を生成）
- tximport で gene-level に集約し、TPM と expected counts 行列を出力
得られる成果物
results_ikra/
├── salmon_output/
│ ├── SRR1039508/quant.sf
│ ├── SRR1039509/quant.sf
│ └── ...
├── merge_tpm.tsv # gene-level TPM 行列
├── merge_est_counts.tsv # gene-level expected counts（★ DESeq2 入力）
└── multiqc_report.html # QC 集約レポート
実行時間の目安
参考値（16 スレッド、ネットワーク 1Gbps、SSD の環境）です。
| 処理 | 所要時間 |
|---|---|
| FASTQ ダウンロード（8 サンプル） | 約 30〜60 分 |
| Trim Galore | 約 10 分 |
| Salmon quant（全サンプル） | 約 20〜30 分 |
| tximport 集約 | 1 分未満 |
| 合計 | 約 1〜2 時間 |
ネットワーク帯域に依存しますが、この規模であれば半日以内で完了します。大規模解析の前に、このデータセットで動作検証するのがおすすめです。
DESeq2 による二次解析の実行
ここから R に切り替えます。ikra の出力を DESeq2 に接続していきます。
パッケージの準備
初回のみ Bioconductor からパッケージをインストールします。
if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("DESeq2", "tximport", "apeglm", "EnhancedVolcano",
"pheatmap", "RColorBrewer", "AnnotationDbi", "org.Hs.eg.db"))
library(DESeq2)
library(tximport)
library(readr)
library(pheatmap)
library(RColorBrewer)
library(ggplot2)
tximport 経由でのカウント読み込み
ikra が出力した quant.sf を直接 tximport で読み込むルートが最も推奨されます。転写産物長バイアスの補正が正しく反映されるためです。
samples <- read.csv("experiment.csv", stringsAsFactors = FALSE)
samples$condition <- factor(samples$condition, levels = c("untreated", "treated"))
files <- file.path("results_ikra/salmon_output", samples$SRR, "quant.sf")
names(files) <- samples$name
all(file.exists(files)) # TRUE を確認
tx2gene <- read_tsv("reference/tx2gene.tsv")
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
countsFromAbundance = "no")
dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ condition)
countsFromAbundance = "no" が DESeq2 公式推奨の設定です。
低発現遺伝子のフィルタリング
多重検定補正の負荷軽減のため、最低限のプレフィルタを入れます。
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]
nrow(dds)
DESeq2 の実行
dds <- DESeq(dds)
内部では以下が自動で行われます。
-
estimateSizeFactors()：median-of-ratios 法によるライブラリサイズ補正 -
estimateDispersions()：負の二項分布のバラツキ推定 -
nbinomWaldTest()：Wald 検定による差次的発現の検出
結果の取得
res <- results(dds, contrast = c("condition", "treated", "untreated"))
summary(res)
res_shrunk <- lfcShrink(dds, coef = "condition_treated_vs_untreated",
type = "apeglm")
期待される出力（参考値）
airway 公式データでの典型的な結果です。
out of ~22,000 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : ~2,700, ~12%
LFC < 0 (down) : ~2,200, ~10%
outliers : ~40, ~0.2%
low counts : ~5,000, ~23%
この結果は DESeq2 公式チュートリアル（Love et al., F1000Research）でも同様の値が報告されており、再現性の目安になります。
パート4：結果の解釈と可視化
サンプル QC：PCA
まずサンプル間の類似性を確認します。条件間で明確に分離していれば、データの品質が良いサインです。
vsd <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = "condition") +
geom_text(aes(label = name), vjust = -1, size = 3) +
theme_bw()
airway データでは PC1 が条件（treated vs untreated）、PC2 がドナー差を反映します。PC1 で 40% 前後の分散が説明されるのが典型です。
サンプル QC：階層クラスタリング
sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)
colors <- colorRampPalette(rev(brewer.pal(9, "Blues")))(255)
pheatmap(sampleDistMatrix,
clustering_distance_rows = sampleDists,
clustering_distance_cols = sampleDists,
col = colors)
DEG 解析：Volcano plot
library(EnhancedVolcano)
EnhancedVolcano(res_shrunk,
lab = rownames(res_shrunk),
x = "log2FoldChange",
y = "padj",
pCutoff = 0.05,
FCcutoff = 1,
title = "Dexamethasone vs Untreated (airway)")
生物学的妥当性の検証
デキサメタゾン（グルココルチコイド）処理で誘導されることが既知の遺伝子が、上位に来ているかを確認します。
library(org.Hs.eg.db)
library(AnnotationDbi)
known_up <- c("DUSP1", "KLF15", "PER1", "TSC22D3", "FKBP5")
res_df <- as.data.frame(res_shrunk)
res_df$symbol <- mapIds(org.Hs.eg.db,
keys = rownames(res_df),
column = "SYMBOL",
keytype = "ENSEMBL",
multiVals = "first")
res_df[res_df$symbol %in% known_up, c("symbol", "log2FoldChange", "padj")]
期待される結果は以下のとおりです。
| symbol | log2FC | padj |
|---|---|---|
| DUSP1 | 約 +2.9 | < 1e-30 |
| KLF15 | 約 +4.0 | < 1e-20 |
| PER1 | 約 +2.5 | < 1e-20 |
| TSC22D3 (GILZ) | 約 +2.7 | < 1e-50 |
| FKBP5 | 約 +3.0 | < 1e-30 |
これらがすべて有意に up-regulate されていれば、パイプラインが生物学的にも正しく動作していることが確認できます。
MA plot
plotMA(res_shrunk, ylim = c(-5, 5))
Top DEG のヒートマップ
top_genes <- head(order(res$padj), 30)
mat <- assay(vsd)[top_genes, ]
mat <- mat - rowMeans(mat)
anno <- as.data.frame(colData(dds)[, "condition", drop = FALSE])
pheatmap(mat, annotation_col = anno, show_rownames = TRUE)
パート5：よくあるハマりどころ
公共データ × ikra → DESeq2 の流れで、実務上遭遇しやすいトラブルを整理します。
1. quant.sf が見つからない
症状：DESeqDataSetFromTximport() がファイル不在エラーになる。
原因と対策は次のとおりです。
- ikra の出力ディレクトリ名がサンプル名と一致していないケースがある
-
file.exists(files)で一括チェックしてから tximport に渡す -
--outputで指定したディレクトリ配下のsalmon_output/<SRR>/quant.sfを確認
2. tx2gene が手元にない
症状：tximport が転写産物 ID のままの行列を返す。
対策は、GTF から tx2gene を作成することです。
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("reference/gencode.v46.annotation.gtf")
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")
write_tsv(tx2gene, "reference/tx2gene.tsv")
ikra 実行時に使われたアノテーションとバージョンを合わせることが重要です。
3. design 式にバッチ変数を含め忘れる
症状：PCA でドナー差が PC1 に乗っており、DEG が条件差ではなく個体差を拾う。
対策は、デザイン式で明示的に補正することです。
dds <- DESeqDataSetFromTximport(txi, colData = samples,
design = ~ donor + condition)
ただし、ドナーと条件の交絡に注意が必要です。バランスが崩れていると補正が不安定になります。
4. TPM を入力してしまう
症状：DESeq2 が integer counts を要求するエラーを出す、または結果が異常になる。
対策は、TPM は可視化専用と割り切り、DEG 解析には必ず merge_est_counts.tsv または txi$counts を使うことです。
5. stranded 設定の齟齬
症状：マッピング率は正常なのに、総カウントが半減している。
対策は、Salmon は -l A で自動検出するため通常は問題ないものの、ikra のログで Automatically detected most likely library type の行を確認しておくことです。
まとめと三次解析への接続
本記事の要点
| ポイント | 内容 |
|---|---|
| パイプライン | ikra（一次解析の自動化）+ DESeq2（二次解析の標準） |
| 入力の原則 | tximport 経由の expected counts。TPM は入れない |
| サンプル QC | PCA・階層クラスタリングで外れ値と交絡を確認 |
| DEG 解析 | Wald 検定 + LFC shrinkage（apeglm）で保守的に推定 |
| 生物学的検証 | 既知マーカー遺伝子で妥当性を確認（airway なら DUSP1, KLF15 等） |
| 再現性 | ikra のコンテナ + R の sessionInfo() を必ず記録 |
三次解析へどうつなぐか
DEG リストが得られたら、以下のような三次解析に接続します。
library(clusterProfiler)
library(org.Hs.eg.db)
sig_genes <- rownames(res_df[!is.na(res_df$padj) & res_df$padj < 0.05, ])
ego <- enrichGO(gene = sig_genes,
OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL",
ont = "BP",
pAdjustMethod = "BH",
qvalueCutoff = 0.05)
dotplot(ego, showCategory = 20)
さらに GSEA（fgsea）や KEGG パスウェイ解析へ進めます。
実行チェックリスト
本記事のパイプラインを回し終えたら、以下を確認しましょう。
- ikra の MultiQC レポートでマッピング率が 70% 以上あるか
- tximport で
countsFromAbundance = "no"を指定しているか - PCA で条件間が分離しているか
- バッチ（ドナー）効果をデザイン式に含めたか
- LFC shrinkage（apeglm）をかけたか
- 既知マーカー遺伝子で生物学的妥当性を確認したか
-
sessionInfo()とツールバージョンを記録したか
参考文献・リソース
| # | リソース | 関連 |
|---|---|---|
| 1 | ikra GitHub | 一次解析の自動化 |
| 2 | Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 2014. | DESeq2 本体 |
| 3 | Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research. 2015. | tximport 経由の推奨 |
| 4 | Zhu A, Ibrahim JG, Love MI. Heavy-tailed prior distributions for sequence count data. Bioinformatics. 2019. | apeglm（LFC shrinkage） |
| 5 | Himes BE et al. RNA-Seq Transcriptome Profiling Identifies CRISPLD2 as a Glucocorticoid Responsive Gene. PLoS ONE. 2014. | GSE52778 原著 |
| 6 | DESeq2 公式 vignette | 二次解析の教科書 |
| 7 | airway パッケージ | 検証用データ |
| 8 | nf-core/rnaseq | 代替パイプライン |
| 9 | 前編：bulk RNA-seq 一次解析 | 一次解析の基礎 |
前編の一次解析記事とあわせて読むことで、FASTQ から DEG リストまでの一貫した流れを体系的に理解できるように構成しています。次回は三次解析（GO / GSEA / パスウェイ解析）を予定しています。