Zenn
🧬

IGVにゲノムアッセンブリーとアノテーションをマニュアルで登録しかつ遺伝子名を検索可能にする

2024/06/22に公開
Tips
バイオインフォマティクス
bioinformatics
ゲノム解析
IGV
tech

IGV version 2.17.4 (Windows) https://igv.org/doc/desktop

IGVでゲノムとアノテーションを自動で取得する

Integrative Genomics Viewer(IGV)は次世代シーケンサーのマッピング結果やアノテーション情報を可視化するツールである。RNA-seqのマッピング結果から遺伝子発現の度合いを確認したり、wgsの結果からSNPを探したり、ChiP-seqやATACの結果からシス制御配列の候補領域を絞ったり等、様々な研究に用いられている。

IGVでゲノム情報を確認するには、まず対象となるゲノムを用意する必要がある。これがヒトやマウスなどのモデル生物のリファレンスゲノムであれば既にIGV側がサーバーにゲノムやアノテーションを用意しており、プルダウンのメニューから選択するだけでDNA配列や遺伝子の位置がわかる状態になる。

Genomes -> Select Hosted Genome...

また最近のバージョンでは"Genomes"タブに"Load Genome from UCSC GenArk"が追加されている。これを使うことでNCBIのAssemblyデータベースに登録されていれば非モデル動物でもゲノムアッセンブリーとアノテーション一式を自動で読み込めるようになっている。

Genomes -> Load Genome from UCSC GenArk...

このようにIGVではモデル動物・非モデル動物に限らず手軽にリファレンスゲノムを読み込み、解析ができる。一方で、近年の次世代シーケンサー技術の発達と低価格化により、自分でゲノムアッセンブリーや遺伝子アノテーションを構築し解析する機会もあるだろう。このような場合、マニュアルでゲノムアッセンブリーとアノテーションをIGVに読み込ませることになる。マニュアル操作での読み込みはGUIで感覚的にできるようになっているが、いくつか問題がある。特に、アノテーションに載っている遺伝子を検索してそのlocusに飛びたいときにはIGVの仕様を知らないと検索ができないようになっている。こうした問題を解決するためには少しCUIの力を借りないといけない。

GUI操作でゲノムアッセンブリーとアノテーションを読み込む

データセット

ゲノムアッセンブリーとアノテーションを読み込む実際の操作を見ていく。今回はシーラカンスのゲノムアッセンブリー(.fna)と遺伝子アノテーション(.gtf)を実例として用いる。

ゲノムアッセンブリーの読み込み

ゲノムアッセンブリーのfastaファイルは"Genomes"のメニューから"Load Genome from File..."から読み込むことができる。

Genomes -> Load Genome from File...
本来fastaファイルを読み込むにはfasta indexファイルである.faiをsamtools等で作成する必要があるが、このようにGUIから読み込む場合には前もって.faiを生成する必要はない。IGVが自動でfastaファイルの直下に.faiを作る。

一度ゲノムアッセンブリーを読み込ませればゲノム一覧のメニューに登録され、再度読み出せるようになる。しかしこのようにゲノム一覧で表示される名前がfastaファイルの名前そのままで表示されて少々カッコ悪い。この登録されているアッセンブリーの名前は変えられないため、登録した名前が気に入らないのならば"Remove Genomes..."で一度消さなければならない。

Load Genome from File...でアッセンブリーを選択した直後

アノテーションファイルの読み込み

bedファイル、gff/gtfファイル、またはbamファイルも広義にはアノテーションと呼んでも良いかもしれない。アノテーションファイルには様々な形式がある。しかし大抵のものはドラッグ&ドロップで読み込めるようになっている。読み込ませるファイルが軽い場合、そのまま何事もなく表示される。しかしgtfのような全遺伝子アノテーション等の重いファイルに対しては以下のようなウィンドウが表示されてidxファイルを作ることを推奨される。

ドラッグ&ドロップ直後。ちなみにgff/gtfはソートされている必要があるので、ソートされていないのなら"Tools -> Run igvtools -> Sort"でソートする。

普通の感覚であれば何も考えずにGoを押すだろう。idxファイルを作成すると、次からの読み込み時間が格段に速くなる。

idxファイルを作成したアノテーションファイルには検索がかからない

しかしidxファイルを作成するデメリットは解析のときにある。我々はゲノムを漠然と眺めるのではなく、興味のある遺伝子周辺を調べる。前置きで書いた方法でリファレンスゲノムを自動で取得すれば検索欄で遺伝子の名前を入力すればその遺伝子のlocusまで移動する。だがマニュアルでアノテーションを読み込ませたとき、idxファイルを作成するとこの遺伝子名検索機能は使うことができない。


遺伝子名が検索可能であればこのように遺伝子名を入力すると候補名が吹き出しが出現し、遺伝子名を選択すればその遺伝子のlocusに飛ぶ

idxファイルが作成されていると、このように正しい遺伝子名を入力しても"Cannot find"となる

このバグについてgithubで制作者も回答しているが、これはバグではなく仕様のようである。おそらくアノテーションの読み込みを軽くするための仕様か何かと噛み合ってるせいだと思うが、しばらくこの仕様が変更されることはないだろう。
IGV imports features for searching from unsorted (and un-indexed) GFFs but not from indexed GFFs #244

アノテーションファイルに記載された遺伝子名で検索をかけたいのならば、idxファイルの作成を促す選択のときに"Cancel"を押せばこの仕様は適用されない。idxファイルが作成されないアノテーションファイルでも通常通りIGVで表示することができるが、ファイルの読み込みはかなり遅く、全ゲノムの遺伝子であればgff/gtfの読み込みに数分かかってしまう。

遺伝子アノテーションの表示

gff/gtfのアノテーションの読み込みに成功すればこのように遺伝子ごとのエキソン・イントロン構造が表示される。

今回のアノテーションファイルの場合、gtfの3列目の"feature"を示す列には

grep -v "#" GCF_037176945.1_fLatCha1.pri_genomic.gtf|awk '{print $3}'|sort |uniq
CDS
exon
gene
start_codon
stop_codon
transcript

の6つのfeatureが存在するが、このうち遺伝子のエキソン・イントロン構造の表示に関連するのはgene以外の5つのfeatureである。一方、geneのfeatureを持つ行も図中の真ん中の1本線のように表示されてしまうため邪魔である。

現在表示されている画面を保存する

読み込ませたアノテーションは保存されないため、一度IGVを閉じたり他のゲノムに切り替えると再度読み込ませなければならない。いちいち読み込ませるのは面倒であるが、現在読み込んでいるアノテーションや編集内容、現在表示しているlocus等を.xmlファイルとして保存する方法が実装されている。


File -> Save Session...

一度.xmlファイルとして保存すれば、"Open Session..."で.xmlファイルを読み込んだときに現在表示されている画面がそのまま再現される。

GUI操作の問題点

以上、GUI操作でゲノムアッセンブリーとアノテーションを読み込ませたときの問題点を３つ挙げた。

  • アッセンブリー名を編集できない
  • idxファイルの仕様
    • idxファイルを作成すると遺伝子名検索ができない
    • idxファイルを作成しないと読み込みに時間がかかる
  • gff/gtfに余計なfeatureが含まれている

これらの問題を一度に簡単に解決する方法はないが、個別に対処可能である。

JSON形式でアッセンブリーとアノテーションのセットを作る

JSON形式の設定ファイルの概要

GUI操作による読み込みの他に、IGVではJSON形式でゲノムの設定ファイルを書く方法も用意されている。IGVが用意しているゲノムとアノテーションのセットもJSON形式で書かれている。例えば最初からIGVに入っているhg38の設定ファイルは以下のように書かれている。

{
  "id": "hg38",
  "name": "Human (GRCh38/hg38)",
  "fastaURL": "https://s3.amazonaws.com/igv.broadinstitute.org/genomes/seq/hg38/hg38.fa",
  "indexURL": "https://s3.amazonaws.com/igv.broadinstitute.org/genomes/seq/hg38/hg38.fa.fai",
  "cytobandURL": "https://s3.amazonaws.com/igv.org.genomes/hg38/annotations/cytoBandIdeo.txt.gz",
  "aliasURL": "https://s3.amazonaws.com/igv.org.genomes/hg38/hg38_alias.tab",
  "tracks": [
    {
      "name": "Refseq Genes",
      "format": "refgene",
      "url": "https://s3.amazonaws.com/igv.org.genomes/hg38/ncbiRefSeq.sorted.txt.gz",
      "indexURL": "https://s3.amazonaws.com/igv.org.genomes/hg38/ncbiRefSeq.sorted.txt.gz.tbi",
      "visibilityWindow": -1,
      "removable": false,
      "order": 1000000
    }
  ],
  "chromosomeOrder": "chr1, chr2, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chr10, chr11, chr12, chr13, chr14, chr15, chr16, chr17, chr18, chr19, chr20, chr21, chr22, chrX, chrY"
}

自分で設定ファイルを作るときはこれを真似すればよい。Pythonのdictのように中括弧{}の中に{"field":"value", "field":"value", ...}とそれぞれのfield名とvalueを入れていくだけである。もちろん上記hg38の設定ファイルにあるすべてのfieldを使う必要はない。fieldとして必須なのは"id", "name", "fastaURL", "indexURL"の4項目のみである。詳しくはIGVのgithubで説明されている。
JSONファイルに関する説明:JSON Genome Format
fieldとして使える項目と意味:Reference Genome

ゲノムアッセンブリーの設定法

項目として必須な"id", "name", "fastaURL", "indexURL"の4項目はゲノムアッセンブリーを規定するためのものである。
id (たぶん)内部的な管理にのみ用いられるのでなんでもよい
name ゲノム一覧で実際に表示される名称
fastaURL ゲノムアッセンブリーの場所
indexURL fasta indexの場所

"name"でアッセンブリーに自由に名前をつけられるので、GUI操作の問題点として挙げたアッセンブリーの名称が変更できない問題はここで解決できる。"fastaURL"と"indexURL"はURLではなくPathでも良い。絶対Pathでも(JSONファイルのあるディレクトリから見た)相対Pathでも書ける。

例えば今用いているシーラカンスゲノムの登録は以下の通りである。

{
  "id":"GCF_037176945.1",
  "name":"Coelacanth (fLatCha1.pri)",
  "fastaURL": "./GCF_037176945.1_fLatCha1.pri_genomic.fna",
  "indexURL": "./GCF_037176945.1_fLatCha1.pri_genomic.fna.fai"
}

これを適当にcoelacanth.jsonと名前をつけてIGVに読み込ませる。
ファイルを読み込ませるときは"Genomes"->"Load Genome from File..."からcoelacanth.jsonを指定する。

GUI操作でゲノムアッセンブリーを読み込んだときと同様に、アッセンブリーに含まれている染色体やcontigが表示された。GUI操作時と違うのはゲノム名の表示が自分の指定した"Coelacanth (fLatCha1.pri)"になっているところである。

注意点としてはまずJSON形式なので"field":"value"同士を区切るコンマを忘れないこと。そして"}"で閉じる前の最後の"field":"value"にはコンマをつけないことである。これは非常に見落としやすい。そしてfastaファイルは圧縮したものは受け付けない。indexファイル(.fai)がない場合はsamtools faidx $fastaで予めindexを作っておく必要がある。

アノテーションの設定法

アノテーションは"tracks"のfieldで指定する。アノテーションには"name"や"format"、"url"等複数の項目を指定しなければならない。またアノテーションは１つとは限らない。そのため、"tracks"が指定する"value"は"tracks":[{アノテーション1の項目},{アノテーション2の項目},...]のように大括弧[]の中にアノテーションごとの情報を中括弧{}で区切った包含構造になる。アノテーションが1つであっても大括弧と中括弧で囲まなければならない。

上のシーラカンスゲノムのJSONファイルにアノテーションの情報を加えると以下のようになる。

{
  "id":"GCF_037176945.1",
  "name":"Coelacanth (fLatCha1.pri)",
  "fastaURL":"./GCF_037176945.1_fLatCha1.pri_genomic.fna",
  "indexURL":"./GCF_037176945.1_fLatCha1.pri_genomic.fna.fai",
  "tracks": [
    {
      "name":"RefSeq genes",
      "format":"gtf",
      "url":"./GCF_037176945.1_fLatCha1.pri_genomic.gtf",
      "indexURL":"./GCF_037176945.1_fLatCha1.pri_genomic.gtf.idx"
    }
  ]
}

tracks内の各アノテーションを規定する項目の説明:Tracks 2.0

アノテーションファイルに関して必須な項目は"name"と"url"のみである。"format"は無くても良い。無い場合は拡張子から自動判別される。

上記JSONファイルをIGVに読み込ませるとGUIでゲノムアッセンブリーとアノテーションを読み込ませたのと同じ状態となる。そのため問題点として挙げたアノテーションに関する問題は解決されない。

gff/gtfファイルの編集

アノテーションファイルの問題点として挙げたのは以下の2点だった。

  • idxファイルの仕様
    • idxファイルを作成すると遺伝子名検索ができない
    • idxファイルを作成しないと読み込みに時間がかかる
  • gff/gtfに余計なfeatureが含まれている

結論を言えばアノテーションファイルをエキソン・イントロン構造表示用のものと検索用のものとに分割し、前者にidxをつければ読み込みも速く遺伝子名の検索が可能となる。

今回用いたgtfファイルには遺伝子の"feature"としてCDS, exon, gene, start_codon, stop_codon, transcriptの6つが用いられており、エキソン・イントロン構造に用いられているのは"gene"以外の"feature"であった。また、"gene"のfeatureは表示上邪魔である。そこでfeatureが"gene"である行を取り除き、残りのfeatureのいずれかを持つ行をエキソン・イントロン構造表示用として残す。取り除いた"gene"の行には遺伝子名が記述されているため、これを検索用のファイルとして再利用する。ファイルサイズの圧倒的大きい前者はidxファイルが必要だが、後者は軽いので必要ない。

まずはfeatureが"gene"の行を取り除いた新しいgtfファイルを作成しidxを生成する。

awk '$3!="gene" {print}' GCF_037176945.1_fLatCha1.pri_genomic.gtf \
> GCF_037176945.1_fLatCha1.pri_genomic.new.gtf

igvtools index GCF_037176945.1_fLatCha1.pri_genomic.new.gtf

igvtoolsはIGVの"Tools -> Run igvtools..."にある機能をコマンドライン上で行うプログラムである。Download IGVのCommand line IGV and igvtoolsに入っているが、使わない人はGUI上でigvtoolsを走らせてindexをつければ良い。

次にfeatureが"gene"の行を抽出する。

awk '$3=="gene" {print}' GCF_037176945.1_fLatCha1.pri_genomic.gtf \
> GCF_037176945.1_fLatCha1.pri_genomic.genes.gtf

これらを新たにJSONファイルに書き込むと以下のようになる。

{
  "id":"GCF_037176945.1",
  "name":"Coelacanth (fLatCha1.pri)",
  "fastaURL":"./GCF_037176945.1_fLatCha1.pri_genomic.fna",
  "indexURL":"./GCF_037176945.1_fLatCha1.pri_genomic.fna.fai",
  "tracks": [
    {
      "name":"RefSeq genes",
      "format":"gtf",
      "url":"./GCF_037176945.1_fLatCha1.pri_genomic.new.gtf",
      "indexURL":"./GCF_037176945.1_fLatCha1.pri_genomic.new.gtf.idx"
    },
    {
      "name":"Genes",
      "format":"gtf",
      "url":"./GCF_037176945.1_fLatCha1.pri_genomic.genes.gtf",
      "hidden":true
    }
  ]
}

Genesのファイルは表示には邪魔なので"hidden":trueをつけると、画面上には表示されないが検索にはかかるようになる。検索用だとわかりやすいように"searchable":trueを明示しても良いが、これはデフォルトでtrueになる。

これでアノテーションファイルの読み込みが軽くなり遺伝子名も検索できるようになった。

Discussion