VCF(Variant Call Format)をRで読みましょう！
Rパッケージ "VariantAnnotation" をご紹介します。
VariantAnnotationはVCFを読み込み、更にアノテーションパッケージを利用してアノテーションを付与することが可能です。
列の分類やPOS表記に少しクセがあり、慣れるのに少し時間がかかります。
ですので、単にVCFファイルの中身を見たいだけならばこのパッケージを使う必要はありません。
VCFやアノテーションの情報から変異を分類して、表やグラフを描きたいときには便利ですね。
今回は、インストール方法、VCFの読み込みと中身の閲覧の方法をご紹介します。 公式説明書はこちらです。
インストール
source("https://bioconductor.org/biocLite.R")
biocLite("VariantAnnotation")
# 以下は依存パッケージ
biocLite("GenomeInfoDb")
biocLite("BiocGenerics")
biocLite("GenomicRanges")
biocLite("SummarizedExperiment")
biocLite("Rsamtools")
VCFの読み込み
# ライブラリの読み込み
library(VariantAnnotation)
# VCFの読み込み
vcf <- readVcf("/path/to/ERX10000.vcf", "hg19") # readVcf(file, genome, param, ...)
読み込んだVCFの中身の閲覧 ~ヘッダー
読み込んだVCFは、R内では列ごとの表（DataFrame）が組み合わされたようなクラス CollapsedVCF として格納されています。
クラスやヘッダー（列名）を確認してみましょう。
基本情報 'ID, REF, ALT, QUAL, FILTER' や、INFO列に'AC ... Allele Count' などがあることが分かります。
> class(vcf) # クラスを確認
[1] "CollapsedVCF"
attr(,"package")
[1] "VariantAnnotation"
> vcf # クラスやヘッダーを確認
class: CollapsedVCF
dim: 82979 1
rowRanges(vcf):
GRanges with 5 metadata columns: paramRangeID, REF, ALT, QUAL, FILTER
info(vcf):
DataFrame with 18 columns: AC, AF, AN, BaseQRankSum, ClippingRankSum, D...
info(header(vcf)):
Number Type Description
AC A Integer Allele count in genotypes, for each ALT...
AF A Float Allele Frequency, for each ALT allele, ...
:
geno(vcf):
SimpleList of length 5: GT, AD, DP, GQ, PL
geno(header(vcf)):
Number Type Description
GT 1 String Genotype
AD R Integer Allelic depths for the ref and alt alleles in the or...
DP 1 Integer Approximate read depth (reads with MQ=255 or with ba...
:
上記の geno(vcf) などを使うと、
ヘッダー、INFO列やサンプル列項目の詳細を個別に見ることもできます。
> header(vcf) # ヘッダー一覧 > info(header(vcf)) # INFO列の項目詳細 > geno(vcf) # サンプル列の項目一覧 > geno(header(vcf)) # サンプル列の項目詳細
読み込んだVCFの中身の閲覧 ~中身
列名（ref, infoなど）を指定すると、
それぞれはリストや表なので、head() や[1:5,] などおなじみのRコマンドで好きな変異情報をみることができます。
例えば [n] で n番目の変異を指定できます。
> alt(vcf)[3] #3行目の変異のALT > ref(vcf)[3] #3行目の変異のREF > rowRanges(vcf)[3] # 3行目の変異のparamRangeID, REF, ALT, QUAL, FILTER
INFO列は、変異が行に、各項目が列の表（DataFrame）になっています。
こちらは列名（項目名）や列番号も指定しましょう。
> info(vcf) DataFrame with 82979 rows and 18 columns ... > info(vcf)$AF # INFOのAF列 > info(vcf)$AF[3] # INFOのAF列 3行目 > info(vcf)[1:4,1:3] # INFOの1~4列、1~3行
サンプル列はリストになっています。 これも項目名と変異の番号を指定しましょう。
> geno(vcf) List of length 5 names(5): GT AD DP GQ PL > geno(vcf)$GQ[1:5] #GQの1~5番目の変異
列によって形式がリストやデータフレームと異なっているため、冒頭に書いた通り、慣れるのに時間がかかります。
覚えるのは難しいですが、 > geno(vcf) のように一度見たい情報を画面に打ってみればそのデータ形式が分かります。
大元のファイルは変化しないので、心配なく何でも試してみましょう！
今後、VariantAnnotationを利用した、アノテーション付加などを紹介していきます。
