RでVCFファイルを扱う場合、CRANにある vcfR
または Bioconductor の VariantAnnotation
というパッケージの2つの選択肢があると思います。今回は Bioconductor の VariantAnnotation
を触ってみたのでその記録を残しておきたいと思います。
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 のページを参照しましょう。