ngs-quality-control

NGS 品質管理

概要

このスキルでは、Illumina ショートリードシーケンシングデータの全ライフサイクル（生リード → トリミング済みリード → アライン済み BAM）にわたる、再現性のある QC ワークフローを実行します。WGS、WES、ターゲットパネルのプロジェクトで QC 判定に利用できる、決定論的なメトリクスを生成します。

パイプラインの各段階：

1. 生リード QC — FASTQ ごとに FastQC を実行し、MultiQC で集約します。
2. アダプター／品質トリミング — ペアエンドのアダプター検出を有効にして fastp を実行します。
3. トリミング後の QC — トリミング済みリードに対して FastQC と MultiQC を実行します。
4. アラインメント QC — samtools stats／flagstat、Picard メトリクススイートを使用します。
5. カバレッジ QC — mosdepth（ゲノムウィンドウまたはターゲット BED）を使用します。
6. 集約 — すべてのツール出力を MultiQC で集約し、単一の HTML サマリーを作成します。

使用方法

THREADS=8、OUT=qc_out、ペア FASTQ sample_R1.fastq.gz／sample_R2.fastq.gz、インデックス付きの座標ソート済み sample.bam、リファレンス ref.fasta、および（WES／パネルの場合は）ターゲット区間があることを前提とします。

1. 生 FASTQ の QC — FastQC

mkdir -p $OUT/fastqc_raw
fastqc -t $THREADS -o $OUT/fastqc_raw *.fastq.gz

確認すべき主な FastQC モジュール：

• 塩基位置ごとの配列品質 — 3′末端に向かって Q30 が低下している場合、トリミングが必要です。
• 配列ごとの品質スコア — 二峰性分布は、タイル／レーンの不良を示します。
• 塩基位置ごとの配列組成 — 最初の約 12 bp の偏りは正常です（ランダムヘキサマープライミング）。それ以外の大きな偏りは、アダプターやコンタミネーションを示します。
• 配列ごとの GC 含有量 — 予想される生物種の GC 含有量からの逸脱は、コンタミネーションまたはバイアスを示します。
• 配列重複レベル — アンプリコンライブラリ以外で重複率が高い場合、ライブラリ複雑度に問題があります。
• 過剰出現配列 — アダプターのリードスルー、rRNA、PhiX のキャリーオーバーを示します。
• アダプター含有量 — トリミング対象のアダプターセットを確認できます。

2. アダプターと品質のトリミング — fastp

mkdir -p $OUT/fastp
fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o $OUT/fastp/sample_R1.trim.fastq.gz -O $OUT/fastp/sample_R2.trim.fastq.gz \
  --detect_adapter_for_pe \
  --thread $THREADS \
  --html $OUT/fastp/sample.fastp.html \
  --json $OUT/fastp/sample.fastp.json

便利な追加フラグ：

• --qualified_quality_phred 20 — 塩基品質のしきい値（デフォルトは 15）です。
• --length_required 50 — トリミング後に 50 bp 未満のリードを除外します。
• --cut_tail --cut_tail_window_size 4 --cut_tail_mean_quality 20 — スライディングウィンドウによる 3′品質トリミングです。
• --dedup — メモリ内で PCR 重複を除去します（リードペアの同一性に基づきます。多くのワークフローでは、アラインメント後に Picard／MarkDuplicates を使う方法を推奨します）。

3. トリミング後の QC と MultiQC

mkdir -p $OUT/fastqc_trim
fastqc -t $THREADS -o $OUT/fastqc_trim $OUT/fastp/*.trim.fastq.gz

multiqc $OUT -o $OUT/multiqc

MultiQC は、スキャン対象ディレクトリ内の FastQC、fastp、Picard、samtools、mosdepth の出力を自動検出します。

4. アラインメント QC — samtools

samtools flagstat -@ $THREADS sample.bam > $OUT/samtools/sample.flagstat.txt
samtools stats    -@ $THREADS sample.bam > $OUT/samtools/sample.stats.txt
samtools idxstats sample.bam               > $OUT/samtools/sample.idxstats.txt

samtools stats から取得する項目：raw total sequences、reads mapped、error rate、insert size average、insert size standard deviation。

5. Picard メトリクススイート

Picard を使うには、ref.fasta と同じ場所にシーケンス辞書（ref.dict。picard CreateSequenceDictionary R=ref.fasta で作成）と、ハイブリッドキャプチャメトリクス用の Picard 形式の区間リスト（.interval_list、ヘッダーは BAM の @SQ 行。単なる BED ではありません）が必要です。

# Convert BED → interval_list (required for CollectHsMetrics)
picard BedToIntervalList \
  I=targets.bed O=targets.interval_list SD=ref.dict

# General alignment metrics (all library types)
picard CollectAlignmentSummaryMetrics \
  R=ref.fasta I=sample.bam O=$OUT/picard/sample.alignment_metrics.txt

# Insert size distribution (paired-end)
picard CollectInsertSizeMetrics \
  I=sample.bam O=$OUT/picard/sample.insert_size.txt \
  H=$OUT/picard/sample.insert_size.pdf

# WGS coverage metrics
picard CollectWgsMetrics \
  R=ref.fasta I=sample.bam O=$OUT/picard/sample.wgs_metrics.txt

# WES / targeted panel hybrid-selection metrics
picard CollectHsMetrics \
  R=ref.fasta I=sample.bam O=$OUT/picard/sample.hs_metrics.txt \
  BAIT_INTERVALS=baits.interval_list \
  TARGET_INTERVALS=targets.interval_list

6. カバレッジ — mosdepth

# WGS: 500 bp fixed windows
mosdepth --by 500 --threads $THREADS --fast-mode --no-per-base \
  $OUT/mosdepth/sample_wgs sample.bam

# WES / panel: per-target coverage
mosdepth --by targets.bed --threads $THREADS --no-per-base \
  $OUT/mosdepth/sample_wes sample.bam

注目すべき出力：

• *.mosdepth.summary.txt — 染色体ごとおよび全体の平均深度です。
• *.mosdepth.region.dist.txt — カバレッジが X 以上の領域の累積割合です（% targets ≥20x に使用します）。
• *.regions.bed.gz — ウィンドウ／ターゲットごとの平均深度です。

層別カバレッジビンを出力するには --quantize 0:1:10:30: を、指定した深度で領域ごとの合格数を出力するには --thresholds 10,20,30,50 を追加します。

応答形式

• まずユーザーが必要とするコマンドまたはコードを提示し、その後に説明します。
• 次の構成にします：入力の確認 → 動作するコード → 主要パラメーターの説明 → 注意点
• 各タスクにつき、動作する完全な例を 1 つ示します。あらゆる代替案を並べないでください。
• コードコメントは最小限かつ実用的なものにします（存在理由ではなく、何をするかを記述します）。
• 目標：前後に簡潔な説明を添えた 50〜100 行のコード

基本概念

シーケンス種別ごとのQCツールの選択

シーケンス種別	カバレッジツール	カバレッジ指標	Picard スイート	主な閾値
WGS	mosdepth --by 500	CollectWgsMetrics	PCT_20X ≥ 95%	平均 ≥ 30×
WES	mosdepth --by targets.bed	CollectHsMetrics	PCT_TARGET_BASES_20X ≥ 80%	オンターゲット ≥ 60%
ターゲットパネル	mosdepth --by targets.bed	CollectHsMetrics	PCT_TARGET_BASES_20X ≥ 95%	平均 ≥ 500×（体細胞）
RNA-seq	—（featureCounts/STAR のログを使用）	CollectRnaSeqMetrics	PCT_CODING_BASES	rRNA < 10%

すべての種別： FastQC（トリミング前＋後）→ fastp（トリミング）→ samtools flagstat → MultiQC（集約）。

クオリティスコアのエンコーディング検出

最新の Illumina（CASAVA ≥1.8、2011年以降）では Phred+33（Sanger）が出力されます。旧形式のデータ（Illumina 1.3–1.7、Solexa）では Phred+64 が使用されていました。アライナーでは Phred+33 が想定されています。

• FastQC では Basic Statistics → Encoding にエンコーディングが表示されます。Sanger / Illumina 1.9 を確認してください。

• 簡単な CLI チェック — クオリティ行の ASCII 範囲を調べます：

  zcat sample_R1.fastq.gz | head -4000 \
    | awk 'NR%4==0' \
    | od -An -c | tr -s ' \n' '\n' \
    | awk 'length==1' | sort -u
  

  !–I の範囲の文字 → Phred+33、@–h の文字 → Phred+64。

• 旧形式のデータを変換します：seqtk seq -Q64 -V in.fq > out.fq（+64 → +33 にシフト）。

トリミングパラメーターの意図

• アダプタートリミングは、アダプターまで読み込んでいる場合（インサート長 < リード長）に行います。--detect_adapter_for_pe はリードの重なりを利用するため、ペアエンドでは配列ベースの検出より信頼性が高くなります。
• クオリティトリミングは、2チャネルのケミストリーにおける 3′ 末端のクオリティ低下に対して行います（NextSeq/NovaSeq では暗サイクルで poly-G が出力されます。これらのプラットフォームでは --trim_poly_g を使用してください。fastp はプラットフォーム検出によって自動的に有効にします）。
• ライブラリ設計（UMI、インラインバーコード）が分かっていない限り、固定長でハードクリップしないでください。

重複：光学的重複と PCR 重複

• samtools flagstat の「duplicates」は BAM に設定されたフラグを反映します（通常は Picard MarkDuplicates または samtools markdup によって設定）。MarkDuplicates はトリミング済み FASTQ ではなく、アラインメントの後に実行してください。
• CollectWgsMetrics / CollectHsMetrics は PCT_EXC_DUPE（重複により除外された塩基の割合）を報告します。これはカバレッジベース QC における標準的な重複率です。

カバレッジ：平均値と均一性

平均深度だけでは不十分です。次の両方を報告してください。

• 平均ターゲットカバレッジ（HsMetrics の MEAN_TARGET_COVERAGE、WgsMetrics の MEAN_COVERAGE）。
• カバレッジ幅 — % targets ≥20x / % genome ≥20x（Picard の PCT_TARGET_BASES_20X、PCT_20X、または mosdepth の region.dist.txt）。

平均 200x でも、20x に達するターゲットが 70% のパネルでは、カバレッジ均一性に問題があります（キャプチャ効率、GC バイアス）。

クイックリファレンス

クオリティ閾値（一般的なゲート判定）

指標	WGS	WES	ターゲットパネル	ソース
平均カバレッジ	≥ 30×	オンターゲット ≥ 50×	≥ 500×（体細胞）、≥ 100×（生殖細胞系列）	WgsMetrics / HsMetrics
20× でのカバレッジ幅	≥ 95%（ゲノム）	≥ 80%（ターゲット）	≥ 95%（ターゲット）	Picard PCT_20X / mosdepth
10× でのカバレッジ幅（生殖細胞系列）	≥ 90%	≥ 90%（ターゲット）	≥ 98%（ターゲット）	Picard PCT_10X
重複率	< 30%	< 40%	< 50%（UMI を推奨）	MarkDuplicates / Picard PCT_EXC_DUPE
マッピング率	> 95%	> 95%	> 95%	samtools flagstat / AlignmentSummaryMetrics
平均塩基クオリティ	≥ 80% の塩基で Q30 以上	同左	同左	FastQC / fastp
インサートサイズ（PE150）	300–500 bp、SD < 平均の 20%	150–300 bp	ライブラリ固有	CollectInsertSizeMetrics
オンターゲット率（キャプチャ）	該当なし	≥ 60%（PCT_SELECTED_BASES）	≥ 70%	CollectHsMetrics
コンタミネーション	< 1%	< 1%	< 1%	VerifyBamID（このスキルの対象外）

ワンライナー

# Full raw → trimmed → FastQC → MultiQC
fastqc -t 8 -o qc/raw *.fastq.gz && \
fastp -i R1.fq.gz -I R2.fq.gz -o R1.trim.fq.gz -O R2.trim.fq.gz \
      --detect_adapter_for_pe --thread 8 --html fastp.html --json fastp.json && \
fastqc -t 8 -o qc/trim *.trim.fq.gz && \
multiqc qc -o qc/multiqc

# BAM QC bundle
samtools flagstat sample.bam > qc/sample.flagstat.txt
picard CollectAlignmentSummaryMetrics R=ref.fasta I=sample.bam O=qc/aln.txt
picard CollectInsertSizeMetrics I=sample.bam O=qc/ins.txt H=qc/ins.pdf
mosdepth --by 500 --threads 8 --fast-mode --no-per-base qc/sample sample.bam

ツールの役割一覧

ツール	入力	主な出力	用途
FastQC	FASTQ	ファイルごとのHTML/zip	トリミング前後のリード品質
fastp	FASTQペア	トリミング済みFASTQ + JSON/HTML	アダプター除去 + 品質トリミング
samtools stats/flagstat	BAM	テキスト	マッピング数、エラー率
Picard AlignmentSummary	BAM	テキスト	鎖のバランス、PFリード、Q20/Q30のアラインメント済み塩基
Picard InsertSize	BAM	テキスト + PDF	PEインサートサイズ分布
Picard WgsMetrics	BAM	テキスト	WGSカバレッジ + カバー率
Picard HsMetrics	BAM	テキスト	WES/パネルのオンターゲット率とカバレッジ
mosdepth	BAM	BED + dist	高速なウィンドウ別／ターゲット別深度
MultiQC	ディレクトリ	単一のHTML	コホートの集約

よくある間違い

• 誤り： トリミング後にのみFastQCを実行する 正しい方法： 生リードとトリミング済みリードの両方に対してFastQCを実行する 理由： トリミング前のQCでライブラリを診断し、トリミング後のQCで過剰なクリッピングをせずにトリミングできたことを確認できます

• 誤り： 生のBEDファイルをPicard CollectHsMetricsに渡す 正しい方法： picard BedToIntervalList I=x.bed O=x.interval_list SD=ref.dict でPicard形式の.interval_listに変換する 理由： BAIT_INTERVALS/TARGET_INTERVALSには、完全な@SQヘッダーを含む.interval_listが必要です。生のBEDを渡すと失敗するか、誤ったメトリクスが警告なしに生成されることがあります

• 誤り： コンティグ名やソート順が一致しないインターバルリストファイルを使用する 正しい方法： コンティグ名（chr1と1など）と@SQの順序がBAMヘッダーと完全に一致することを確認する 理由： 不一致があると、メトリクスが警告なしに誤った値になったり、ツールが失敗したりします

• 誤り： トリミング済みFASTQまたは名前順ソート済みBAMに対してMarkDuplicatesを実行する 正しい方法： アラインメント済みの座標順ソート済みBAMに対してのみMarkDuplicatesを実行する 理由： MarkDuplicatesでは、光学的重複とPCR重複を特定するために、アラインメント済みリードが座標順にソートされている必要があります

• 誤り： QCメトリクスとして平均カバレッジのみを報告する 正しい方法： 平均深度にカバー率（% ≥20x）と重複率を必ず併記する 理由： 均一性が低い場合、高い平均値によって失敗しているターゲットが見えなくなります

• 誤り： 2012年以前のFASTQデータで自動エンコーディング検出を信用する 正しい方法： 古いデータでは、アラインメント前にPhredエンコーディングを明示的に確認する 理由： Phred+64データを現代のアライナーに渡すと、塩基品質値が誤ったまま警告なしに処理されます

• 誤り： レーンごとのQCを実行する前に、レーン／サンプルをマージする 正しい方法： まずレーン／FASTQごとにFastQCを実行し、レーンごとのQCに合格した後にのみマージする 理由： マージ後は、問題のあるレーンを特定できなくなります

• 誤り： fastpの--dedupを主な重複除去ステップとして使用する 正しい方法： アラインメント後の重複除去にはPicard MarkDuplicatesまたはsamtools markdupを使用する 理由： fastpの--dedupはアラインメント前に完全一致する配列の重複のみを検出するため、アラインメントに基づく重複を見逃します

• 誤り： --no-per-baseを指定せずに大きなゲノムでmosdepthを実行する 正しい方法： 大きなゲノムのQCでは、常に--no-per-baseを使用する 理由： 塩基ごとの出力は膨大になり、必要になることはほとんどありません。QCにはウィンドウ別／領域別の出力で十分です

参考資料

• FastQC：https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
• fastp：Chen et al. Bioinformatics 2018, https://doi.org/10.1093/bioinformatics/bty560
• Picard：https://broadinstitute.github.io/picard/
• mosdepth：Pedersen & Quinlan, Bioinformatics 2018, https://doi.org/10.1093/bioinformatics/btx699
• カバレッジに関する推奨事項：Sims et al. Nat Rev Genet 2014, https://doi.org/10.1038/nrg3642