RでVCFファイルを操作する

RでVCFファイルを扱う場合、CRANにある vcfR または BioconductorVariantAnnotation というパッケージの2つの選択肢があると思います。今回は BioconductorVariantAnnotation を触ってみたのでその記録を残しておきたいと思います。

Bioconductor パッケージなのでインストールは以下のようにしました。

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("VariantAnnotation")

まずはパッケージをロードします。

library(VariantAnnotation)

今回はパッケージ付属のVCFファイルを使ってみます。

# パッケージ付属のVCFファイルに何があるか確認
dir(system.file("extdata", package = "VariantAnnotation"))
#>  [1] "chr22.vcf.gz"                  "chr22.vcf.gz.tbi"              "chr7-sub.vcf.gz"               "chr7-sub.vcf.gz.tbi"          
#>  [5] "ex2.vcf"                       "gl_chr1.vcf"                   "h1187-10k.vcf.gz"              "h1187-10k.vcf.gz.tbi"         
#>  [9] "hapmap_exome_chr22.vcf.gz"     "hapmap_exome_chr22.vcf.gz.tbi" "structural.vcf"

# chr22.vcf.gzを使う(今回のコードの一部は同じディレクトリに.tbiという拡張子のインデックスファイルがあることが必要)
fl <- system.file("extdata", "chr22.vcf.gz", package = "VariantAnnotation")

まずはファイルの読み込み。

vcf <- readVcf(fl)
vcf
#> class: CollapsedVCF 
#> dim: 10376 5 
#> rowRanges(vcf):
#>   GRanges with 5 metadata columns: paramRangeID, REF, ALT, QUAL, FILTER
#> info(vcf):
#>   DataFrame with 22 columns: LDAF, AVGPOST, RSQ, ERATE, THETA, CIEND, CIPOS, END, HOMLEN, HOMSEQ, SVLEN, SVTYPE, AC, AN, AA, AF, ...
#> info(header(vcf)):
#>              Number Type    Description                                                                                             
#>    LDAF      1      Float   MLE Allele Frequency Accounting for LD                                                                  
#>    AVGPOST   1      Float   Average posterior probability from MaCH/Thunder                                                         
#>    RSQ       1      Float   Genotype imputation quality from MaCH/Thunder                                                           
#>    ERATE     1      Float   Per-marker Mutation rate from MaCH/Thunder                                                              
#>    THETA     1      Float   Per-marker Transition rate from MaCH/Thunder                                                            
#>    CIEND     2      Integer Confidence interval around END for imprecise variants                                                   
#>    CIPOS     2      Integer Confidence interval around POS for imprecise variants                                                   
#>    END       1      Integer End position of the variant described in this record                                                    
#>    HOMLEN    .      Integer Length of base pair identical micro-homology at event breakpoints                                       
#>    HOMSEQ    .      String  Sequence of base pair identical micro-homology at event breakpoints                                     
#>    SVLEN     1      Integer Difference in length between REF and ALT alleles                                                        
#>    SVTYPE    1      String  Type of structural variant                                                                              
#>    AC        .      Integer Alternate Allele Count                                                                                  
#>    AN        1      Integer Total Allele Count                                                                                      
#>    AA        1      String  Ancestral Allele, ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/technical/reference/ancestral_a...
#>    AF        1      Float   Global Allele Frequency based on AC/AN                                                                  
#>    AMR_AF    1      Float   Allele Frequency for samples from AMR based on AC/AN                                                    
#>    ASN_AF    1      Float   Allele Frequency for samples from ASN based on AC/AN                                                    
#>    AFR_AF    1      Float   Allele Frequency for samples from AFR based on AC/AN                                                    
#>    EUR_AF    1      Float   Allele Frequency for samples from EUR based on AC/AN                                                    
#>    VT        1      String  indicates what type of variant the line represents                                                      
#>    SNPSOURCE .      String  indicates if a snp was called when analysing the low coverage or exome alignment data                   
#> geno(vcf):
#>   SimpleList of length 3: GT, DS, GL
#> geno(header(vcf)):
#>       Number Type   Description                      
#>    GT 1      String Genotype                         
#>    DS 1      Float  Genotype dosage from MaCH/Thunder
#>    GL G      Float  Genotype Likelihoods

続いて header 情報を見てみます。

header(vcf)
#> class: VCFHeader 
#> samples(5): HG00096 HG00097 HG00099 HG00100 HG00101
#> meta(1): fileformat
#> fixed(2): FILTER ALT
#> info(22): LDAF AVGPOST ... VT SNPSOURCE
#> geno(3): GT DS GL
samples(header(vcf))
#> [1] "HG00096" "HG00097" "HG00099" "HG00100" "HG00101"

VCFの1から7列目の情報を見てみます。

rowRanges(vcf)
#> GRanges object with 10376 ranges and 5 metadata columns:
#>               seqnames    ranges strand | paramRangeID            REF                ALT      QUAL      FILTER
#>                  <Rle> <IRanges>  <Rle> |     <factor> <DNAStringSet> <DNAStringSetList> <numeric> <character>
#>     rs7410291       22  50300078      * |         <NA>              A                  G       100        PASS
#>   rs147922003       22  50300086      * |         <NA>              C                  T       100        PASS
#>   rs114143073       22  50300101      * |         <NA>              G                  A       100        PASS
#>   rs141778433       22  50300113      * |         <NA>              C                  T       100        PASS
#>   rs182170314       22  50300166      * |         <NA>              C                  T       100        PASS
#>           ...      ...       ...    ... .          ...            ...                ...       ...         ...
#>   rs187302552       22  50999536      * |         <NA>              A                  G       100        PASS
#>     rs9628178       22  50999538      * |         <NA>              A                  G       100        PASS
#>     rs5770892       22  50999681      * |         <NA>              A                  G       100        PASS
#>   rs144055359       22  50999830      * |         <NA>              G                  A       100        PASS
#>   rs114526001       22  50999964      * |         <NA>              G                  C       100        PASS
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

最初の6つの Variant の Genotype を見てみます。

head(geno(vcf)$GT)
#>             HG00096 HG00097 HG00099 HG00100 HG00101
#> rs7410291   "0|0"   "0|0"   "1|0"   "0|0"   "0|0"  
#> rs147922003 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"  
#> rs114143073 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"  
#> rs141778433 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"  
#> rs182170314 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"  
#> rs115145310 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"

VCFファイルのサイズが大きすぎてすべてをメモリに収められないときは少しずつ読み込んでいくということもできます。

# VCFを10行ずつ読み込む
tf <- TabixFile(fl, yieldSize = 1000)
open(tf)
while (TRUE) {
  vcf <- readVcf(tf)
  # 何かしらのコード
  if (nrow(vcf) == 0) break
}
close(tf)

他にも色々機能があると思いますが、今回はここまでにしておきます。 詳細は Bioconductor のページを参照しましょう。