bioinformatics-fundamentals

バイオインフォマティクスの基礎

ゲノミクスおよびバイオインフォマティクスのワークフローに関する基礎知識です。ファイル形式、シーケンシング技術、一般的なデータ処理パターンを理解するために必要な知識を提供します。

このスキルを使用する場面

• シーケンシングデータ（PacBio HiFi、Hi-C、Illumina）を扱う場合
• SAM/BAM のアライメントやフィルタリングに関する問題をデバッグする場合
• ゲノムアセンブリのキュレーションで AGP ファイルを処理する場合
• AGP の座標系と unloc の割り当てを検証する場合
• ペアエンドデータとシングルエンドデータの違いを理解する場合
• 品質指標（MAPQ、PHRED スコア）を解釈する場合
• 空の出力や破損したリードペアのトラブルシューティングを行う場合
• GenomeArk の QC データ（GenomeScope、BUSCO、Merqury）にアクセスする場合
• 核型データのキュレーションや染色体数の解析を行う場合
• 一般的なバイオインフォマティクスデータ解析を行う場合

SAM/BAM 形式の基本

SAM フラグ（ビット単位）

フラグは加算式です。つまり、1 つのリードに複数のフラグを同時に設定できます。

一般的なフラグ：

• 0x0001 (1)：リードがシーケンシングでペアになっている
• 0x0002 (2)：各セグメントが正しくアライメントされている（適正なペア）
• 0x0004 (4)：リードがマッピングされていない
• 0x0008 (8)：メイトがマッピングされていない
• 0x0010 (16)：リードが逆鎖にマッピングされている
• 0x0020 (32)：メイトが逆鎖にマッピングされている
• 0x0040 (64)：ペアの 1 番目（R1／フォワード）
• 0x0080 (128)：ペアの 2 番目（R2／リバース）
• 0x0100 (256)：セカンダリアライメント
• 0x0400 (1024)：PCR または光学的重複
• 0x0800 (2048)：補足アライメント

フラグの組み合わせ：

• 適正にペアリングされた R1：99 (0x63 = 1 + 2 + 32 + 64)
• 適正にペアリングされた R2：147 (0x93 = 1 + 2 + 16 + 128)
• マッピングされていないリード：4
• メイトがマッピングされていない：8

フラグの一覧、CIGAR 操作、オプションタグ、SAM の必須フィールドについては、reference.md を参照してください。

適正ペアフラグ（0x0002）

「適正なペア」の意味：

• R1 と R2 の両方がマッピングされている
• マッピングの向きが正しい（通常は R1 がフォワード、R2 がリバース）
• インサートサイズがライブラリとして妥当である
• ペアがアライナーの想定に合致している

重要： 適正なペアの判定基準は、アライナーによって異なります。

MAPQ（マッピング品質）

計算式： MAPQ = -10 * log10(P(mapping is wrong))

一般的な閾値：

• MAPQ >= 60：高信頼度（誤りの確率 < 0.0001%）
• MAPQ >= 30：良好な品質（誤りの確率 < 0.1%）
• MAPQ >= 20：許容可能（誤りの確率 < 1%）
• MAPQ >= 10：低信頼度（誤りの確率 < 10%）
• MAPQ = 0：マルチマッパー、または未マッピング

注： MAPQ=0 は、未マッピングの場合と、同等に良好なマッピング先が複数ある場合のどちらも意味することがあります。

CIGAR 文字列

リードとリファレンスのアライメントを表します。

• M：一致または不一致（アライメント上の一致）
• I：リファレンスに対するリードの挿入
• D：リファレンスに対するリードの欠失
• S：ソフトクリッピング（リード内でアライメントされていない塩基）
• H：ハードクリッピング（リード配列に含まれない塩基）
• N：スキップ領域（RNA-seq のスプライシングに使用）

例： 100M = 100 bp の完全一致
例： 50M5I45M = 50 bp の一致、5 bp の挿入、45 bp の一致

シーケンシング技術

PacBio HiFi（高忠実度）

特徴：

• ロングリード：通常 10～25 kb
• 高精度：>99.9%（Q20 以上）
• サーキュラーコンセンサスシーケンシング（CCS）
• シングルエンドデータ（ただし環状分子由来）
• de novo アセンブリに最適

最適なマッパー：

• minimap2 のプリセット：map-pb、map-hifi
• BWA-MEM2 も使用できますが、ショートリード向けに最適化されています

一般的な用途：

• de novo ゲノムアセンブリ
• 構造変異の検出
• アイソフォームシーケンシング（Iso-Seq）
• ハプロタイプのフェージング

Hi-C（クロマチンコンフォメーションキャプチャ）

特徴：

• ペアエンドのショートリード（通常 100～150 bp）
• リードペアがクロマチン相互作用を捉えます
• R1 と R2 は別々のスキャフォールド／染色体にマッピングされることがよくあります
• ペアの適切な処理が必要です
• スキャフォールディングと 3D ゲノム構造の解析に使用されます

推奨マッパー：

• BWA-MEM2（ペアエンドモード）
• BWA-MEM（ペアエンドモード）

重要な概念： Hi-C のリードペアは、意図的に離れた遺伝子座にマッピングされます。領域フィルタリングによって、ペアが簡単に分断されることがあります！

一般的な用途：

• ゲノムスキャフォールディング（コンティグの連結）
• 3D クロマチン構造解析
• ハプロタイプフェージング
• アセンブリ品質評価

Illumina ショートリード

特徴：

• ショートリード：50～300 bp
• ペアエンドまたはシングルエンド
• 高スループット
• 確立された品質スコア

推奨マッパー：

• BWA-MEM2、BWA-MEM（汎用）
• Bowtie2（高速、ローカルアライメント）
• STAR（RNA-seq のスプライスアライメント）

各テクノロジーの詳細な技術仕様、エラープロファイル、品質指標については、reference.md を参照してください。

一般的なツールとその動作

samtools view

目的： SAM/BAM ファイルのフィルタリング、変換、表示

主なフラグ：

• -b：BAM 形式で出力
• -h：ヘッダーを含める
• -f INT：フラグを必須にする（これらのフラグを持つリードを保持）
• -F INT：フラグでフィルタリングする（これらのフラグを持つリードを除外）
• -q INT：MAPQ の最小しきい値
• -L FILE：BED ファイル内の領域と重なるリードを保持

重要な動作：

• -L（領域フィルタリング）は、ペア単位ではなく個々のリードをチェックします
• 片方のリードが別の領域にマッピングされている場合、リードペアが分断されることがあります
• フラグフィルター（-f、-F）は、領域フィルター（-L）の前に適用されます

例 — 領域内の適切なペア（正しい順序）：

samtools view -b -f 2 -L regions.bed input.bam > proper_pairs_in_regions.bam

bamtools filter

目的： 複雑な条件による高度なフィルタリング

一般的なフィルター：

• isPaired: true - リードがペアエンドシーケンシング由来である
• isProperPair: true - リードが適切なペアの一部である
• isMapped: true - リードがマッピングされている
• mapQuality: >=30 - マッピング品質のしきい値

samtools との重要な違い：

• isProperPair は samtools の -f 2 より厳密です
• ペアの妥当性をより詳細にチェックします

samtools fastx

目的： SAM/BAM を FASTQ/FASTA に変換

重要： R1/R2 ファイルが一致するよう、適切なフィルターを使用してください！

すべてのオプションと例を含む、ツールコマンドの完全なリファレンスについては、reference.md を参照してください。

一般的なパターンとベストプラクティス

パターン 1：領域によるペアエンドデータのフィルタリング

誤った方法（ペアが分断されます）：

# Region filter first -> breaks pairs when mates are in different regions
samtools view -b -L regions.bed input.bam | bamtools filter -isPaired -isProperPair
# Result: Empty output (all pairs broken)

正しい方法（ペアが保持されます）：

# Proper pair filter FIRST, then region filter
samtools view -b -f 2 -L regions.bed input.bam > output.bam

パターン 2：フィルタリング済み BAM からの FASTQ 抽出

ペアエンドの場合：

samtools fastx -1 R1.fq.gz -2 R2.fq.gz \
  --i1-flags 2 \  # Require proper pair
  input.bam

シングルエンドの場合：

samtools fastx -0 output.fq.gz input.bam

パターン 3：品質フィルタリング

保守的（高品質）：

samtools view -b -q 30 -f 2 -F 256 -F 2048 input.bam
# MAPQ >= 30, proper pairs, no secondary/supplementary

寛容（低カバレッジデータ向け）：

samtools view -b -q 10 -F 4 input.bam
# MAPQ >= 10, mapped reads

よくある問題のまとめ

問題 1：領域フィルタリング後に出力が空になる（Hi-C データ）

領域フィルター（samtools view -L）によってリードペアが分断されます。片方のリードは領域内、もう片方は領域外にあります。適切なペアのフラグが失われます。領域フィルタリングの前に、適切なペアのフィルターを適用してください：

samtools view -b -f 2 -L regions.bed input.bam > output.bam

問題 2：R1 ファイルと R2 ファイルのリード数が異なる

不適切なフィルタリングによって、一部のペアが分断されています。抽出時に適切なペアを必須にしてください：

samtools fastx -1 R1.fq -2 R2.fq --i1-flags 2 input.bam

問題 3：Hi-C データのマッピング率が低い

キメラリードがあるため、Hi-C では正常です。Hi-C 専用のパイプライン（HiC-Pro、Juicer）を使用してください。MAPQ によるフィルタリングを厳しくしすぎないでください。

問題 4：マッピング後に適切なペアが失われる

インサートサイズの分布、リファレンスの不一致、または誤った配向フラグを確認してください。

samtools stats input.bam | grep "insert size"
samtools flagstat input.bam

AGP の処理に関する問題、HiFi 固有の問題、診断コマンドなど、詳細な解決策を含む包括的なトラブルシューティングについては、common-issues.md を参照してください。

品質指標

N50 と関連指標

N50： アセンブリ全体の 50％がその長さ以上のコンティグに含まれるとき、該当する最短コンティグの長さ

関連指標：

• L50： N50 に達するまでに必要なコンティグ数
• N90： N50 より厳しい指標（90％のカバレッジ）
• NG50： ゲノムサイズを基準とした N50（比較に適しています）

カバレッジと深度

カバレッジ： 少なくとも 1 本のリードでカバーされているリファレンス塩基の割合 深度： 各塩基をカバーするリード数の平均

推奨深度：

• ゲノムアセンブリ（HiFi）：30〜50x
• バリアントコーリング：最低 30x
• RNA-seq：2,000 万〜4,000 万リード
• Hi-C スキャフォールディング：ゲノムカバレッジ 50〜100x

カバレッジの計算方法、BUSCO の解釈、QV スコア、アセンブリ品質指標の詳細については、reference.md を参照してください。

ファイル形式のクイックリファレンス

FASTA

>sequence_id description
ATCGATCGATCG

• ヘッダー行は > で始まります
• 品質スコアはありません

FASTQ

@read_id
ATCGATCGATCG
+
IIIIIIIIIIII

• 1 リードにつき 4 行です
• 品質スコアがあります（通常は Phred+33 エンコーディング）

BED

chr1    1000    2000    feature_name    score    +

• 0 ベース座標
• 半開区間 [start, end)

AGP

chr1    1    5000    1    W    contig_1    1    5000    +
chr1    5001 5100    2    U    100    scaffold    yes    proximity_ligation

• タブ区切りのゲノムアセンブリ形式
• 1 ベースの閉区間 [start, end]
• オブジェクトとコンポーネントの長さは一致する必要があります：obj_end - obj_beg + 1 == comp_end - comp_beg + 1

AGP の完全な仕様、座標系、検証ルール、処理パターンについては、reference.md を参照してください。 AGP 座標のデバッグと unloc の処理に関する問題については、common-issues.md を参照してください。

座標系

1 ベース（SAM、VCF、GFF、AGP）： 最初の塩基は位置 1 です。区間 [2,5] には位置 2、3、4、5 が含まれます。 0 ベース（BED、BAM バイナリ）： 最初の塩基は位置 0 です。区間 [2,5) には位置 2、3、4 が含まれます（5 は含まれません）。

変換： BED_start = SAM_start - 1；BED_end = SAM_end

ベストプラクティス

一般

1. データ型を必ず確認する： ペアエンドかシングルエンドかによって、フィルタリング方法が決まります
2. シーケンシング技術を理解する： Hi-C は HiFi とは異なる挙動をします
3. 適切な順序でフィルタリングする： 領域でのフィルタリングより先に適正ペアを処理します
4. 出力を検証する： ファイルサイズ、リード数、flagstat を確認します
5. 適切な MAPQ しきい値を使用する： 厳しすぎるとデータが失われ、緩すぎるとノイズが増えます

Hi-C データの場合

1. 離れたリードペアを想定する： 異なるスキャフォールドにマッピングされても驚かないでください
2. 適正ペアを保持する： 下流のスキャフォールディングに不可欠です
3. ペアを考慮したツールを使用する： 標準的なフィルターではペアが壊れる場合があります
4. MAPQ で過剰にフィルタリングしない： Hi-C は DNA-seq より MAPQ が低いことがよくあります

HiFi データの場合

1. シングルエンド処理： ペアに関する懸念はありません
2. 高品質を想定する： 厳しいフィルターを使用できます
3. 適切なプリセットを使用する： minimap2 map-hifi または map-pb
4. リード長の分布を考慮する： HiFi リードの長さにはばらつきがあります

ツールのテストの場合

1. 自己完結型のデータセットを作成する： 選択した領域に両方のメイトを含めます
2. 適正ペアを維持する： 現実的なテストに不可欠です
3. 代表性のあるデータを使用する： ランダムではなく、割合を保ってサブサンプリングします
4. ファイルサイズを確認する： 小さすぎる場合は、過剰にフィルタリングされています

関連スキル

• vgp-pipeline - VGP ワークフローでは Hi-C と HiFi のデータを処理します
• galaxy-tool-wrapping - Galaxy ツールは SAM/BAM やシーケンシングデータ形式を扱います
• galaxy-workflow-development - ワークフローでシーケンシングデータを処理します

補足ドキュメント

• reference.md： 詳細な形式仕様（SAM/BAM の完全なリファレンス、CIGAR 操作、AGP 形式、FASTQ エンコーディング、ツールコマンドのリファレンス、シーケンシング技術の仕様、アセンブリ品質指標、カバレッジ計算、座標系）
• common-issues.md： 包括的なトラブルシューティングガイド（出力が空になる問題、ペアエンドの問題、品質やマッピングの問題、形式変換、Hi-C と HiFi 固有の問題、AGP 処理エラー、診断コマンド）
• genomeark-data-access.md： GenomeArk の AWS S3 データアクセスパターン（ディレクトリ構造の変遷、GenomeScope/BUSCO/Merqury の QC データの場所、取得方法、パスの正規化、アセンブリ日の抽出）
• genomic-analysis-patterns.md： 分野固有の解析パターン（核型データのキュレーション、半数体と二倍体の染色体数、系統樹における種のマッピング、BED/テロメア解析、NCBI データ統合の戦略）

バージョン履歴

• v1.2.0： 保守性を高めるため、SKILL.md と補足ファイルに分割。GenomeArk へのアクセス、核型／染色体解析、系統マッピング、AGP の詳細、テロメア／NCBI のパターンを補足ファイルに移動
• v1.1.1： テロメア解析のための BED ファイル処理パターンと、NCBI データ統合の戦略を追加
• v1.1.0： 座標検証、unloc 処理、一般的なエラーパターンを含む、包括的な AGP 形式のドキュメントを追加
• v1.0.0： SAM/BAM、Hi-C、HiFi、一般的なフィルタリングパターンを含む初回リリース