はじめに
この記事は、HTSlibのリリースノートを Claude や Gemini に入力して得られた複数の出力を基に、人間の kojix2 が全文を目を通し、再構成することで作成されています。生成AIを利用していますが、読み応えのある内容になっていると思います。
htslib は SAM、BAM、CRAM、VCFといったNGSデータの基本フォーマットを扱うためのライブラリです。
このライブラリのリリースノートの履歴を眺めていると、「今、研究の現場で何が新しいか」「どういう解析が普及しつつあるのか」という動向が見えてくることがあります。
この記事では、2022年から2024年にかけてのhtslibのリリース内容を基に、そこから見えてくるゲノム科学の潮流について考察します。
ロングリードと修飾塩基
塩基修飾
-
v1.18:
MNタグの追加（ハードクリップによるデータ破損の検証用） - v1.18: 塩基修飾APIの改善（未チェック塩基と未修飾塩基を区別できるように）
- v1.20: タグの妥当性チェック強化
- v1.23: 空のMMタグで落ちるバグ修正
最近のhtslibは、SAM/BAMファイルの塩基修飾タグ（MM/ML）の周りを継続的に改善しています。
- MMタグ：Nanoporeなどで推定した塩基修飾（例：m⁶A, 5mC）の位置と種類を記録するタグ。
- MLタグ：MMタグで示された各修飾についての修飾である確率（信頼度；0–255）を数値で表したタグ。
これまで、DNAのメチル化を調べるには、薬品でDNAを化学変換し、未メチル化シトシンをウラシルに変えてから読むバイサルファイト処理が必要でした。これはDNAに損傷が発生し手間もかかります。NanoporeやPacBioなどのロングリードシーケンサーは、DNAを化学変換せずに、配列（A/T/G/C）と修飾（メチル化など）を直接読み取ります。これにより、実験データとしてではなく、リード情報の一部として修飾塩基が出力されます。
FASTQには修飾塩基情報を掲載する場所がないので、生の信号データである pod5 形式や、uBAM/modBAM（未アライメントBAM）といった形式も見かけるようになりました。
htslibで、メチル化など修飾塩基の情報をBAMファイルから扱えるようにするコード。
// 状態の確保と初期化
hts_base_mod_state *m = hts_base_mod_state_alloc();
bam_parse_basemod2(b, m, HTS_MOD_REPORT_UNCHECKED);
// 修飾タイプ一覧の取得
int ntypes;
int *types = bam_mods_recorded(m, &ntypes);
// 各塩基位置での修飾情報を取得
hts_base_mod mods[5];
for (int i = 0; i < b->core.l_qseq; i++) {
int n = bam_mods_at_next_pos(b, m, mods, 5);
// 修飾情報を処理
}
// 後片付け
hts_base_mod_state_free(m);
構造変異の表現力の向上
ロングリードのもう一つの強みは、構造変異の検出です。
従来のショートリード（短い断片を読む技術）では、染色体が大きく入れ替わるような構造変異を正確に捉えるのは苦手でした。 ロングリードで数万塩基を一度に読めるようになり、複雑な転座や逆位の観察が可能になってきました。また、ヒトゲノムの完全なリファレンス配列（T2T）も登場しました。このようなケースでVCF形式ではBreakend記法が使われます。htslibでもBreakend表記を正しく検出できるようになりました。
-
v1.22:
T > A[chr15:12345[のような複雑なbreakend表記をサポート
// VCFレコード読み込み後にbreakendを検出
bcf_set_variant_types(rec);
if (bcf_has_variant_type(rec, 1, VCF_BND)) {
// breakend特有の処理を実行
}
シングルセル解析と分子バーコード
UMI（分子バーコード）の公式サポート
- v1.23: FASTQ入出力でのUMIサポート
UMI（Unique Molecular Identifier）は、PCRバイアスを補正して分子数を正確にカウントするための技術です。
1つの細胞から取れるRNAは微量なので、PCRで何倍にも増幅します。すると、元々たくさんあった分子なのか、PCRでたまたま増えた分子なのか区別がつかなくなります。 そこで、増幅前にランダムなバーコード（UMI）を付けておき、「同じバーコードを持つ配列は、コピーとみなして1つにまとめる」ことで、正確な分子数を数えます。
htslib に追加されたのは、リード名からUMIを抽出して、それをSAM/BAMのRXタグに保存したり、逆にRXタグからリード名にUMIを追加するような単純な機能です。しかし、これがhtslib本体に入ったということは、シングルセル解析が以前より一般的なアプローチになっていることを意味します。
// この1行で以下の処理が自動的に行われる
// - 読み込み時: `fastq_parse1()` がリード名からUMIを抽出し、RXタグに保存
// - 書き出し時: `fastq_format1()` がRXタグからUMIを読み取り、リード名に追加
hts_set_opt(in, FASTQ_OPT_UMI, "RX");
集団遺伝学と大規模コホート研究
数万人〜数十万人規模のコホート研究（UK Biobankなど）では、ストレージコストが高額になります。また私たちの日々の解析でも、ゲノムデータは容易にディスクを圧迫します。
CRAMファイル形式の改善
CRAMは、参照配列（リファレンス）との差分だけを記録することでBAMのサイズを縮小するフォーマットです。CRAM 3.1では、GZIP, BZIP2, LZMA, RANS などの従来の圧縮アルゴリズムに加えて、RANSPR、ARITH、FQZ、TOK3 などの専用のコーデックが追加されました。コーデックは、htscodecs に実装されており、統計情報に基づいてカラム（データ系列）ごとに最適なものが利用されます。たとえば、CRAM3.1では品質スコアに FQZ が、リード名に TOK3 が使われます。
以前はHTSlibはデフォルトで欧州分子生物学研究所EBIのCRAM Reference Registryから参照配列を取得してました。しかし、サーバー・ネットワークへの負担が大きすぎるとしてEBIからのリファレンス配列の自動取得は廃止されました。今では各研究機関で ref-cache という付属のミニサーバーを起動してキャッシュ取得できるようにする方法が推奨されているようです。
# ref-cache を起動
ref-cache -b -d /data/ref_cache -l /data/logs -p 8080 \
-u https://www.ebi.ac.uk/ena/cram/md5/
サンプル数の増大に伴うIDのハッシュ衝突に対処する
- v1.21: ハッシュアルゴリズムをFNV1aに変更
サンプルIDなどを管理する際、プログラム内部ではハッシュ計算でIDを整理します。HTSlibはC言語で実装されており、ハッシュテーブルは Klib の khash で実装されています。しかし、サンプル数が数十万を超えたことで、偶然同じ整理番号になってしまう（ハッシュ衝突）事例が多発し、処理が重くなるケースがありました。このためアルゴリズムが X31 から FNV1a に変更されました。
クラウド（AWS S3）のための最適化
- v1.23: S3読み取りのチャンク化など
HTSファイルの多くは、BGZF圧縮によりランダムアクセスが可能になっています。クラウド上に置いたデータも同様にランダムアクセスできれば、通信量の減少につながります。htslibでは、AWS S3 などのオブジェクトストレージのデータの転送コストを下げるため、チャンク化しています。HTTP Rangeヘッダーを使用して必要なデータ範囲のみを取得します。
# 4MBチャンクに設定する
export HTS_S3_READ_PART_SIZE=4
samtools view s3://bucket/sample.cram
個人的にちょっと思うこと
HTSlibの開発の中心イギリスではクラウド上で大規模コホートを研究するためのCRAM推進が精力的に行われているようです。しかしローカルPCで孤軍奮闘する日本の零細現場ではCRAMは不便と考える人も少なくないのではないかと思います。
フェージングと3次元ゲノム
フェージング
ヒトは父由来・母由来の2本の染色体を持っています。従来の解析はこれを区別できないケースが多く、ごちゃまぜにして読んでいました。 Phasing（フェージング）とは、ある変異は父由来の染色体に乗っている、別の変異は母由来の染色体に乗っていると区別することです。
Prefixed Phasing（プレフィックス付きフェージング）はVCF 4.4で導入された機能で、ジェノタイプの最初のアレルにも明示的なフェージング記号（|または/）を付ける仕様です。それ以前は、最初のアレルのフェージング状態は他のアレルから推測されており曖昧でした。
#include <htslib/vcf.h>
// ジェノタイプ配列の準備（3サンプル、2倍体）
int32_t *gt_array = malloc(3 * 2 * sizeof(int32_t));
// サンプル1: |0|1 (両アレルがフェーズ済み)
gt_array[0] = bcf_gt_phased(0); // 最初のアレル: |0
gt_array[1] = bcf_gt_phased(1); // 2番目のアレル: |1
// サンプル2: /0/1 (両アレルが非フェーズ)
gt_array[2] = bcf_gt_unphased(0); // 最初のアレル: /0
gt_array[3] = bcf_gt_unphased(1); // 2番目のアレル: /1
// サンプル3: |0/1 (混合フェーズ)
gt_array[4] = bcf_gt_phased(0); // 最初のアレル: |0
gt_array[5] = bcf_gt_unphased(1); // 2番目のアレル: /1
// BCFレコードに更新
bcf_update_genotypes(hdr, record, gt_array, 3 * 2);
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE1 SAMPLE2 SAMPLE3
1 100 . A G . . . GT |0|1 /0/1 |0/1
四倍体であっても同じように対処できるそうです。
#include <htslib/vcf.h>
// 四倍体のジェノタイプ配列（2サンプル）
int32_t *gt_array = malloc(2 * 4 * sizeof(int32_t));
// サンプル1: |0|1|2|3 (すべてフェーズ済み)
gt_array[0] = bcf_gt_phased(0); // |0
gt_array[1] = bcf_gt_phased(1); // |1
gt_array[2] = bcf_gt_phased(2); // |2
gt_array[3] = bcf_gt_phased(3); // |3
// サンプル2: /0/1/2/3 (すべて非フェーズ)
gt_array[4] = bcf_gt_unphased(0); // /0
gt_array[5] = bcf_gt_unphased(1); // /1
gt_array[6] = bcf_gt_unphased(2); // /2
gt_array[7] = bcf_gt_unphased(3); // /3
// BCFレコードに更新
bcf_update_genotypes(hdr, record, gt_array, 2 * 4);
Number 記法
- v1.22: VCFv4.4/4.5のNumberフィールドサポート
htslibでは VCF/BCF の FORMAT の中の GT フィールドでジェノタイプを表します。しかし GT フィールド は特別扱いです。VCFヘッダーではString型として定義されていますが、BCFではInt型として格納されていたりしています。
しかし、そのジェノタイプに、さらにメタデータを付加したい場合があります。そのときに、以下のフィールドでメタデータの「区分」と「長さ」を指定してから、別のフィールドにメタデータを記録します。
-
BCF_VL_P(Number=P): GTで定義された各アレル値ごと -
BCF_VL_LA(Number=LA): LAAでリストされたaltアレルのみを考慮 -
BCF_VL_LG(Number=LG): LAAに基づいたgenotypeごと -
BCF_VL_LR(Number=LR): LAAに基づいた全アレルごと
※ LAA（Local ALT Alleles） とは「理論上あり得る ALT／Genotype の中から、このレコードの GT に 実際に現れたものだけ を抜き出した集合」のこと。そういう専用のタグがあるわけではない。
これは、Hi-Cなどの3次元ゲノム解析（3D Genomics）と関係しています。Hi-CはPhasingにも有効です。「父方・母方」という区別だけでなく、染色体上の長距離にわたる相互作用を含めた複雑な情報が含まれます。その結果をVCFに記述する際に、まず先ほどの BCF_VL_P (Number=P) でアレルの長さを記載して、以下の3つの FORMAT を記録します。
-
PSL(Phasing String List): 所属するハプロタイプブロックの名前（ID）。 -
PSO(Offset): そのブロックの中での位置。 -
PSQ(Quality): その情報の信頼度。
bcf_hdr_id2length() // Number記法の取得
bcf_hdr_id2number() // Number値の取得
6. トランスクリプトーム解析とスプライシング
RNA-seqも成熟した分野ですが、ロングリードの影響が出ています。
- v1.20: SAM配列タイプBタグの高速化（ONTやPacBioデータで頻出）
遺伝子（DNA）からRNAができるとき、不要な部分を切り貼りする「スプライシング」が起きます。 ショートリードでは、切り刻まれたRNA断片を読むため、元の全貌を推測する必要がありました。 しかし、ロングリードを使えば、RNAの端から端まで（全長）を一気に読めます。たとえばPacBio社の「Iso-Seq」という技術があります。がん細胞特有の異常なスプライシングなどがわかるようになりました。その際に、スプライス境界ごとのスコアを記録する必要があります。
B はSAM/BAMのタグの型指定子の一つですが、異なる複数のデータ型の配列を格納できます。ONTやPacBioのロングリードデータでシグナル品質や修正情報を格納するためによく使用されます。タイプ B を指定して異なる型を扱う処理は重くなります。htslibにこれを高速化する改善が行われました。
7. がんゲノムと臨床応用
多くのセキュリティが認められます
- v1.20-23: バッファオーバーフローやTLS周りの脆弱性修正が多数
これは、ゲノム解析が研究室のネットワークから、病院の検査システムやクラウド環境へ移行していることを物語っています。 患者さんの個人情報を含むゲノムデータを扱う以上、データの漏洩や改ざんを防がなければなりません。また、クリニカルシーケンスではツールの動作安定性が診断結果に直結します。
8. パフォーマンス・チューニング
最後に低レイヤーのパフォーマンスの改善です。
SIMD命令を使った低レイヤーな最適化や、VCFパーサーの高速化などが行われています。
まとめ
htslibの更新履歴を追ってみると、次のような傾向が見えてきます。
- ロングリードシーケンス技術・シングルセルシーケンス技術への対応の跡が見られます
- 増大するゲノムデータに対処するために、データ圧縮、ハッシュアルゴリズム、クラウド対応といった地道な改善が続けられています。
- 医療データを実用的なインフラとして利用するためセキュリティ修正が多く行われています。
ポエム
私たちは、どうしても BAM や BCF のようなファイルフォーマットは海の向こうの賢い人たちが策定した動かない仕様であり、天から与えられた完成品かのように考えがちです。
しかし、実際にはSAMやBAMといったファイル形式は、ある時代の研究グループが、当時の技術を背景に策定したファイルフォーマットに過ぎません。彼らは私たちと同じ人間、科学者なので限界もあります。
たとえば、今となっては修飾塩基情報は、塩基品質ように別のカラムとして扱う方が自然に感じられるかもしれません。また、uBAM/modBAMのようなBAMの使用法はマッピングされてないのですから、本来のBAMファイルの利用方法とは大きく異なります。こういった予想できない未来に対して、弾力性をもって対処できるように、HTSのファイル群の多くは、TAGで拡張できるようになっています。
また、もしも無制限の計算速度と計算資源があれば、そもそもHTSファイルを使わず、SQLiteやSQL言語などを使ったり、ARROW形式にしてメモリに乗せればいいわけですが、科学という営みは常に人間社会の限界に挑戦するものなので、いつでも特殊用途に最適化されたファイルフォーマットが普及します。この点については、これからも変わらないと思います。
人間は自然界から見たら本当に小さな存在で、個人はさらに小さいです。それでも決して探求をやめようとはしないでしょう。
この記事は以上です。