第4章: ゲノム解析技術

4. ゲノム解析技術

この章を読む理由

ゲノム解析は、FASTQ から BAM / VCF へ至る成果物の流れを最初に固定する章です。 第5章以降の発現解析や臨床応用を理解するためにも、 「どこで品質が落ちるか」「どの成果物を次工程へ渡すか」をここで整理しておく必要があります。

本章で使う共通題材

• 題材A: SARS-CoV-2 公開データ
  • ラン: SRR11140744
  • 参照配列: MN908947.3
  • 本章では、最小の QC → マッピング → 変異検出の流れを確認するために使います。

学習目標

• NGS 解析パイプライン（QC→アライメント→後処理→変異検出→アノテーション→要約）と、代表的な入出力（FASTQ/BAM/VCF）を説明できる
• QC / マッピング / 変異検出で確認すべき品質指標と、その確認手段（FastQC/MultiQC、samtools flagstat、bcftools stats 等）を説明できる
• 参照ゲノム、パラメータ、ツール / DB のバージョン、ログ記録をどう残せば第三者が再実行できるかを整理できる

前提知識

• 第3章の k-mer、seed-and-extend、BWT / FM-index、minimizer、de Bruijn graph の考え方。ツール名よりも、どの処理が「探索」「整列」「圧縮」「グラフ化」に相当するかを説明できること。
• FASTQ の read / quality score / paired-end、FASTA の参照配列、SAM/BAM/CRAM のアラインメント、VCF/BCF の変異表現、BED/GFF/GTF の座標範囲を区別できること。
• 解析前に sample sheet、library type、platform、read length、reference accession / version、checksum、tool version、container / conda environment を記録すること。
• WGS、WES、targeted panel、amplicon、viral sequencing、short-read、long-read では、期待する variant type、coverage、失敗しやすい点が異なること。
• 本章の例は研究・教育用途の最小構成であり、臨床報告、診断、治療方針、感染症サーベイランス上の公式判断にはそのまま使わないこと。

基本概念

ゲノム解析では、read、reference、alignment、variant、annotation を別々の成果物として管理する。FASTQ は観測された配列と quality score、BAM / CRAM は参照配列への配置、VCF / BCF は変異候補と genotype / allele frequency、annotation table は外部 DB に基づく意味付けを表す。

重要なのは、これらを単一の「正解ファイル」として扱わず、参照配列、ツール、パラメータ、DB version、実行ログ、確認日とセットで読むことである。variant annotation は病的意義を自動確定するものではなく、研究・教育用途の候補整理として扱う。

入力データと出力データ

標準的なワークフロー

1. 解析目的、サンプル、library type、参照配列、利用条件を manifest に固定する。
2. FASTQ の品質、adapter、read length、contamination、sample sheet を確認する。
3. 参照配列の accession / version / checksum を記録し、index を作成する。
4. mapping、sort、index、duplicate / coverage 指標、variant calling を実行し、ログを保存する。
5. VCF normalization、filter、annotation を行い、DB version と確認日を残す。
6. QC summary、variant summary、限界、再実行手順をレポートにまとめる。
7. annotation 結果を臨床判断として断定せず、必要に応じて専門家レビューや追加検証へ接続する。

実務上の落とし穴

4.1 配列解析パイプライン

コマンドライン例（最小ワークフロー）: この手順は、Ubuntu 22.04/24.04 相当のLinux環境で、BWA / SAMtools / bcftools を導入済みであることを前提にした環境依存例です。実行時は bwa --version、samtools --version、bcftools --version をログに残し、再現できない場合は公式手順またはコンテナ / conda 環境で同じ入出力を確認してください。 ⚠️ 環境依存（BWA/SAMtools/bcftools 手順例）

# 1) 参照インデックス作成（例: BWA）
bwa index ref.fa

# 2) マッピング → SAM
bwa mem -t 8 ref.fa sample_R1.fq.gz sample_R2.fq.gz > sample.sam

# 3) SAM→BAM & ソート
samtools view -bS sample.sam | samtools sort -o sample.sorted.bam
samtools index sample.sorted.bam

# 4) 変異検出（例: bcftools）
bcftools mpileup -Ou -f ref.fa sample.sorted.bam | \
  bcftools call -mv -Ob -o sample.bcf
bcftools index sample.bcf
bcftools view -Oz -o sample.vcf.gz sample.bcf
bcftools index -t sample.vcf.gz

# 5) 簡易QC/統計
samtools flagstat sample.sorted.bam > sample.flagstat.txt
bcftools stats sample.bcf > sample.bcftools.stats.txt

この例では、変異ファイルのインデックスとして sample.bcf.csi と、必要に応じて sample.vcf.gz.tbi を期待成果物に含めます。plain VCF だけを出力する運用にする場合は、最小入出力のインデックス要件も合わせて緩めます。

QCの基本（FastQC/MultiQC）: FastQC / MultiQC はインストール方法やレポート項目が版により変わるため、実行ログに fastqc --version と multiqc --version を残します。手元で実行できない場合は、公式ドキュメントやサンプルレポートを使い、どのQC指標を読むかを代替確認してください。 ⚠️ 環境依存（FastQC/MultiQC 手順例）

# 入力FASTQの品質評価（並列実行例）
mkdir -p qc/fastqc
fastqc -t 8 -o qc/fastqc data/*.fq.gz

# レポート統合
multiqc -o qc qc/fastqc

前処理フェーズ: 🔁 疑似コード（前処理順序の説明例）

def preprocess_sequence(raw_sequence):
    # 品質管理
    quality_filtered = quality_control(raw_sequence)

    # アダプター除去
    trimmed = adapter_trimming(quality_filtered)

    # 重複除去
    deduplicated = remove_duplicates(trimmed)

    return deduplicated

実装例（要環境調整）: このコードは構造の例です。依存ライブラリ（BioPython / pysam / NumPy 等）や入出力、例外設計は環境に合わせて調整してください（CIの実行対象ではありません）。実行する場合は Python と各ライブラリの版を固定し、最小FASTQ/BAMを使った単体確認または公式チュートリアルで代替確認します。 ⚠️ 環境依存（BioPython / pysam 依存の構造例）

import logging
import os
from typing import List, Dict, Optional, Tuple
from dataclasses import dataclass
from Bio import SeqIO
import pysam
import numpy as np

# ロギング設定
logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s - %(name)s - %(levelname)s - %(message)s'
)
logger = logging.getLogger(__name__)

@dataclass
class QualityMetrics:
    """品質管理メトリクス"""
    total_reads: int
    passed_reads: int
    mean_quality: float
    q20_percentage: float
    q30_percentage: float
    gc_content: float

class SequenceAnalysisPipeline:
    """エラーハンドリングを含む配列解析パイプライン"""

    def __init__(self, min_quality: int = 20, min_length: int = 50):
        self.min_quality = min_quality
        self.min_length = min_length

    def quality_control(self, fastq_file: str) -> Tuple[str, QualityMetrics]:
        """
        品質管理とフィルタリング

        Args:
            fastq_file: 入力FASTQファイル

        Returns:
            tuple: (フィルタ済みファイル名, 品質メトリクス)
        """
        try:
            output_file = fastq_file.replace('.fastq', '_filtered.fastq')

            total_reads = 0
            passed_reads = 0
            quality_scores = []
            gc_counts = []

            with open(output_file, 'w') as out_handle:
                for record in SeqIO.parse(fastq_file, 'fastq'):
                    total_reads += 1

                    # 品質スコアの計算
                    quals = record.letter_annotations['phred_quality']
                    mean_qual = np.mean(quals)

                    # 長さチェック
                    if len(record.seq) < self.min_length:
                        logger.debug(f"Read {record.id} too short: {len(record.seq)}")
                        continue

                    # 品質チェック
                    if mean_qual < self.min_quality:
                        logger.debug(f"Read {record.id} low quality: {mean_qual:.2f}")
                        continue

                    # GC含量計算
                    gc_count = record.seq.count('G') + record.seq.count('C')
                    gc_content = gc_count / len(record.seq)

                    # フィルタを通過
                    passed_reads += 1
                    quality_scores.extend(quals)
                    gc_counts.append(gc_content)

                    # 出力
                    SeqIO.write(record, out_handle, 'fastq')

            # メトリクス計算
            if passed_reads == 0:
                raise ValueError("No reads passed quality filters")

            quality_array = np.array(quality_scores)
            metrics = QualityMetrics(
                total_reads=total_reads,
                passed_reads=passed_reads,
                mean_quality=np.mean(quality_array),
                q20_percentage=np.sum(quality_array >= 20) / len(quality_array) * 100,
                q30_percentage=np.sum(quality_array >= 30) / len(quality_array) * 100,
                gc_content=np.mean(gc_counts)
            )

            logger.info(f"Quality control completed: {passed_reads}/{total_reads} reads passed")
            return output_file, metrics

        except FileNotFoundError:
            logger.error(f"Input file not found: {fastq_file}")
            raise
        except Exception as e:
            logger.error(f"Quality control failed: {e}")
            raise

マッピングフェーズ:

• BWA: Burrows-Wheeler Alignerによる高速マッピング
• SAMtools: マッピング結果の操作・統計
• GATK: 変異検出パイプライン

後処理フェーズ:

• 変異アノテーション: 機能予測・既存DBとの照合（臨床解釈では専門家レビューが必要）
• 品質管理: 偽陽性率の制御
• 統計解析: 有意性検定・多重検定補正

4.1.1 参照ゲノムの選択（GRCh38 / T2T-CHM13 / HPRC pangenome）

参照ゲノムは「新しいものへ単純に置き換える」対象ではなく、解析目的、既存DBとの互換性、座標系、annotation、臨床報告の前提に応じて選びます。GRCh38 は現在も annotation、既存ツール、臨床・研究DB連携で中心的に使われます。一方、T2T-CHM13 は従来の未解読領域や反復領域を扱う研究に有用で、HPRC pangenome / graph genome は多数の haplotype を含む参照表現として、構造多型や reference bias の評価に向きます。

参照変更時は、次の項目を manifest または reference_version テーブルへ必ず残します。

Graph reference は、多数の haplotype や構造差分を含む参照表現であり、read mapping と variant genotyping / discovery に使われます。GFA は segment / link / path / walk で graph 構造を表すテキスト形式、GAF は read や contig が graph 上のどの path に対応したかを表す alignment 形式です。初学者向けには HLA や KIR のような複雑な実データではなく、小さな人工配列で線形参照と graph path の違いを確認します。

🔁 疑似コード（実行不可: toy graph で線形参照と graph path を比較する例）

linear_ref: A-C-G-T
sample_a:   A-C-G-T
sample_b:   A-G-G-T

graph_segments: A, C, G1, G2, T
graph_paths:
  ref_path: A -> C  -> G1 -> T
  alt_path: A -> G2 -> G1 -> T

目的: graph path が allele を表現できることを理解する。
注意: 実ヒトデータの臨床解釈やHLA解析をこの toy graph から一般化しない。

解析パイプラインの3系統（概要）:

• short-read: QC → 線形参照へのマッピング → 変異コーリング → 注釈付与
• long-read: QC → ロングリードマッピング/アセンブリ → 構造変異検出 → ハプロタイプ解析
• pangenome/graph: graph build / release確認 → graph alignment（GAF等） → variant genotyping / discovery → VCF等への変換 → annotation対応assemblyの確認

研究解析と臨床実装の境界: 標準化されたパイプラインは、FASTQ / BAM / VCF と QC 指標を再現可能に整理し、研究・教育用の変異候補や品質サマリを専門家レビューへ渡すための基盤になる。ただし、本章のコマンド例や概念例だけで診断、治療方針、投薬、緊急症例の判断を行うものではない。臨床実装では、検査法の妥当性確認、施設SOP、承認済みパイプライン、専門職レビュー、監査ログ、規制・倫理手順が別途必要である。

4.2 変異検出と機能解析

体細胞変異検出: 🔁 疑似コード（変異検出処理順序の説明例）

class SomaticVariantCaller:
    """腫瘍-正常ペアからの体細胞変異検出"""

    def __init__(self, tumor_bam, normal_bam, reference):
        self.tumor_bam = tumor_bam
        self.normal_bam = normal_bam
        self.reference = reference

    def call_variants(self, min_vaf=0.05, min_depth=20):
        """
        体細胞変異の検出

        Args:
            min_vaf: 最小変異アレル頻度
            min_depth: 最小リード深度
        """
        # MuTect2スタイルの変異検出
        variants = []

        # 腫瘍特異的変異の統計的検定
        # Fisher's exact testによる有意性評価

        return variants

構造変異（SV）検出:

• 欠失・挿入・逆位・転座の検出
• Split-read/Paired-end情報の統合
• CNV（コピー数変異）解析

深層学習ベースの変異コーリング:

• 体細胞変異の文脈では、DeepSomaticやNeuSomaticなどの深層学習ベースcallerが提案されている。
• 生殖細胞系列では、DeepVariant等が画像化したリード情報を用いて変異を推定する。
• これらのツールは参照ゲノムのバージョン、学習データの前提、カバレッジ条件に依存するため、適用範囲の確認が必要。

変異の機能的影響予測:

• SIFT/PolyPhen-2などによる機能影響予測（単独で病原性を確定しない）
• VEP（Variant Effect Predictor）
• スプライシング影響予測

4.3 ロングリードシークエンシング

PacBio/Oxford Nanopore技術:

• リード長: 10kb-100kb+
• エラー率・精度: chemistry、basecaller、read type（HiFi / ONT など）、model に依存するため、固定値として扱わない
• リアルタイムシークエンシング

4.3.1 Long-read WGS の最小ワークフロー

Long-read WGS は、short-read WGS を置き換える単一の上位互換ではなく、目的に応じて使い分ける解析系です。SNV/short indel を高スループットに扱う既存解析では short-read が成熟しています。一方、PacBio HiFi や Oxford Nanopore は、SV、repeat expansion、挿入配列、phasing、複雑領域、T2T / pangenome 参照で変わる領域を調べるときに有効です。caller や basecaller、chemistry、学習 model、参照ゲノムが結果に影響するため、解析単位ごとに固定して記録します。

Long-read WGS の教育用最小ワークフローは次のように切り分けます。

この表は臨床実装の標準手順ではありません。研究・教育用途でも、DeepVariant / Clair3 などの small variant caller は入力プラットフォーム・model・reference に依存し、Sniffles2 / pbsv / cuteSV などの SV caller は SV の種類、repeat、coverage、filter 条件で結果が変わります。したがって、caller_model、reference_version、annotation_release、database_release、container_digest、workflow_commit、retrieved_at を第2章・第10章の provenance 記録へ渡します。

ハイブリッド解析（short+long）:

• short-readの精度とlong-readの連続性を組み合わせ、構造変異や反復領域の解析精度を高める。
• 解析設計としては、short-readを基準としつつ、long-readで難所を補完する形が一般的。

アセンブリアルゴリズム: 🔁 疑似コード（アセンブリ処理順序の説明例）

class LongReadAssembler:
    """ロングリードアセンブリ"""

    def overlap_layout_consensus(self, reads):
        # 1. オーバーラップグラフ構築
        overlap_graph = self.build_overlap_graph(reads)

        # 2. レイアウト（グラフ簡約化）
        layout = self.simplify_graph(overlap_graph)

        # 3. コンセンサス配列生成
        consensus = self.generate_consensus(layout)

        return consensus

ハプロタイプ分離アセンブリ:

• 父母由来の染色体を分離
• HiFiリードによる高精度化
• 構造変異の正確な検出

小さな実例

SRR11140744 と MN908947.3 を使う最小例では、FASTQ 取得、FastQC / MultiQC、mapping、BAM index、variant calling、VCF / BCF summary までを小さく通す。目的は高精度な解析結果を得ることではなく、どの段階でどのファイルとログが発生するかを確認することである。

結果の読み方

mapping rate、coverage、variant count は、サンプル品質、library type、参照配列、target region に依存する。単一の閾値で成功・失敗を決めず、期待する解析目的と比較対象を明記して読む。

VCF / BCF は候補変異の表現であり、annotation は外部 DB の時点・submitter・population・evidence level に依存する。ClinVar、gnomAD、VEP / SnpEff などの結果は、研究・教育用途の候補整理として扱い、診断・治療判断へ直接接続しない。

限界と注意点

short-read SNV pipeline だけでは、SV、repeat expansion、複雑な indel、methylation、pangenome / T2T で変わる領域を十分に扱えない場合がある。long-read、orthogonal validation、liftover / remap、pangenome graph などが必要な条件をレポートに明記する。

🎯 認定試験ポイント

重要概念チェックリスト

ゲノム解析基礎 ⭐⭐⭐

• ゲノムサイズ・遺伝子数・非コード領域の割合をヒトで把握している
• SNP・Indel・CNV・構造変異の定義と検出手法を理解している
• 遺伝率とは何か、その推定方法を説明できる
• GWAS（ゲノムワイド関連解析）の基本原理と限界を理解している

配列解析技術 ⭐⭐⭐

• Sanger法と次世代シークエンシング技術の特徴比較ができる
• リードマッピング・アセンブリの基本戦略を理解している
• ペアエンドリード・ロングリードの利点と使い分けを把握している
• カバレッジ・精度・N50等の品質指標の意味を説明できる

変異解析 ⭐⭐

• SNP・Indelコーリングのアルゴリズムと品質フィルタリング
• アレル頻度・ハーディー・ワインベルグ平衡の計算
• 連鎖不平衡（LD）とハプロタイプの概念
• 集団遺伝学パラメータ（FST、Tajima’s D等）の解釈

アセンブリ手法 ⭐⭐

• グリーディアルゴリズムとグラフベースアセンブリの違い
• De Bruijnグラフの構築とパス探索
• アセンブリ品質評価（N50、コンティグサイズ分布等）
• ギャップクローニングとスキャフォールディング

典型的な出題パターン

【遺伝学計算】 🧪 概念例

問題例: 日本人集団でのSNPアレル頻度がA=0.7、a=0.3の場合、
ハーディー・ワインベルグ平衡における各遺伝型頻度は？

解答: AA = 0.7² = 0.49 (49%)
      Aa = 2×0.7×0.3 = 0.42 (42%)
      aa = 0.3² = 0.09 (9%)

【技術選択】 🧪 概念例

問題例: 新規ゲノムのde novoアセンブリを行う際、
ショートリードとロングリードの使い分けを説明せよ。

解答: ショートリード（Illumina）: 高精度、反復配列に課題
      ロングリード（PacBio/ONT）: 反復配列解決、精度向上が課題
      → ハイブリッドアプローチで両者の利点を活用

【統計・品質管理】 🧪 概念例

問題例: GWAS解析で多重検定補正が必要な理由と、
一般的な補正手法を2つ挙げよ。

解答: 理由: 数百万SNPを同時検定するため偽陽性率が急増
      手法: 1) Bonferroni補正（保守的）
           2) FDR補正（Benjamini-Hochberg法）

【データ解釈】 🧪 概念例

問題例: RNA-seqデータで発現変動遺伝子（DEG）を検出した後、
生物学的意味を調べるための解析手法を3つ挙げよ。

解答: 1) Gene Ontology (GO) enrichment解析
      2) KEGGパスウェイ解析
      3) Gene Set Enrichment Analysis (GSEA)

関連する付録・章

• 付録G: JSBi認定試験の理解補助
• 付録H: プログラム・ツール・データベース一覧（ゲノム解析）
• 第3章: データ構造とアルゴリズム（技術基盤）
• 第5章: トランスクリプトーム解析（発現解析）
• 第6章: エピゲノム・マルチオミクス解析（制御解析）

まとめ

• 第4章で固定すべきなのは、FASTQ / BAM / VCF の役割と、各工程で確認する QC 指標です。
• 最初から大規模データを扱うのではなく、SRR11140744 のような小規模公開データで最小パイプラインを再現すると理解しやすくなります。
• 後続章でも、ここで作った成果物とログの残し方が再現性の基盤になります。
• 本章の例は研究・教育用の最小構成であり、臨床報告や診療判断には使いません。

Source notes / 次の一歩

• 2026年版出典監査メモ: GRCh38 / T2T / pangenome / reference bias（GRC、NCBI、HPRC、HPRC Data Use の確認日）
• GRC Human Overview / NCBI T2T-CHM13 / HPRC Data Release 2 / HPRC release timeline / GFA specification / minigraph GAF/rGFA docs（2026-06-05確認）
• minimap2 / pbmm2 / mosdepth / DeepVariant / Clair3 / Sniffles2 / pbsv / cuteSV / WhatsHap（2026-06-05確認）
• 付録F: 第4章の次の一歩
• 付録J: 実在アクセッション一覧
• 付録H: プログラム・ツール・データベース一覧

2026年時点の更新メモ

• GRCh38、T2T-CHM13、HPRC pangenome graph などは座標系と利用条件が異なるため、reference accession、assembly release、checksum、alt / decoy / PAR、liftover / remap の有無、graph build ID、GFA / GAF / VCF version、annotation release を記録する。
• long-read WGS は SV、repeat、phasing、複雑領域の解析で有用だが、basecaller、chemistry、caller model、reference、coverage、filter に依存する。small variant caller、SV caller、phasing tool、annotation DB の version と container digest を manifest に残す。
• variant annotation は DB version、submitter、population、evidence level、確認日に依存する。臨床的意義を本文の概念例だけで確定しない。
• 更新根拠は 2026年版出典監査メモの GRCh38 / T2T / pangenome / reference bias 項目と、付録Fの第4章項目へ集約する（確認日: 2026-05-12）。

次章への橋渡し

ゲノム解析で「配列から変異を取り出す」流れを固定したら、次は 「その変化が発現量や機能にどう現れるか」を追う必要があります。 第5章では、同じ再現性の考え方を保ちながら RNA-seq と発現解析へ進みます。

最小入出力（期待成果物/期待ログ）

• 入力: 参照ゲノム（FASTA）とリード（FASTQ）（小規模データで可）
• 出力（期待成果物）: ソート済みBAM＋インデックス、変異（BCF＋インデックス、必要に応じて圧縮VCF＋タビックスインデックス）、QCレポート（FastQC/MultiQC）、研究・教育用の品質サマリ
• 期待ログ（例）: accession.version、reference_version、checksum、マッピング率・read数・変異数などのサマリ（例: samtools flagstat, bcftools stats）

演習

1. SRR11140744 と MN908947.3 を用い、QC→アライメント→BAM 整形→変異検出までの最小パイプラインを実行し、成果物（BAM / VCF / QC レポート）を作成せよ。
2. トリミング閾値、最小品質、最小深度を 1 つずつ変え、マッピング率・変異数・主要 QC 指標の変化を比較せよ。
3. 使用したコマンド、パラメータ、ツール / 参照 DB のバージョンを README 形式で残し、第三者が再実行できるよう整理せよ。
4. 小さな人工配列で、線形参照にしか存在しない allele と graph path として表現できる allele を比較し、GFA / GAF / VCF の役割差と annotation へ渡す前提を説明せよ。

具体課題例

• samtools flagstat と bcftools stats を使い、最小限の品質サマリを 1 ページで説明する。
• マッピングツールまたは変異検出ツールを 1 つ差し替え、結果差分と採用判断を文章でまとめる。
• 第10章で再利用できるよう、出力ファイル名・サンプル ID・参照配列 ID の命名規則を決める。

Tiny end-to-end ラボ

examples/labs/genome/ には、合成FASTAを入力としてQCと成果物検証までを通すネットワーク不要の最小ラボがあります。実行入口は make lab-smoke です。これは本格的なalignment、variant calling、臨床解釈を代替しません。