GVFとVCFを編集する

## GVFとVCFとは

- `GVF` = `Genome Variant Format`
- `VCF` = `Variant Call Format`
- どちらもリファレンスゲノム配列に対して、種内変異の場所や種類を記載したフォーマット。
- どちらもよく使われるフォーマットだが、データベースやソフトウェアによって片方しか使われていないことがあるので、両方の編集ができるようになるといい。
- [NCBIの説明](https://www.ncbi.nlm.nih.gov/variation/tools/reporter/docs/examples)が詳しい。

## GVFを見てみよう

- 脊椎動物ゲノムに特化したデータベースである[Ensembl](https://www.ensembl.org)はGVFファイルで変異を記載している。
- ヒトに最も近縁な哺乳類である[チンパンジー](https://www.ensembl.org/Pan_troglodytes/Info/Index)のゲノムの種内変異を確認してみる。
- [Download all variants (GVF)](http://ftp.ensembl.org/pub/release-105/variation/gvf/pan_troglodytes/)のリンクから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](https://gatk.broadinstitute.org/)は、GVFファイルをサポートしていない。そのため、GATKなどで、EnsemblなどからダウンロードしたGVFの情報を使う場合には、事前にVCFに変換しないといけない。
- GVFとVCFを変換するツールはいろいろあるが、[RSAT](http://rsat.sb-roscoff.fr/convert-variations_form.cgi)が便利。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`　変異の座標
    - `ID`　ID（アクセッション番号など）
    - `REF`　リファレンスの塩基配列
    - `ALT`　変異の塩基配列
    - `QUAL`　変異の確からしさ（クオリティ）
    - `FILTER`　フィルター
    - `INFO`　変異に関する情報

## Indexの作成

- VCFファイルをGATKで使いたいときは、事前にIndexを作成しておこう。

```
gatk --java-options "-Xmx64G -Xms64G" \  # メモリサイズの指定
IndexFeatureFile -I pan_troglodytes.vcf  # Indexの作成
```

以上。（2022年1月5日、了。）