genomic-coordinates

ゲノム座標

使用する場面

ファイル形式間、ツール間、アセンブリ間、またはゲノムと転写産物の間で座標が境界を越えるときは、いつでも使用します。

ルール

座標は1つの値ではなく、数値、表記に使われている規約、そして測定対象のアセンブリという3つの情報です。 3つすべてを保持しなければ、その数値は解釈できません。

座標の誤りは、ゲノミクスにおける最も見つけにくいバグの一種です。オフバイワンのBEDファイルは、エラーなく読み込まれ、ソートされ、インターセクト処理されます。GRCh37のVCFをGRCh38のアノテーションと結合しても、行が返されます。右にシフトしたインデルはClinVarの該当エントリと一致せず、その結果、変異は新規と報告されます。エラーは何も発生しません。答えが単に間違っているだけで、しかももっともらしく見える方向に間違っています。

したがって、記憶に頼らず表を使って変換し、参照配列がある場合は必ず照合してください。

2つの変換

1-based inclusive  ->  0-based half-open :  start - 1,  end
0-based half-open  ->  1-based inclusive :  start + 1,  end

終了座標は決して変わりません。変換によって両方の数値が変わった場合は、誤りです。

各形式の座標規約

0ベース、半開区間	1ベース、両端を含む
BED、bedGraph、bigWig、narrowPeak	GFF3、GTF、VCF
BAM/CRAM（バイナリのPOS）	SAM（テキストのPOS）
PSL、genePred、refFlat	WIG、Picard interval_list
MAF（UCSCの複数アラインメント）	MAF（TCGAの変異アノテーション）
PyRanges、pybedtools	GRanges/IRanges、samtools、UCSC、Ensemblの領域文字列

「MAF」形式はどちらも実在しますが、それぞれ意味が異なり、互いに一致しません。UCSCは0ベースのファイルを、1ベースのブラウザ入力欄を通じて提供します。形式ごとの詳細を含む完全な表は、references/format-conventions.mdにあります。

cd skills/genomic-coordinates/scripts

python3 convert_coords.py --list                          # the table
python3 convert_coords.py --from bed --to gff chr1 999 1000
python3 convert_coords.py --from ucsc --to bed "chr7:5,530,601-5,530,625"
python3 convert_coords.py --from granges --to pyranges --input regions.tsv

contig  input                 output           length  status  detail
chr7    chr7:5530601-5530625  5530600-5530625  25      ok

ゼロ長BEDフィーチャー（chromStart == chromEnd。挿入位置として有効）は、end = start - 1に変換せず、unrepresentableとして報告されます。退化または無効な区間が1つでもある場合、終了コードは1になります。

変異は区間ではありません

インデルのVCF POSはアンカーベースであり、イベントの前にある、変化していない塩基です。また、同じ変化を複数の方法で表記できます。 chr1:7:CAC:C、chr1:3:CAC:C、chr1:2:GCA:Gは、いずれも同じ欠失です。正規化の前に変異を結合、重複排除、または検索すると、実際に一致するものが気付かれないまま失われます。インデルが集中する反復配列では、とりわけ多く失われます。

比較の前に正規化してください。最小表現になるようにトリミングし、その後、参照配列に対して左寄せします。

python3 normalize_variant.py --fasta ref.fa chr1 7 CAC C
python3 normalize_variant.py --fasta ref.fa --split --input cohort.vcf
python3 normalize_variant.py --fasta ref.fa --compare chr1:7:CAC:C chr1:2:GCA:G

input         normalized    type      pos_shift  ref_check  changed
chr1:7:CAC:C  chr1:2:GCA:G  deletion  5          ok         yes

各レコードのREFは、まずFASTAと照合されます。MISMATCHは、変異と参照配列のアセンブリが異なることを意味します。座標を調整せずに処理を止め、check_contigs.pyを実行してください。多対立遺伝子レコードは、正規化の前に--splitで分割する必要があります。後から分割してはいけません。

HGVSでは、インデルを転写産物に沿って反対方向、つまり3'側へシフトします。マイナス鎖の遺伝子では、これはVCFの左寄せとはゲノム上で反対の方向になります。詳細と完全な手順は、references/variant-representation.mdを参照してください。

結合を信頼する前にアセンブリを確認する

python3 check_contigs.py --identify unknown.fa.fai
python3 check_contigs.py variants.vcf annotation.gtf --genome GRCh38.fa.fai

file          kind    contigs  naming        assembly  detail
ref.fa.fai    sizes   25       plain         GRCh37    24/24 primary chromosome lengths match;
                                                       chrM is 16569 bp, i.e. GRCh37/38 (rCRS MT)

このスクリプトは.fai、.chrom.sizes、VCFヘッダー、SAMヘッダー、FASTA、BED、GTF/GFFを読み込み、主要染色体の長さからアセンブリを特定します。また、2つのファイルを結合したときに問題となる理由をすべて報告します。たとえば、名前の不一致、長さの不一致、コンティグ末端を超える座標、一方のファイルにしか存在しないコンティグなどです。互換性がない場合、終了コードは1になります。

GRCh37とhg19の違いはミトコンドリアだけです。長さは16,569 bp（rCRS）と16,571 bpです。核の座標は同一のため、混在したパイプラインも問題なく実行され、mtDNAの結果だけが誤ります。どちらが検出されたかはcheck_contigs.pyが報告します。ビルド、命名規則、ALTコンティグ、リフトオーバーに関する注意点は、references/reference-builds.mdを参照してください。

ファイルの形式に照らして監査する

python3 audit_intervals.py peaks.bed
python3 audit_intervals.py gencode.gtf --genome hg38.chrom.sizes
python3 audit_intervals.py cohort.vcf --genome GRCh38.fa.fai

座標の誤りが残す証拠を探します。

検出項目	何を示すか
GFF/GTF の start_below_one	1-based のファイルに 0-based データが含まれています。すべての座標が 1 塩基左にずれています
BED の many_zero_length	0-based のファイルに 1-based の1塩基フィーチャーが記述されています
past_contig_end	アセンブリが違うか、コンティグ末端で1塩基ずれています
mixed_contig_naming	結合すると、一部だけが黙って一致します
first_block_offset	BED12 の blockStarts に絶対座標が記述されています
not_parsimonious	アレルがトリミングされていません。結合前に正規化してください
bad_alt_allele	VCF におけるアンカーベースのない Ensembl/VEP の - 表記です

致命的な検出項目が1つでもあれば終了コード 1 を返すため、データディレクトリに対する CI ゲートとして機能します。

トランスクリプト、CDS、タンパク質の位置

c.742 と chr17:7,674,220 はどちらも「位置」ですが、算術演算で相互に変換することはできません。トランスクリプト座標は、スプライス後の塩基を転写方向に数えます。マイナス鎖ではゲノム座標が減少します。また、c.1 はトランスクリプトの開始位置ではなく、開始コドン ATG の A です。

記憶違いが起こりやすいルールは次のとおりです。c.0 はありません。5' UTR の位置は負数で、3' UTR の位置には * を付けます。GFF のフェーズは次のコドンに到達するために取り除く塩基数であり、start % 3 ではありません。また、同じバリアントでもトランスクリプトごとに番号が異なるため、バージョン付きトランスクリプトアクセッションがなければ c. 記述に意味はありません。変換手順と境界ケースについては、references/transcript-coordinates.md を参照してください。

トランスクリプトモデルを保持するツール（VEP、bcftools csq、Mutalyzer、hgvs パッケージ）を使って変換してください。手作業で変換しないでください。

結果の報告

座標の横には、毎回アセンブリを明記してください。 chr7:5,530,601-5,530,625 だけでは位置を特定できませんが、chr7:5,530,601-5,530,625 (GRCh38) なら特定できます。座標列のヘッダーまたはファイルのドキュメントに、座標の規約を明記してください。変換によって結果が得られた場合は、どちら向きに変換したかも明記してください。

参考資料

• references/format-conventions.md — 各形式の規約と形式ごとの詳細、BED12 のブロックルール、領域文字列の構文、ツールの動作。
• references/variant-representation.md — VCF のアレル規約、正規化アルゴリズム、同等性の確認、多アレルの分割、HGVS と VCF の不一致。
• references/reference-builds.md — ビルドの特徴、GRCh37 と hg19 の違い、ALT コンティグ、命名方式、リフトオーバーの失敗パターン。
• references/transcript-coordinates.md — ゲノム座標 ↔ トランスクリプト座標 ↔ CDS 座標 ↔ タンパク質座標、HGVS の番号付け、フェーズ、トランスクリプトの選択。

エージェントの作業手順

1. 環境を確認します。 ネットワークにアクセスできること、API キー（必要な場合）、データベースの最新の API ドキュメントとレート制限を確認します。
2. 入力を明確にします。 データまたはユーザーから、形式、識別子、パラメーターを確認します。結果に影響する値を推測せず、質問してください。
3. まず小規模に実行します。 既知のレコードを 1 件取得し、一括クエリの前に応答形式を確認します。
4. 上記の手順と参考情報を使って、タスク全体を実行します。
5. 結果を検証します。 識別子が解決できること、レコード数が報告されていること、データベースのバージョンまたはアクセス日が記録されていることを確認します。
6. 報告します。 実行した内容（バージョン、コマンド、パラメーター）、確認した内容、まだ不確かな点を説明します。

発生した場合	対応
HTTP 429 または 5xx エラーが発生した場合	バックオフを行い、リクエストをバッチ処理するなどしてレート制限を守り、部分的な結果を正直に報告します。
これらの手順に記載された関数、フラグ、エンドポイントがインストール済みバージョンに存在しない場合	インストール済みバージョンのドキュメント（help()、--help、公式ドキュメント）を確認して対応し、ユーザーに伝えます。API を捏造してはいけません。
必須の入力、識別子、パラメーターが曖昧な場合	ユーザーに質問するか、実行前に仮定を明示します。

整合性に関するルール

• 結果、パラメーター、識別子、引用、統計情報を決して捏造しないでください。実行または検証できないことは、はっきりと伝えてください。
• アクセッション番号、ID、レコードを捏造してはいけません。見つからない場合は「見つかりません」と報告します。
• バージョン固有の情報は古くなっている可能性があるものとして扱い、インストール済みバージョンの公式ドキュメントで確認してください。
• 金銭が発生する、共有 GPU やクラウドのクォータを消費する、個人データや患者データを扱う、または取り消せない操作を行う前には、ユーザーに確認してください。

関連スキル

• polars-bio：Polars DataFrame 上でゲノム領域の操作やバイオインフォマティクスのファイル入出力を行う Python ライブラリ polars-bio。Arrow と DataFusion を基盤としています。
• gtars：ゲノム領域向け Rust／Python／CLI ツールキット Gtars を使った作業を調査・計画します。BED RegionSet の集合演算（reduce、setdiff、intersec…）。
• geniml：Geniml（0.8.4）と Gtars を使ったローカルなゲノム領域機械学習ワークフローを計画・監査します。BED ファイルを染色体 si…に照らして検証します。