アメリエフのブログ

Amelieff Staff Blog

(Rで)VCFの読み込み "VariantAnnotation" ①

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を利用した、アノテーション付加などを紹介していきます。