ねらい: Bioconductorの airway RNA-seq データで変動解析の基本〜応用を一気に実装し、遺伝子発現以外（プロテオーム、メタボローム、メチローム等）にも一般化できる考え方・枠組みを整理する。

0. セットアップ

knitr::opts_chunk$set(message = FALSE, warning = FALSE, fig.width = 8, fig.height = 5)
set.seed(42)

必要パッケージの確認と（必要なら）導入

## ---- パッケージ：不足分のみ導入 -------------------------------------------
# Bioconductor側
bioc_pkgs <- c(
  "SummarizedExperiment","DESeq2","airway","limma","edgeR",
  "sva","AnnotationDbi","org.Hs.eg.db","apeglm"
)
# CRAN側
cran_pkgs <- c(
  "ggplot2","pheatmap","patchwork","tibble","dplyr","stringr",
  "matrixStats","ggrepel","msigdbr"
)

need <- unique(c(bioc_pkgs, cran_pkgs))
missing <- need[!suppressWarnings(sapply(need, requireNamespace, quietly = TRUE))]

if (length(missing) > 0) {
  if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
  bioc_missing <- intersect(bioc_pkgs, missing)
  cran_missing <- intersect(cran_pkgs, missing)
  if (length(bioc_missing) > 0) BiocManager::install(bioc_missing, ask = FALSE, update = FALSE)
  if (length(cran_missing) > 0) install.packages(cran_missing)
}

## ---- 読み込み（Startupメッセージ抑制） ------------------------------------
suppressPackageStartupMessages({
  library(SummarizedExperiment); library(DESeq2); library(limma); library(edgeR)
  library(airway); library(sva); library(GSVA); library(AnnotationDbi); library(org.Hs.eg.db); library(apeglm)
  library(ggplot2); library(pheatmap); library(patchwork); library(tibble); library(dplyr); library(stringr)
  library(matrixStats); library(ggrepel); library(msigdbr)
})

## ---- ggplot2 のデフォルトテーマ -------------------------------------------
theme_set(ggplot2::theme_bw(base_size = 12))

パッケージ構成と環境管理の方針

Bioconductor（SummarizedExperiment, DESeq2, limma, edgeR, sva, AnnotationDbi, org.Hs.eg.db, apeglm 等）と、CRAN（ggplot2, pheatmap, patchwork, dplyr, matrixStats, ggrepel, msigdbr 等）を 用途別に説明 すると：

  • データ構造：SummarizedExperiment は 行=特徴量（遺伝子）×列=サンプル の行列と rowData/colData をひとまとめにした中核のコンテナ。

  • 変動解析：DESeq2 と edgeR/limma は カウント と 連続量 の両輪。負の二項（NB）/精密化 t 検定（moderated t）の代表実装。

  • 正規化・可視化：matrixStats はベクトル化統計、ggplot2/pheatmap は図化、patchwork は図のレイアウト統合に用いる。

  • アノテーション注釈：AnnotationDbi と org.Hs.eg.db は ID 変換・遺伝子注釈。msigdbr は経路・遺伝子集合の定義源となる。

  • 効果量縮約：apeglm は小標本での log2 fold change の過大評価を抑制 するベイズ的縮約器である。

Bioconductor はリリースごとに挙動が変わるため、R と Bioconductor のバージョンを必ず記録 するとよい。

1. コンセプト：オミクス横断の変動解析ワークフロー

  • 目標: 条件A vs 条件Bの平均の違い（効果量）を、ノイズやバッチの影響を抑えつつ、多数の特徴量（遺伝子/タンパク/代謝物/ピーク/部位）に対して検出する。

  • 抽象モデル：

    • カウント系: \(\text{counts}\times{gi} \sim \text{NB}(\mu\times{g_i}, \alpha_g)\)、\(\log \mu_{gi} = \log s_i + X_i^\top\beta_g\)（DESeq2/edgeR）

    • 連続量系: \[y_{gi} = X_i^\top\beta_g + \varepsilon_{gi}\]、\(\operatorname{Var}(\varepsilon_{gi})\)は平均依存→voomで重み化しlimmaの経験ベイズで分散を安定化。

  • 論点（どのオミクスでも共通）:

    1. QC/前処理（欠損・外れ値・低発現除去・正規化）

    2. 設計行列（交絡の明示、ブロック/繰返しの扱い）

    3. 尤度/分散モデルの選択（NB vs 正規）と分散の借用（moderation）

    4. 多重検定（FDR制御）と効果量の縮約（小標本の過大評価抑制）

    5. 生物学的解釈（経路/集合の単位に持ち上げる）

データタイプとモデリングの違い

遺伝子発現データ（RNA‑seq）は主にカウントデータであり、読取り数はポアソンや負の二項などの離散分布に従う。

カウントでは平均と分散が密接に関連し平均が大きいほど分散も大きくなる（平均–分散依存）性質をもつ。ライブラリサイズや構成の違いに応じた正規化と分散安定化が重要である。

たとえば RNA‑seq のクラスタ解析では TPM のように遺伝子長とシーケンス深度を補正した値を用い、変動解析では raw count に基づいた NB‑GLM を採用する。

ポアソン/負の二項モデルはカウントの自然な出発点で、DESeq2::vst() は NB モデルに基づく 遺伝子ごとの変換 により、平均に依存しがちな分散をならし、PCA や距離計算といったユークリッド幾何の前提を満たしやすくする。

一方、プロテオームやメタボロームでは測定されるのはイオン強度などの連続値であり、対数正規 に近い分布で、正規近似が働きやすいが、電離効率や保持時間に起因するばらつきを含む。

この違いが quantile正規化と分散安定化（VST） を必要とする理由となる。VSN は技術反復間のばらつきを抑えるのに有効であると報告されている。

下図は、シミュレートしたカウントデータ（負の二項／ポアソン）と正規分布に従う連続強度データのヒストグラムの比較例である。

離散カウントでは 0 が多く裾が重い一方、連続強度はおおむね対称で中心付近に集中していることがわかる。

測定タイプと代表的な解析手法

オミクス層 測定タイプ モデル 正規化の例
RNA-seq/ATAC-seq カウント（離散） NB-GLM (DESeq2 / edgeR) RLE/TMM/UpperQ + VST
プロテオーム (DDA/DIA) 連続強度（イオン） 線形モデル (voom+limma) 参照チャネル/内部標準, PQN, VSN
メタボローム 連続強度（ピーク面積） 線形モデル/ロジ回帰/加重回帰など PQN, ログ変換, VS/N
DNAメチル化/ChIP など カウント（バイサルファイト）または強度 NB-GLM や β-binomial モデルなど RLE/TMM、β-mixture、ログ変換

オミクス横断のワークフロー（概要）

以下のように、原データから変動解析解釈までの流れを段階的に整理すると全体像が掴みやすくなる。

  • QC：欠損・外れ値・低発現を検査し、不良サンプルや特徴を除外する。

  • 正規化：ライブラリサイズや電離効率などスケール差を調整する（RLE/TMM/UpperQ、参照チャネルなど）。

  • 設計行列の定義：群とブロック、交互作用を明示する。ペアデザインやブロック共変量が含まれる場合は設計行列で吸収する。

  • モデル適合：離散カウントには NB-GLM (DESeq2/edgeR)、連続強度には線形モデル (voom+limma) を用いる。

  • 効果量の縮約と多重検定：経験ベイズや事前分布を使い効果量を安定化し、FDR を制御する。

  • 解釈と集合解析：MA/Volcano/PCA/ヒートマップで変動解析を可視化し、GSVA/fgsea 等で経路レベルにまとめる。

この図では、原データ（counts や強度）から QC → 正規化 → 設計 → モデル → 効果量縮約＋FDR → 可視化 → 集合解析という流れを示している。これらの手順は遺伝子発現以外のオミクスにも共通する。

重要なのは、ライブラリ調製 → シーケンス → マッピング → カウント集計 という流れのどこでもバイアスは入り得る点である。

特に、

  1. ライブラリサイズ（総読取り数）

  2. 組成（高発現遺伝子の支配）

  3. アラインメントの曖昧さ（多重割当）

これらは下流の変動解析検定に直接影響する。

以降の QC では、可視化→最小限の補正 の順に、解析に不要なばらつきを落としていくことにする。

実例：生データと正規化後の比較

以下では、Airway データから最初の数サンプルの raw counts、RLE サイズ係数による正規化、VST 変換の違いを簡単に比較する。prepare_data チャンクで se と cts は読み込まれている。

##                 SRR1039508 SRR1039509 SRR1039512
## ENSG00000000003        679        448        873
## ENSG00000000005          0          0          0
## ENSG00000000419        467        515        621

このように、正規化と分散安定化によってデータのスケールと分布が整うことがわかる。実務では各 P の手順を繰り返し適用して、最終的な変動解析に備える。

2. QCと前処理 — 原理・チェックリスト

本節では 「何を、なぜ、どこで」 チェック・制御・処理するのかを、オミクス横断の原理としてまとめ、そのうえで例として シークエンシング（RNA-seq想定） と 質量分析（プロテオーム/メタボローム） の要点に落とし込む。

実際には”装置→測定→前処理→統計”の各層で 系統誤差（bias） と 偶然誤差（variance） を切り分け、シグナル保存 と ノイズ抑制 のトレードオフを常に意識する。

2.1 QCの設計思想（共通の5原則）

  1. 予防が最良のQC：実験計画・ランダム化・ブロッキング・QCサンプル/スパイクインの設計で、統計補正に頼りすぎない。

  2. 層で考える：\(サンプル層（採取/処理）→ 計測層（装置/ラン）→ データ層（行列の欠損・外れ）→ 推定層（モデル/分散）\) のどこに問題があるかを同定。

  3. 外部基準＋内部整合：外部基準（iRT/標準品/スパイクイン/遺伝子体カバレッジ等）と、内部整合（反復・相関・ペア一致）を併用。

  4. ルールを先に決める：フィルタ閾値・外れ値規則・バッチ補正の要否は、解析前にドキュメント化（“研究計画に準拠”）。

  5. リークを避ける：変動解析検定の目的変数に触れる操作（例：群ラベルを使った過学習的なノーマライズ/除外）は厳禁。

原則 狙い 代表的な定量指標
P1. 同定・整合性（ID/メタデータ） サンプル取り違え・属性不一致の早期検知 メタデータ整合率、性別推定一致、重複/近縁性（IBD/IBS）、バーコード衝突率
P2. 信号品質・S/N 技術ノイズ・低品質測定の除外 Qスコア分布・%Q30、マッピング率、重複率、欠損率、TSS/FRiP、TIC安定性、質量誤差ppm
P3. 系統バイアス・汚染 バッチ・ドリフト・外来汚染の検知/補正 バッチ寄与（R²/ANOVA）、ドリフト指標（QC列のRT/強度トレンド）、decontam率、交差汚染推定
P4. 正規化・較正 機器/ロット差のスケール調整 内部標準比、QCベースLOESS、ライブラリサイズ補正、VSN/量子化、Noob/FunNorm
P5. 再現性・妥当性/外れ値 QCリピートの再現・生物学的妥当性確認 QCサンプルCV（%CV）、相関（Pearson/Spearman）、既知ポジコン回収率、クラスタ安定度、外れ値指標（RLE/NUSE等）

以降の各表では項目末尾に [P#] で該当原則を付記。

2.1a シークエンシング系：コアスキル

総論（本質）

配列読取りは「分子の離散サンプリング（counts）」であり、読取り数の総量差（ライブラリサイズ・構成バイアス）、断片化と配列特性に由来する位置・塩基依存バイアス、割当ての不確実性（多重割当・救済割当）を必然的に生む。

したがって、盲目的（blind）に(i)身元整合、(ii)信号健全性、(iii)系統偏り、(iv)総量/分布整合、(v)再現性、の順で可視化→最小補正→共通尺度化し、解析用アーティファクトを凍結して引き渡すのが核となる。

5つのコアスキル（工程別）

RNA‑seq などの配列データでは、以下の 5 つの観点を順に確認し、目的と手段をセットで意識することで QC と前処理を体系化できる。

  • P1 同定・整合性：目的はサンプル ID と実験設計（群・ブロック）の整合を確認することです。手段として sample swap や index hopping の有無を検査し、性別や既知マーカーで照合して配列方向と参照・注釈の版を決めて記録する。

  • P2 FASTQ 品質：目的は FASTQ リードが統計的に扱える品質かどうかを見極めることである。手段として Q スコアの分布やアダプタ残存、PCR 重複、汚染などを監査し、必要に応じてトリミングや重複除去（UMI 誤り補正付き）を行う。

  • P3 MAP/COUNT：目的は参照への割当てと集計規約を固定し、割当ての不確実性や構成バイアスを抑えることである。手段としてユニーク／EM 割当、外れ長さや低品質断片の除外、集計単位（遺伝子 vs 転写産物、曖昧リード・重複・アンチセンス、長さ補正）と注釈・参照版を事前に規定する。

  • P4 スケーリングと分散安定化：目的はライブラリサイズ差と平均–分散依存といった離散カウント特有のアーティファクトを緩和することである。手段として RLE／TMM／UpperQ などのサイズ係数でスケールを合わせ、VST・rlog・voom などで分散を安定化します（TPM/RPKM は記述用であり比較には用いない）。

  • P5 再現性と外れ値：目的は技術的再現性の確保と外れサンプルの検出である。手段として QC 指標や RLE／MA プロット、Cook’s 距離などを評価し、合否基準を SOP として文書化した上で counts とメタデータ、QC レポートを凍結・引き渡す。

P 主たる工程 核となる問い 代表指標（例） 推奨処置（最小） Rパッケージ／関数（例）
P1 同定・整合性 FASTQ / デマルチ＋初期アライン設定 サンプルと設計情報は整合しているか（sample swap / index hopping / strandness） ラベル整合、性別・既知マーカー照合、配列方向 照合・監査を先行、方向/参照/注釈の版を固定して記録 SummarizedExperiment::colData／dplyr（設計確認）
RSeQC infer_experiment.py（外部；方向性）
SNPRelate::snpgdsIBD（遺伝的照合がある場合）
P2 信号品質 FASTQ → MAP 読取りは統計的に扱える品質か（Q30／アダプタ／重複／被覆／位置決め） %Q30、アダプタ残存、重複率、マップ率、遺伝子体被覆（5′/3′） 最小トリム、UMI誤り訂正＋重複除去、極端断片の除外、再調製判断 ShortRead::qa／Rqc（FASTQ品質）
fastp／cutadapt（外部；トリム）
Picard（外部；重複）
Rsamtools::idxstatsBam（配置統計）
RSeQC geneBody_coverage.py（外部；被覆）
Rsubread::featureCounts（rRNA率など）
P3 系統的偏り MAP / COUNT 設計要因や集計規約・救済割当（EM/rescue）による構成バイアスはないか PCAのバッチ分離、PC1–ライブラリ規模相関、救済割当の影響 ランダム化設計、サイズ係数の適用（RLE/TMM）→必要最小限のバッチ補正（ComBat/RUV） edgeR::cpm／edgeR::calcNormFactors("TMM")
DESeq2::estimateSizeFactors
stats::prcomp＋ggplot2（PCA）
sva::ComBat、RUVSeq::RUVg/RUVs
limma::removeBatchEffect
P4 正規化・較正 SCALE 離散カウント特有の総量差と平均–分散依存をどう整えるか RLE中央値≈0、MAのA依存の緩和、PCに規模が乗らない サイズ係数（RLE/TMM/UpperQ）→分散安定化（VST/rlog/voom）※blind／TPMは記述用 DESeq2::estimateSizeFactors（RLE）
`edgeR::calcNormFactors(“TMM” “upperquartile”)<br>DESeq2::vst／DESeq2::rlog<br>limma::voom（重み付け）<br>matrixStats::rowMedians（RLE計算）<br>DESeq2::plotMA`
P5 再現性・外れ REPRO（引き渡し前） 技術再現と外れの扱いをどう確証するか サンプル間相関・距離、QC%CV、Cook’s距離、既知陽性回収 合否SOPで判定・記録、外れの方針（除外/別扱い）を適用→アーティファクト凍結 DESeq2::plotPCA、pheatmap::pheatmap（距離）
DESeq2::plotMA（俯瞰）
assays(dds)[["cooks"]]（Cook’s）
stats::cor／ggplot2（相関可視化）

ワークフローの例：QC（P1–P3）→ 正規化/安定化（P4；blind）→ 再確認（P2–P4）→ 再現性確認（P5）→ 引き渡し（解析用アーティファクト凍結）。

2.1b 質量分析系：コアスキル

総論（本質）

質量分析は「連続強度のサンプリング（イオン強度）」であり、電離効率（matrix effect）・分離時間の揺らぎ（保持時間ドリフト, RT drift）・質量校正誤差（ppm error）・測定順の劣化（run-order drift）・同定の確率性（DDA）やライブラリ依存性（DIA）・チャネル比の縮小（ratio compression, TMT）・欠測の機構（MNAR/MAR）を必然的に伴う。

したがって、盲検的（blind）に

(i)身元整合

(ii)信号健全性

(iii)系統偏り

(iv)総量/分布整合

(v)再現性

の順で可視化→最小補正→共通尺度化し、解析用アーティファクトを凍結して引き渡すのが核となる。

5つのコアスキル（工程別）

質量分析データでも、以下の 5 つの観点に沿って QC と前処理を進める。各項目で「目的」と「手段」を明確に区別すると初心者でも流れを理解しやすくなる。

  • P1 同定・整合性：目的はラベル割付や内部標準、参照チャネルが設計どおりに反映されているか確かめることである。手段として割付表や参照比、内部標準の回収率を照合し、版情報をログに残す。

  • P2 信号健全性：目的は測定の健全性をモニタリングすることである。手段として総イオン量 (TIC)、保持時間 (RT)、質量精度 (ppm)、同定率、ピーク形状を時系列で監視し、QC 混合物をラン全体に挿入して必要に応じて再校正やソース／カラムの整備、パラメータ最適化を行う。

  • P3 系統的偏りと汚染：目的は測定順ドリフト、キャリーオーバー、バッチ差などの構成バイアスを抑制することである。手段として QC 基盤の時系列補正 (QC‑RLSC)、ブランク監視、ランダム化配置、残差バッチ補正 (ComBat 等) を組み合わせる。

  • P4 スケーリングと分布調整：目的は注入量や電離効率の差、比圧縮に起因するスケール不一致や分布歪みを解消することである。手段として内部標準や参照チャネル、PQN、分位合わせ／VSN、log/glog などを用いて正規化する（TMT は参照チャネル基準）。

  • P5 再現性と欠測管理：目的は QC 反復の再現性を確認し、欠測データを適切に扱うことである。手段として %CV や相関、既知スパイクの回収率を指標に合否を判定し、欠測機構 (MNAR/MAR) に応じて MinProb/QRILC や kNN などの補完法を選び、アーティファクトを凍結する。

質量分析系コアスキル（P1–P5対応＋Rパッケージ）

P 主たる工程 核となる問い 代表指標（例） 推奨処置（最小） Rパッケージ／関数（例）
P1 同定・整合性 RAW/設計＋ラベル/標準確認 サンプル設計は整合しているか（TMT割付・内部標準・参照/iRT） 参照比の逸脱、内部標準回収、iRT整合 割付表照合、参照/内標の回収監査、設計ログ固定 MSstatsTMT（設計/集約）／QFeatures（メタ連結）／PTXQC（外部；MaxQuantログ）
P2 信号品質 RAW/CHROM 測定は健全か（TIC/RT/ppm/ID率/ピーク形状） TIC安定、RTシフト、質量誤差、ID数、ピーク幅 QC列による時系列監視、再校正・ソース/カラム整備、パラメータ最適化 MSstatsQC2（時系列QC），Spectra/mzR（TIC, ppm），xcms（クロマト処理），IPO（パラメタ最適化）
P3 系統的偏り ALIGN/ID/QUANT ドリフト・キャリー・バッチ差はないか QCトレンド、blankピーク、PCAでバッチ分離 QC-RLSCで時系列補正、ブランク監視、残差バッチ補正（ComBat/proBatch） statTarget::shiftCor（QC-RLSC），proBatch，sva::ComBat，ggplot2＋stats::prcomp（PCA）
P4 正規化・較正 SCALE/NORM スケール不一致・分布差・比圧縮をどう整えるか 分位プロファイル、VSN適合、チャネル比の安定 内標/参照・PQN・quantile・VSN・log/glog（TMTは参照基準） NormalizeMets::pqnorm，limma::normalizeBetweenArrays（quantile），vsn::vsn2，MSstatsTMT::dataProcess（TMT比補正）
P5 再現性・外れ REPRO/引き渡し 技術再現と欠測の扱いは妥当か QC%CV、rep相関、欠測パターン（MNAR/MAR） 合否SOPで判定、左側代入（MinProb/QRILC）／kNNなど機構別補完、アーティファクト凍結 MSstats（QC評価/要約），DEP（可視化・前処理），imputeLCMD::impQRILC/MinDet，DMwR::knnImputation（代替）

ワークフロー例：ピーク検出/整列 → QC時系列補正（P3） → スケール/分布補正（P4；blind） → 再確認（P2–P4） → 再現性確認（P5） → 引き渡し（アーティファクト凍結）。 DDA と DIA：DDA は取得の確率性・欠測（MNAR）に配慮、DIA はライブラリ品質・整列（DIAlignR）の妥当性が鍵。

TMT/iTRAQ：参照チャネル基準で比を安定化し、比圧縮に注意（MSstatsTMT）。LFQ：分位/VSNと外れ頑健な集約（MSstats）。メタボロミクス：QC-RLSC→PQN/内標→log/glogが定石。

2.2 共通チェックリスト

  • サンプル識別：ラベル取り違え/重複。遺伝的指紋（SNP）や性染色体マーカー、メタデータ整合で検出。

  • 量の健全性：総シグナル（ライブラリサイズ/TIC）、非欠損率、上位特徴への偏り（トップ100比率）。

  • 技術一貫性：QCリプリケートのCV、バッチ/ラン間のドリフト（RT/質量精度/GC含量/3’バイアス）。

  • 構造的交絡：群とバッチ・施設・日付が完全交絡していないか（→推定不能）。

  • 外れ値：PCA/距離ヒートマップ/ロバスト指標（MAD）で検出、原因特定→処置（やむを得ない場合のみ除外）。

各論

RNA-seqの場合

段階 原理（なぜ） データ現象 検出・管理の要点
P1 同定 取り違え・混線 属性と不一致、近縁性の破綻 生物学的指標で照合、読み取り方向の推定
P2 信号 読み取り誤差・重複・端への偏り 低品質・不要断片・位置決め不良・被覆偏り 品質・被覆・重複・不要断片の監査と是正
P3 かたより 設計要因差・構成バイアス バッチ集団化、全体量が主因 量合わせ→残差的な設計要因の補正
P4 正規化 全体量差・平均–分散依存 スケール不一致、MAの歪み 基準比でのスケール合わせ＋分散安定化
P5 再現 劣化・希少事象 繰り返し不一致、外れ 相関・ばらつき・対照回収で合否判断

P1 同定・整合性

  • 原理：試料の取り違えや多重化バーコードの混線が起こりうる（用語：sample swap, index hopping, ストランド性）。

  • データ現象：生物学的属性とメタ情報が不一致／近縁な試料同士の距離が不自然（例：XIST, RPS4Y1 の不一致、IBS/IBDの異常）。

  • 検出・管理：性別・既知マーカー・遺伝的指標の照合、配列方向の推定を初期に行う（sex-check, SNPRelate, infer_experiment）。

P2 信号品質

  • 原理：読み取り誤差・PCR重複・不要断片・増幅効率の差・端への偏りが生じる（Qスコア, PCR duplicates, adapter contamination, GC bias, 5′/3′ bias）。

  • データ現象：%Q30低下、アダプタ残存、位置決め率低下、遺伝子体被覆の片寄り、重複率上昇、rRNA比高値（mapping rate, gene body coverage, duplication rate, rRNA mapping）。

  • 検出・管理：品質監査→必要なトリミング・重複処理・再調製（FastQC/fastp, Picard, RSeQC, UMI dedup）。

P3 系統的偏り・汚染

  • 原理：測定日・試薬・機器差などの設計要因と、相対量の歪み（batch effect, compositional bias）。

  • データ現象：PCAが群ではなくバッチで分離、第一主成分とライブラリ量が強相関（PC1–library size correlation）。

  • 検出・管理：実験計画での分散配置＋サイズ係数での構成補正→必要なら統計的補正（RLE, TMM, ComBat, RUVg/RUVs）。

P4 正規化・較正

  • 原理：総量差と平均–分散依存が推定を歪める（library size difference, mean–variance relationship）。

  • データ現象：スケール不一致、MA依存の歪み、RLE中央値の逸脱（MA plot, RLE）。

  • 検出・管理：基準比にもとづくスケーリング＋分散安定化（median-of-ratios (RLE), TMM, upper-quartile, VST/rlog, voom）。

P5 再現性・外れ値

  • 原理：劣化・操作差・希少事象が技術ばらつきを増やす（technical variability, degradation）。

  • データ現象：繰り返しの相関低下、外れ値点、既知陽性の未回収（Cook’s distance, distance heatmap）。

  • 検出・管理：相関・CV・外れ値診断で合否判定し、原因に応じて除外／再測定／ロバスト化（QC replicate, positive controls）。

どのように実装するか（RNA-seqの実装例）

この演習で身につけるRのコア操作（QC／前処理）

到達目標（3つ）

  1. QCの可視化と基礎統計：サンプル単位の要約（総量・非欠損・集中度）と、低次元表示（PCA）で品質と偏りを診断できる。

  2. 前処理の最短手順：低カウント除去 → サイズ係数 → 分散安定化を正しく実装・検証できる（blind）。

  3. 再現性の確認と引き渡し：PCA／距離／外れ指標で合否を判定し、解析用アーティファクト（行列・メタ・QC図）として凍結できる。

QC（P1・P2・P3・P5）— コアスキル早見

P 何を確かめる（目的） Rで何をする（操作） 主な関数・パッケージ（例）
P1 同定・整合性 サンプルIDと設計（群・ブロック）の整合 列名の一意性・設計表のクロス集計 duplicated()、SummarizedExperiment::colData、table()
P2 信号品質 総量・非欠損・上位集中（明白な不良の検出） qc_sample_summary()で要約列を作成し可視化 colSums()／colMeans()、自作qc_sample_summary()、ggplot2、ggrepel
P3 系統的偏り 規模やバッチが主なばらつき要因になっていないか edgeR::cpm(log=TRUE)→prcomp()→PC1とログ総量の相関 edgeR::cpm、stats::prcomp、cor()、ggplot2
P5 再現性・外れ 群内一貫性・外れの有無 VST後のPCA／サンプル距離ヒートマップ／Cook’s距離の概観 DESeq2::plotPCA、pheatmap::pheatmap、assays(dds)[["cooks"]]

参考（FASTQ段階のQC）：実データでは FastQC / fastp / RSeQC で Q30・アダプタ・ストランド・gene body coverage を併用（Airway例では省略）。

前処理（P4中心）— コアスキル早見

ステップ 目的 Rで何をする（操作） 主な関数・パッケージ（例）
低カウント除去 検出力の低い特徴を除き誤検出を抑制 デザインを踏まえてフィルタ edgeR::filterByExpr
サイズ係数 総量差（ライブラリサイズ・構成差）を補正 RLEのサイズ係数を推定 DESeq2::estimateSizeFactors（代替：edgeR::calcNormFactors）
分散安定化 平均–分散依存を緩和し可視化・距離計算を安定化 VST（blind）を適用 DESeq2::vst（代替：limma::voom）
妥当性確認 正規化の効きすぎ／足りなさを点検 RLE箱ひげ（中央値≈0）、PCA、距離図 matrixStats::rowMedians＋boxplot()、DESeq2::plotPCA、pheatmap
引き渡し 解析用アーティファクトの凍結 行列・メタ・QC図を保存 saveRDS()、readr::write_csv()

任意：バッチ補正は正規化後に必要最小限（sva::ComBat など）。DE（変動解析）の設計式はQC/前処理の後段で初めて投入する。

最小コマンド列（この順で実行できればOK）

# P2: QC要約
qc_df <- qc_sample_summary(cts)

# P3: 偏りの事前診断
logCPM <- edgeR::cpm(cts, log=TRUE, prior.count=1)
pc_df  <- prcomp(t(logCPM))$x[,1:2]

# P4: 前処理（低カウント→サイズ係数→VST）
keep <- edgeR::filterByExpr(DGEList(cts), group = colData(se)$dex)
dds  <- DESeq2::DESeqDataSet(se[keep,], design = ~ cell + dex) |> DESeq2::estimateSizeFactors()
vsd  <- DESeq2::vst(dds, blind = TRUE)

# P4: 妥当性（RLE中央値≈0）
rle_mat <- SummarizedExperiment::assay(vsd) - matrixStats::rowMedians(SummarizedExperiment::assay(vsd))
boxplot(as.data.frame(rle_mat), outline=FALSE, las=2)

# P5: 再現性（PCA／距離）
DESeq2::plotPCA(vsd, intgroup=c("dex","cell"))
dist_mat <- dist(t(SummarizedExperiment::assay(vsd))); pheatmap::pheatmap(as.matrix(dist_mat))

準備

airway オブジェクトは SummarizedExperiment で提供され、assay(se) に カウント行列、colData(se) に サンプル注釈、rowData(se) に 遺伝子注釈 が入る。

ここで dex を基準化（relevel）するのは、効果の符号と解釈 を一定にするためである。

cell を因子として保持することで、のちの設計行列 ~ cell + dex で ペア差 を適切に差し引けます。

# 補助: サンプル単位QCの簡易要約
qc_sample_summary <- function(counts_mat, topN = 100) {
  stopifnot(is.matrix(counts_mat))
  tibble(
    sample       = colnames(counts_mat),
    libsize      = colSums(counts_mat),
    nonzero_frac = colMeans(counts_mat > 0),
    topN_frac    = apply(counts_mat, 2, function(x){
      x <- x[x > 0]; if(length(x)==0) return(NA_real_)
      s <- sort(x, decreasing = TRUE); sum(head(s, topN))/sum(s)
    })
  )
}
`%||%` <- function(a,b) if (!is.null(a)) a else b

# データ読込（Airway）
data("airway")
se <- airway
colData(se)$dex  <- relevel(factor(colData(se)$dex), "untrt")  # 処理: untrt(基準), trt
colData(se)$cell <- factor(colData(se)$cell)
cts <- SummarizedExperiment::assay(se) |> as.matrix()
head(cts)
##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003        679        448        873        408       1138
## ENSG00000000005          0          0          0          0          0
## ENSG00000000419        467        515        621        365        587
## ENSG00000000457        260        211        263        164        245
## ENSG00000000460         60         55         40         35         78
## ENSG00000000938          0          0          2          0          1
##                 SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003       1047        770        572
## ENSG00000000005          0          0          0
## ENSG00000000419        799        417        508
## ENSG00000000457        331        233        229
## ENSG00000000460         63         76         60
## ENSG00000000938          0          0          0

P1：同定・整合性—サンプルと設計の健全性

解析の出発点は 『誰が誰か』が一意に定まっていることである。重複したサンプル名や、注釈列の命名ゆれ（Run/run/SampleName）は、結合やロング化の段階で サイレントに破綻 を生む。

ここでは列名から sample を生成し、利用可能な列から 実測のラン ID を一貫して採用する。

設計上は ~ cell + dex として ブロック（cell） を明示し、ペア内の差を消した上で dex の主効果を評価できる形に整える。

この段階で、群とブロックに 交絡がない（例：ある cell では trt のみ、などになっていない）ことを必ず確認する。

問題：サンプルの取り違え・ラベル不整合があると以降の解析が破綻する。

目的：設計情報（群・ブロック）がデータと整合しているかを最初に確認する。

ワークフロー：

  1. サンプルIDの重複/欠落チェック → 2) 群×ブロックのバランス確認 → 3)（可能なら）生物指標での照合。 手法：設計表のクロス集計、基本統計（生物指標があれば性別/指標照合）。

  2. 使用パッケージ・関数：SummarizedExperiment::colData、table、duplicated。

# --- P1: 同定・整合性 ------------------------------------------------------
# サンプルIDの一意性
stopifnot(!any(duplicated(colnames(se))))

# 設計の確認：群×ブロック（cell × dex）
#
# Airway データではサンプル名や run の列名がデータセットごとに異なる場合があるため、
# 行名（=サンプル ID）から sample 列を生成し、Run/run/SampleName のいずれかを run 列として採用します。
cdf <- as.data.frame(colData(se))
design_df <- data.frame(
  sample = rownames(cdf),
  run    = cdf$Run %||% cdf$run %||% cdf$SampleName %||% rownames(cdf),
  cell   = cdf$cell,
  dex    = cdf$dex,
  row.names = NULL,
  check.names = FALSE
)
print(design_df)
##       sample        run    cell   dex
## 1 SRR1039508 SRR1039508  N61311 untrt
## 2 SRR1039509 SRR1039509  N61311   trt
## 3 SRR1039512 SRR1039512 N052611 untrt
## 4 SRR1039513 SRR1039513 N052611   trt
## 5 SRR1039516 SRR1039516 N080611 untrt
## 6 SRR1039517 SRR1039517 N080611   trt
## 7 SRR1039520 SRR1039520 N061011 untrt
## 8 SRR1039521 SRR1039521 N061011   trt
print(with(design_df, table(cell, dex)))
##          dex
## cell      untrt trt
##   N052611     1   1
##   N061011     1   1
##   N080611     1   1
##   N61311      1   1

P2 信号品質（測定の健全性：S/N・偏り）

サンプルの健全性を (a) ライブラリサイズ（量）、(b) 非ゼロ比率（広がり）、(c) TopN fraction（集中） の 3 軸で概観する。

  • ライブラリサイズ：極端に小さいサンプルは 有効検出力の低下 と 正規化の不安定 を招きます。群間で著しく偏っていれば要注意である。

  • 非ゼロ比率：多くの遺伝子がゼロのサンプルは、抽出・逆転写・ライブラリ化のトラブルや 過度のフィルタ を疑う。

  • TopN fraction：上位 100 遺伝子の寄与が高すぎる場合、一本勝ちの遺伝子 に総量が吸われ、変動解析検定が不安定になる。

RNA‑seq（バルク）ではミトコンドリア比は参考程度だが、scRNA‑seq では 高ミトコンドリア比＝細胞破砕 の指標として強い意味を持つ。しきい値はデータ依存だが、群間で系統的な差 があるかどうかに注目すること。

問題：読み取りの質や不要断片、増幅の偏りが高いと、下流の変動解析に系統的な歪みが入る。

目的：サンプルごとの「総量」「非欠損」「上位特徴への集中」「（参考）ミトコンドリア比」を把握し、明らかな不良を検出する。

ワークフロー：

  1. サンプル別の総量・非欠損率・トップ集中を算出 → 2) 参考指標（ミト比など） → 3) 可視化で外れ・偏りを見る。

手法：ライブラリサイズ、非ゼロ率、上位寄与率（TopN fraction）。

使用パッケージ・関数：qc_sample_summary（自作）、ggplot2、ggrepel::geom_text_repel。

# --- P2: 信号品質 ----------------------------------------------------------
qc_df <- qc_sample_summary(cts) %>%
  left_join(as.data.frame(colData(se)) %>% tibble::rownames_to_column("sample"),
            by = "sample")

# 参考: ミトコンドリア遺伝子比（Airwayにはgene symbolが無いことが多いのでNA対策）
sym <- rowData(se)$symbol %||% rowData(se)$gene_name %||% rep(NA_character_, nrow(se))
mt_rows <- if (all(is.na(sym))) rep(FALSE, nrow(se)) else grepl("^MT-", sym)
qc_df$mito_frac <- if (any(mt_rows)) colSums(cts[mt_rows,,drop=FALSE]) / qc_df$libsize else NA_real_

# 可視化1: ライブラリサイズ
ggplot(qc_df, aes(x = reorder(sample, libsize), y = libsize/1e6, fill = dex)) +
  geom_col() + coord_flip() +
  labs(x=NULL, y="Library size (million reads)", title="ライブラリサイズ（サンプル別）")

# 可視化2: 非欠損 × 上位集中
ggplot(qc_df, aes(nonzero_frac, topN_frac, color = dex, shape = cell, label = sample)) +
  geom_point(size=3) +
  ggrepel::geom_text_repel(show.legend = FALSE, max.overlaps = 20) +
  labs(x="非ゼロ比率", y="上位100遺伝子の寄与（Top100 fraction）",
       title="信号の密度と上位遺伝子への集中")

# 可視化3: （あれば）ミトコンドリア比
if (!all(is.na(qc_df$mito_frac))) {
  ggplot(qc_df, aes(x = dex, y = mito_frac, fill = dex)) +
    geom_boxplot(width=.6, alpha=.7) + geom_jitter(width=.1, height=0, size=2) +
    labs(x=NULL, y="ミトコンドリア比", title="ミトコンドリア由来リードの割合（参考）")
}

P3 系統的偏り（設計要因・構成バイアス）

粗い正規化（logCPM）で PCA を行い、PC1 が log10(library size) と強く相関しているかを確認する。

高相関であれば 規模効果 が支配している可能性が高く、下流のクラスタや変動解析にも影響する。

この段階で対処できるのは、(i) 適切な正規化（サイズ係数/TMM）、(ii) 分散安定化（VST/voom）、(iii) 必要なら バッチ補正（sva, ComBat など）である。ただし、群と完全に交絡した『バッチ』は 推定不能 で、補正してはならない。

問題：測定日や装置などの設計要因、または総量差に伴う相対量の歪みが、主なばらつき要因になると、生物学的差異が埋もれる。

目的：主要成分（PCA）で「群」ではなく「バッチ/規模」が支配していないかを事前診断する。

ワークフロー：

  1. 粗い正規化（logCPM）→ 2) PCA で低次元表示 → 3) PC1 と総量の相関を評価。

手法：logCPM、PCA、PC1–library size 相関。

使用パッケージ・関数：edgeR::cpm、stats::prcomp、ggplot2。

# --- P3: 系統的偏りの兆候 ---------------------------------------------------
logCPM <- edgeR::cpm(cts, log=TRUE, prior.count = 1)
pca0   <- prcomp(t(logCPM), scale. = FALSE)

pc_df  <- tibble::as_tibble(pca0$x[,1:2], rownames = "sample") %>%
  dplyr::left_join(
    dplyr::select(qc_df, sample, libsize, dex, cell),
    by = "sample"
  )

ggplot(pc_df, aes(PC1, PC2, color = dex, shape = cell, size = libsize)) +
  geom_point() + scale_size_continuous(range = c(2.5,6)) +
  labs(title="PCA（粗い正規化：logCPM）", subtitle="規模・ブロックの影響を事前確認")

cor_PC1_lib <- cor(pc_df$PC1, log10(pc_df$libsize))
message(sprintf("PC1 と log10(library size) の相関: %.3f", cor_PC1_lib))

P4 正規化・分散安定化（ライブラリ補正 → VST）

edgeR::filterByExpr() は設計（群ラベル）を尊重して 検出可能性の低い遺伝子 を除去し、偽陽性率を制御しつつ検出力を上げる。

その後、DESeq2::estimateSizeFactors() により RLE ベースのサイズ係数 を推定し、総量差と構成差を均すことで、VST（vst(blind=TRUE)）で 平均–分散依存 を弱める。

RLE の箱ひげは 中央値≈0 が目安、PCA は群とブロックの分離が妥当か、サンプル距離ヒートマップは 群内高い一貫性・外れの有無 を点検する。

問題：サンプル間の総量差と、平均が大きいほど分散も大きくなる平均–分散依存。

目的：まずサイズ係数で総量を合わせ、次に分散安定化でばらつきを均質化し、比較可能な空間を作る。

ワークフロー：

  1. 低カウント除去（群設計を考慮）→ 2) サイズ係数（RLE）→ 3) 分散安定化（VST）→ 4) RLE/可視化で妥当性確認。

手法：filterByExpr、DESeq2::estimateSizeFactors、DESeq2::vst、RLEプロット、PCA。

使用パッケージ・関数：edgeR、DESeq2、matrixStats::rowMedians、DESeq2::plotPCA。

# --- P4: 正規化・分散安定化 -----------------------------------------------
# 低カウント除去（群デザインを尊重）
y <- DGEList(counts = cts, group = colData(se)$dex)
keep <- edgeR::filterByExpr(y, group = colData(se)$dex)
se_f <- se[keep, ]

# DESeq2データセット & サイズ係数（RLE）
dds <- DESeqDataSet(se_f, design = ~ cell + dex)
dds <- estimateSizeFactors(dds)
sf  <- sizeFactors(dds); print(sf)
## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 
##  1.0219051  0.8963516  1.1732398  0.6692888  1.1744015  1.3995153  0.9163097 
## SRR1039521 
##  0.9452228
# VST（可視化・距離・クラスタ安定化に有用）
vsd <- vst(dds, blind = TRUE)

# RLEプロット（理想は各箱の中央値≈0）
rle_mat <- assay(vsd) - matrixStats::rowMedians(assay(vsd))
boxplot(as.data.frame(rle_mat), outline = FALSE, las = 2,
        main = "RLEプロット（分散安定化後の中央値差）", ylab = "相対対数発現")

# PCA（群・ブロック・過補正の点検）
DESeq2::plotPCA(vsd, intgroup = c("dex","cell"))

# サンプル間距離（外れ・群内一貫性）
dist_mat <- dist(t(assay(vsd)))
pheatmap::pheatmap(as.matrix(dist_mat),
                   annotation_col = as.data.frame(colData(se_f)[,c("dex","cell")]))

P5 再現性・外れ値の点検 ＋ 変動解析（効果量縮約込み）

DESeq() は NB GLM を適合し、results() は FDR（Benjamini–Hochberg） を既定で返す。

lfcShrink(type="apeglm") により 効果量（log2FC）の縮約 を行うと、小標本での過大推定を抑え、ランキングが安定する。

同時に、assays(dds)[["cooks"]] から Cook の距離 を点検し、特定サンプルが推定に過大な影響 を与えていないかを確認する。最終報告では、p 値・FDR に加えて 効果量とその不確実性（区間推定） を主役に据えるのが推奨。

問題：技術的ばらつきや外れ値が変動解析検出を不安定にし、推定のばらつき（とくに小標本・低発現）を過大評価しがち。

目的：距離・相関・影響度で再現性と外れを点検し、縮約推定を用いて効果量のばらつきを抑えた変動解析を行う。

ワークフロー：

  1. 本解析（分散の借用・回帰）→ 2) 効果量縮約 → 3) MAと上位表 → 4) 影響度（Cook’s距離）の概観。

手法：DESeq2::DESeq、lfcShrink(type="apeglm")、MAプロット、Cook’s距離。

使用パッケージ・関数：DESeq2、apeglm。

# --- P5: 再現性・外れ点検 + 変動解析 --------------------------------------
dds <- DESeq(dds)

# 変動解析（trt vs untrt）
res <- results(dds, contrast = c("dex","trt","untrt"))

# 効果量縮約（apeglm）
coef_name <- grep("dex.*trt.*vs.*untrt", resultsNames(dds), value = TRUE)[1]
res_shr   <- lfcShrink(dds, coef = coef_name, type = "apeglm")

summary(res_shr)
## 
## out of 15926 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 2589, 16%
## LFC < 0 (down)     : 2306, 14%
## outliers [1]       : 0, 0%
## low counts [2]     : 618, 3.9%
## (mean count < 13)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
plotMA(res_shr, main = "MA plot (lfcShrink: apeglm)", ylim = c(-4,4))

# 上位ヒット
top_tbl <- as_tibble(res_shr, rownames = "gene_id") %>%
  arrange(padj) %>% filter(!is.na(padj)) %>% slice(1:20)
print(top_tbl)
## # A tibble: 20 × 6
##    gene_id         baseMean log2FoldChange  lfcSE    pvalue      padj
##    <chr>              <dbl>          <dbl>  <dbl>     <dbl>     <dbl>
##  1 ENSG00000152583     998.           4.55 0.186  1.14e-135 8.91e-132
##  2 ENSG00000165995     495.           3.27 0.133  1.16e-135 8.91e-132
##  3 ENSG00000101347   12706.           3.75 0.158  4.81e-128 1.84e-124
##  4 ENSG00000120129    3411.           2.93 0.123  4.54e-128 1.84e-124
##  5 ENSG00000189221    2343.           3.34 0.143  1.44e-122 4.42e-119
##  6 ENSG00000211445   12291.           3.71 0.169  4.98e-110 1.27e-106
##  7 ENSG00000157214    3012.           1.96 0.0908 5.98e-105 1.31e-101
##  8 ENSG00000162614    5397.           2.02 0.0959 4.57e-100 8.74e- 97
##  9 ENSG00000125148    3658.           2.20 0.107  5.06e- 95 8.60e- 92
## 10 ENSG00000154734   30335.           2.33 0.118  2.68e- 88 4.10e- 85
## 11 ENSG00000139132    1224.           2.21 0.113  1.99e- 86 2.77e- 83
## 12 ENSG00000162493    1100.           1.88 0.0963 9.54e- 86 1.22e- 82
## 13 ENSG00000162692     510.          -3.68 0.192  1.31e- 84 1.54e- 81
## 14 ENSG00000179094     777.           3.17 0.165  2.04e- 84 2.23e- 81
## 15 ENSG00000134243    5515.           2.18 0.113  3.23e- 84 3.30e- 81
## 16 ENSG00000163884     561.           4.43 0.237  1.03e- 80 9.82e- 78
## 17 ENSG00000178695    2658.          -2.52 0.137  7.54e- 78 6.79e- 75
## 18 ENSG00000148848    1369.          -1.81 0.101  8.17e- 73 6.95e- 70
## 19 ENSG00000146250     315.          -2.75 0.154  2.11e- 72 1.70e- 69
## 20 ENSG00000107562   25203.          -1.91 0.107  2.66e- 72 2.04e- 69
# 影響度（Cook's距離 99%分位・サンプル別）
ck <- assays(dds)[["cooks"]]
ck_max <- apply(ck, 2, function(v) quantile(v, 0.99, na.rm=TRUE))
barplot(ck_max, las = 2, main = "Cook's distance (99%分位)", ylab = "Cooks 99th percentile")

付録：QC概念を「マルチ」にとらえて理解を深めていこう

QCそのものの原理原則をおさえて本質を捉えることで、それらを一般化して、様々な階層のオミクスへと拡張して、展開していくスキルを身につけていくことも大事である。ここでは、他のオミクス階層の場合のQC・前処理の概要を示す。各自調べて、RNA-seqと同様に具体的課題と解決方法の原理、実装方法を含めて、理解するように。

質量分析（プロテオーム／メタボローム）

段階 原理（なぜ） データ現象 検出・管理の要点
P1 同定 ラベル・参照・標準の誤扱い 参照比逸脱、標準不回収 割付表照合、参照比・回収率の監査
P2 信号 電離効率低下・時間ずれ・校正ずれ 総量の滑り、保持時間シフト、同定数低下 確認用試料の時系列監視と再校正
P3 かたより ドリフト・残留混入・日差 時系列トレンド、空試料ピーク、バッチ分離 時間補正、空試料監視、残差補正
P4 正規化 注入量差・形状差・比圧縮 スケール不一致、右裾、比縮小 内部基準・分位合わせ・分散安定化
P5 再現 検出・整列の不安定、欠測機構 変動増大、群偏在の欠測 確認用CV・相関で判定、機構に応じた補完

P1 同定・整合性

原理：同時定量用ラベルの割付や内部標準・参照測定の取り扱いに誤りが生じる。

データ現象：参照比の逸脱、内部標準の回収不良。

検出・管理：ラベル割付表・参照比・内部標準の回収率を照合し、整合しない試料は再点検する。

P2 信号品質

原理：電気噴霧の効率低下（基質効果）、保持時間のずれ（カラム劣化や温度変動）、質量校正のずれ、確率サンプリングによる取りこぼし。

データ現象：総イオン量の時間的な滑り、保持時間のシフト、質量誤差の増加、同定数の減少、ピーク形状の崩れ。

検出・管理：確認用混合試料を列に挿入して時系列で監視し、必要に応じて校正・部品交換・測定条件の最適化を行う。

P3 系統的偏り・汚染

原理：測定順に沿った機器ドリフト、前試料の残留（キャリーオーバー）、日・装置・オペレーター差。

データ現象：確認用試料の強度が時間とともに一方向に変化、空試料にも特定のピークが出現、主成分空間でバッチごとに分離。

検出・管理：確認用試料にもとづく時間ドリフト補正、空試料の監視、残存するバッチ差の統計的補正。測定順は計画段階で分散配置する。

P4 正規化・較正

原理：注入量やイオン化効率の差が全体量の差を生み、強度分布の形状もサンプル間で異なる。ラベル法では共分離により比が縮むことがある。

データ現象：スケールの不一致、右に裾を引く分布、比の縮小。

検出・管理：内部標準や参照チャンネル、分位合わせや分散安定化により、スケールと形状を統一する。ラベル法では参照比にもとづく補正で比の安定性を高める。

P5 再現性・外れ値

原理：ピーク検出・整列の不安定、検出限界付近の欠測（左側の打ち切り）。

データ現象：繰り返しの変動増大、群ごとの欠測の偏り。

検出・管理：確認用試料のばらつき（変動係数）や繰り返し相関で合否を判定し、欠測は発生機構（偶然か検出限界か）に応じて適切な補完法を用いる。

シングルセル

段階 原理（なぜ） データ現象 検出・管理の要点
P1 同定 二重封入・ラベル交差 分子数過多、相反指標の同時高値 二重検出・ラベル重複の同定と除去
P2 信号 破損・遊離分子の混入 小器官比の上昇、量の極端 品質指標で外れ除去、混入の差し引き
P3 かたより 試薬・日付・周期差 バッチが主分離軸 量合わせ後のデータ統合と事前後比較
P4 正規化 捕捉効率差・ヘテロスケ スケール不一致、強発現の不均一 量合わせ＋確率モデルによる均質化
P5 再現 希少型の取りこぼし等 構成の揺れ、指標の不一致 クラスタ安定性と擬似バルク整合で確認

P1 同定・整合性

原理：一つの容器に複数細胞が同封される、または混合実験のラベルが交差する。

データ現象：単一細胞としては不自然に多い分子数、異なる細胞系譜の指標が同時に高い。

検出・管理：二重封入の確率モデルやラベル重複の検査で同定し、除去する。遺伝型情報があれば分離推定を行う。

P2 信号品質

原理：細胞破損・分解により細胞小器官由来の読み取りが増え、周囲の遊離分子が混入する。

データ現象：小器官関連の読み取り比率の上昇、分子数・遺伝子数の極端な外れ。

検出・管理：品質指標の分布から外れ細胞を除去し、遊離分子の寄与を差し引く。

P3 系統的偏り・汚染

原理：試薬・キット・日付の違い、細胞周期の位相差。

データ現象：低次元空間でバッチが主な分離軸になる。

検出・管理：量の補正の後、データ統合手法でバッチ差を弱め、統合前後で既知マーカーの一貫性と混合度を確認する。

P4 正規化・較正

原理：捕捉効率の差により細胞間の総分子数が異なり、強度が大きいほどばらつきも大きくなる。

データ現象：スケールの不一致、強発現遺伝子での不均一な変動。

検出・管理：細胞ごとの量合わせの後、確率モデルにもとづく分散の均質化を行い、安定に変動する遺伝子を選んで表現空間を構築する。

P5 再現性・外れ値

原理：希少細胞の取りこぼし、過度なクラスタ分割。

データ現象：測定単位ごとに細胞型構成が揺れる、同一細胞型の指標が再現しない。

検出・管理：クラスタの安定性指標や「擬似バルク」（細胞型合算）と外部参照との整合で妥当性を確認する。

メチローム（チップ／全ゲノム）

段階 原理（なぜ） データ現象 検出・管理の要点
P1 同定 取り違え 性別・指標不一致 性別推定・遺伝的指標で照合
P2 信号 結合・色素・深さの不足 非検出点・弱信号の増加 検出性指標での選別
P3 かたより 位置・日・型の差／配列特性 プレート集団化／配列依存の偏り 専用補正・共変量調整
P4 正規化 背景・型差／深さ差 分布の不揃い 背景・型補正／平滑化と重み付け
P5 再現 技術ばらつき 反復の不一致 相関・既知領域の方向性で確認

アレイ（450k/EPIC）

  • P1 原理：取り違え（sample swap）。

    データ現象：性別推定・SNPプローブの不一致（getSex, SNP probes）。

    検出・管理：性別・SNP照合（minfi, ewastools）。

  • P2 原理：結合不良・色素差・低ビーズ数（hybridization, dye bias, beadcount）。

    データ現象：検出p高値・弱信号（detection P）。

    検出・管理：品質レポートと不良プローブ除外（qcReport, detectionP）。

  • P3 原理：プレート/日/プローブ型差・細胞組成差（plate effect, Type I/II, cell composition）。

    データ現象：PCAでプレート分離。

    検出・管理：背景・色素・型補正と組成推定（Noob/FunNorm, BMIQ, estimateCellCounts, ComBat）。

  • P4 原理：スケール差・型差。

    データ現象：β値/M値分布の不揃い（beta/M values）。

    検出・管理：補正後の尺度で解析（M-value modeling）。

  • P5 原理：技術ばらつき。

    データ現象：反復の相関低下。

    検出・管理：反復相関・既知CpGの方向性（smoking CpGs）。

WGBS/oxBS

  • P1：取り違え（SNP/sex chromosome depth）。

  • P2：変換効率・深さ（bisulfite conversion rate, coverage）。

  • P3：配列特性・方向バイアス（GC bias, strand bias）。

  • P4：平滑化・重み付け（BSmooth, depth weighting）。

  • P5：領域差の再現（DMR reproducibility）。

ChIP-seq / ATAC-seq

段階 原理（なぜ） データ現象 検出・管理の要点
P1 同定 標的・対照の誤設定 既知標的で信号不在 既知領域での対照比較
P2 信号 酵素活性・抗体特異性の不足 開始近傍の集積低下（TSS）、ピーク内割合低（FRiP）、断片長の乱れ 主要指標の定量と閾値管理
P3 かたより ドリフト・残留・配列特性 時系列滑り、空試料ピーク、配列依存 時間補正、空試料監視、残差補正
P4 正規化 深さ差／疎な行列 高さ不一致、低深さで消失 深さスケーリング、単一細胞は重み付け後に低次元化
P5 再現 検出の不安定・標的の当たり外れ 独立測定で重なり低 再現性指標（重なり率・IDR 等）で判定

P1 同定・整合性

  • 原理：標的・対照の誤設定（target antibody, input/IgG control）。

  • データ現象：既知標的領域でシグナル不在。

  • 検出・管理：既知座の検証（positive loci）。

P2 信号品質

  • 原理：切断酵素活性・抗体特異性の不足（Tn5 activity, antibody specificity）。

  • データ現象：開始点近傍の集積低下、ピーク内割合低、断片長周期性の崩れ、ミト過多（TSS enrichment, FRiP, fragment length periodicity, mitochondrial fraction）。

  • 検出・管理：主要指標の定量と閾値管理（ATACseqQC, ChIPQC）。

P3 系統的偏り・汚染

  • 原理：測定順ドリフト・残留・配列特性（run-order drift, carryover, mappability/GC）。

  • データ現象：時系列トレンド、空試料ピーク、GCと強度の相関。

  • 検出・管理：時間補正・空試料監視・統計補正（LOESS, blank, ComBat）。

P4 正規化・較正

  • 原理：読み取り深さ差・疎行列（library depth, sparsity）。

  • データ現象：高さ不一致、低深さで消失。

  • 検出・管理：深さスケーリング／単一細胞は重み付け後の低次元化（CPM/RPKM/TMM in peaks, TF-IDF→LSI）。

P5 再現性・外れ値

  • 原理：ピーク検出の不安定・標的の当たり外れ（peak caller variability）。

  • データ現象：独立測定での重なり低下。

  • 検出・管理：再現性評価・相関・クラスタ安定性（IDR, DiffBind）。