GVFとVCFとは
-
GVF=Genome Variant Format -
VCF=Variant Call Format - どちらもリファレンスゲノム配列に対して、種内変異の場所や種類を記載したフォーマット。
- どちらもよく使われるフォーマットだが、データベースやソフトウェアによって片方しか使われていないことがあるので、両方の編集ができるようになるといい。
- NCBIの説明が詳しい。
GVFを見てみよう
- 脊椎動物ゲノムに特化したデータベースであるEnsemblはGVFファイルで変異を記載している。
- ヒトに最も近縁な哺乳類であるチンパンジーのゲノムの種内変異を確認してみる。
- Download all variants (GVF)のリンクからHTTPでもFTPでもダウンロードできる。
wget ftp://ftp.ensembl.org/pub/release-105/variation/gvf/pan_troglodytes/pan_troglodytes.gvf.gz
zcat pan_troglodytes.gvf.gz | less # GZ圧縮されたGVFの中身を見る
$ ##gff-version 3
$ ##gvf-version 1.07
$ ##file-date 2021-08-24
$ ##genome-build ensembl Pan_tro_3.0
$ ##species http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=9598
$ ##feature-ontology $ http://song.cvs.sourceforge.net/viewvc/song/ontology/so.obo?revision=1.283
$ ##data-source Source=ensembl;version=105;url=http://vertebrates.ensembl.org/Pan_troglodytes
$ ##file-version 105
$ ##sequence-region AACZ04000307.1 1 22487
$ ##sequence-region AACZ04000323.1 1 76425
$ ##sequence-region AACZ04000444.1 1 32089
- GVFファイルの冒頭（ヘッダー）には、
##を行頭につけて情報が書かれている。 -
##sequence-regionに、contig/scaffold/chromosomeの配列のアクセッション番号と有効な塩基の開始位置と終了位置が書かれている。 - いくつ配列があるか（＝
##sequence-regionの行数）を数えてみよう。
grep '##sequence-region' pan_troglodytes.gvf \ # 配列情報の行を取り出す
> pan_troglodytes.gvf.list
head -5 pan_troglodytes.gvf.list
$ ##sequence-region AACZ04000307.1 1 22487
$ ##sequence-region AACZ04000323.1 1 76425
$ ##sequence-region AACZ04000444.1 1 32089
$ ##sequence-region AACZ04000445.1 1 1084
$ ##sequence-region AACZ04000446.1 1 30626
tail -5 pan_troglodytes.gvf.list
$ ##sequence-region 4 1 194502333
$ ##sequence-region 12 1 137163284
$ ##sequence-region 7 1 166211670
$ ##sequence-region Y 1 26350515
$ ##sequence-region MT 1 16554
wc -l pan_troglodytes.gvf.list # 行数を数える
$ 44449 pan_troglodytes.gvf.list
- なので、このチンパンジーゲノム
Pan_tro_3.0には、44449個の配列があることがわかる。 - ヘッダーが終わると、1行に1個ずつ、変異情報が記載される。
head -44463 pan_troglodytes.gvf | tail -9
$ ##sequence-region 7 1 166211670
$ ##sequence-region Y 1 26350515
$ ##sequence-region MT 1 16554
$ AACZ04000508.1 dbSNP SNV 855 855 . + . ID=1;Variant_seq=T;Dbxref=dbSNP_136:rs25750188;Reference_seq=C
$ AACZ04000572.1 dbSNP SNV 2004 2004 . + . ID=2;Variant_seq=A;Dbxref=dbSNP_136:rs25750606;Reference_seq=C
$ AACZ04000752.1 dbSNP SNV 690 690 . + . ID=3;Variant_seq=C;Dbxref=dbSNP_136:rs26756649;Reference_seq=A
$ AACZ04000752.1 dbSNP SNV 1044 1044 . + . ID=4;Variant_seq=T;Dbxref=dbSNP_136:rs26166753;Reference_seq=G
$ AACZ04000752.1 dbSNP SNV 1078 1078 . + . ID=5;Variant_seq=G;Dbxref=dbSNP_136:rs26731541;Reference_seq=A
$ AACZ04000752.1 dbSNP SNV 1132 1132 . + . ID=6;Variant_seq=C;Dbxref=dbSNP_136:rs24952187;Reference_seq=T
-
たとえば、変異の一行目である
AACZ04000508.1 dbSNP SNV 855 855 . + . ID=1;Variant_seq=T;Dbxref=dbSNP_136:rs25750188;Reference_seq=Cは次のように解釈する。- 配列名
AACZ04000508.1 - データベース
dbSNP - 変異の種類
SNV（※Single Nucleotide Variant：一塩基多型のこと） - 変異の開始位置
855 - 変異の終了位置
855 - 方向
+鎖 - ID
1 - 変異の塩基配列（Variant_seq）
T - DbxrefのID
dbSNP_136 - 変異のアクセッション番号
rs25750188 - リファレンスの塩基配列（Reference_seq）
C
- 配列名
-
このGVFの終わりの行を見てみる。
tail -5 pan_troglodytes.gvf
$ Y dbSNP SNV 26328961 26328961 . + . ID=1539131;Variant_seq=A;Dbxref=dbSNP_136:rs80465449;Reference_seq=C
$ Y dbSNP SNV 26329102 26329102 . + . ID=1539132;Variant_seq=T;Dbxref=dbSNP_136:rs80477922;Reference_seq=C
$ Y dbSNP SNV 26329180 26329180 . + . ID=1539133;Variant_seq=A;Dbxref=dbSNP_136:rs80474295;Reference_seq=G
$ Y dbSNP SNV 26329184 26329184 . + . ID=1539134;Variant_seq=A;Dbxref=dbSNP_136:rs80445970;Reference_seq=G
$ Y dbSNP SNV 26329213 26329213 . + . ID=1539135;Variant_seq=T;Dbxref=dbSNP_136:rs80553644;Reference_seq=G
- 最後の行はY染色体のSNV（一塩基多型）が書かれており、
ID=1539135で終わりなので、このIDは連番でついているから、1539135個の変異情報がこのGVFに記載されているのがわかる。（ちなみにこのGVFファイルの全行数は1583592で、headerの行数は44457なので、1583592 - 44457 = 1539135で間違いない。）
VCFに変換してみよう
-
VCF（Variant Call Format）もヘッダーで対応するゲノム情報を記載して、そのあとに一行に一個ずつ変異の情報を書いていく。書き方の形式が異なるだけで、GVFと互いに変換可能である。 - 変異解析を行う定番のソフトウェアであるGATKは、GVFファイルをサポートしていない。そのため、GATKなどで、EnsemblなどからダウンロードしたGVFの情報を使う場合には、事前にVCFに変換しないといけない。
- GVFとVCFを変換するツールはいろいろあるが、RSATが便利。Web上で扱えるので、GVFファイルをアップロードすれば、変換後、メールでダウンロードリンクを送ってもらえる。
-
RSATでpan_troglodytes.gvfを変換すると次の通り。
less pan_troglodytes.vcf
$ ##fileformat=VCFv4.1
$ ##RSAT; Phased False
$ ##RSAT; Homozygotes False
$ ##gff-version 3
$ ##gvf-version 1.07
$ ##file-date 2020-12-08
$ ##genome-build ensembl Pan_tro_3.0
$ ##species http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=9598
$ ##feature-ontology http://song.cvs.sourceforge.net/viewvc/song/ontology/so.obo?revision=1.283
$ ##data-source Source=ensembl;version=103;url=http://vertebrates.ensembl.org/Pan_troglodytes
$ ##file-version 103
$ ##sequence-region AACZ04000307.1 1 22487
$ ##sequence-region AACZ04000323.1 1 76425
$ ##sequence-region AACZ04000444.1 1 32089
- ヘッダーから始まるのはGVFと同じである。
head -44466 pan_troglodytes.vcf | tail -9
$ ##sequence-region 7 1 166211670
$ ##sequence-region Y 1 26350515
$ ##sequence-region MT 1 16554
$ #CHROM POS ID REF ALT QUAL FILTER INFO
$ AACZ04000508.1 855 rs25750188 C T . . TSA=SNV
$ AACZ04000572.1 2004 rs25750606 C A . . TSA=SNV
$ AACZ04000752.1 690 rs26756649 A C . . TSA=SNV
$ AACZ04000752.1 1044 rs26166753 G T . . TSA=SNV
$ AACZ04000752.1 1078 rs26731541 A G . . TSA=SNV
- ヘッダーが終わると、一行に一個ずつ変異の情報が書かれる。
- それぞれの列の意味は、ヘッダーの最後の行
#CHROM POS ID REF ALT QUAL FILTER INFOの通りである。具体的には次の通り。-
#CHROM配列名 -
POS変異の座標 -
IDID（アクセッション番号など） -
REFリファレンスの塩基配列 -
ALT変異の塩基配列 -
QUAL変異の確からしさ（クオリティ） -
FILTERフィルター -
INFO変異に関する情報
-
Indexの作成
- VCFファイルをGATKで使いたいときは、事前にIndexを作成しておこう。
gatk --java-options "-Xmx64G -Xms64G" \ # メモリサイズの指定
IndexFeatureFile -I pan_troglodytes.vcf # Indexの作成
以上。（2022年1月5日、了。）