RでVCFファイルを操作する2ということで2ヶ月前の記事に引き続き VariantAnnotation
でのVCFの扱い方について説明していこうと思います。
今回はデータのプロットまでしてみようと思います。
データインポート
前回と同様、パッケージに同梱されているVCFファイルを用います。VCFファイルを readVcf 関数で読み込んでいます。
VCFファイルのサイズがでかいとメモリに収まりきらない場合があるので比較的小さいもしくはフィルタリングされたVCFを扱う場合で説明します。
library(VariantAnnotation) library(ggplot2) 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,... ##> 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 var... ##> CIPOS 2 Integer Confidence interval around POS for imprecise var... ##> END 1 Integer End position of the variant described in this re... ##> HOMLEN . Integer Length of base pair identical micro-homology at ... ##> HOMSEQ . String Sequence of base pair identical micro-homology a... ##> 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.u... ##> AF 1 Float Global Allele Frequency based on AC/AN ##> AMR_AF 1 Float Allele Frequency for samples from AMR based on A... ##> ASN_AF 1 Float Allele Frequency for samples from ASN based on A... ##> AFR_AF 1 Float Allele Frequency for samples from AFR based on A... ##> EUR_AF 1 Float Allele Frequency for samples from EUR based on A... ##> VT 1 String indicates what type of variant the line represents ##> SNPSOURCE . String indicates if a snp was called when analysing the... ##> 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
コンソール上でVCFデータを格納した変数を叩くとVCFの情報を閲覧することができます。
dim: で示されている数字はそれぞれ行と列の大きさで、今回の場合は 10376 のバリアントが 5
サンプルについてコールされたデータとなります。
それ以降は、VCF内の各データにアクセスするための関数とその関数で得られるデータの説明が見られます。
データアクセス
実際にVCFのデータにアクセスしてみましょう。
VCFの中身を見たことがある方はご存知かと思いますがVCFは基本的に表形式のデータです。 各行に変異のデータが記載されています。
左から7列は各変異の座乗位置やID,リファレンスのアレル型とそれに対応する対立のアレル型、クオリティとフィルター情報になります。
これらの情報はどんなVCFにも必ず含まれている情報で、 rowRanges
関数でアクセスすることができます。
また、fixed 関数は座乗位置情報とID（最初の3列）を除く4列のデータを得ることができます。
いずれも表形式のデータとして取得することになります。
また、各列のみを抽出するような関数（seqnames, ref, alt,
qual)もあります。これらにより返されるオブジェクトのクラスは
様々なのでそれぞれの方法で対応する必要があります。得にALT列は1行に複数のデータが含まれている場合があり、それも踏まえた形でデータが返されます。
rowRanges(vcf) ##> GRanges object with 10376 ranges and 5 metadata columns: ##> seqnames ranges strand | paramRangeID REF ##> <Rle> <IRanges> <Rle> | <factor> <DNAStringSet> ##> rs7410291 22 50300078 * | <NA> A ##> rs147922003 22 50300086 * | <NA> C ##> rs114143073 22 50300101 * | <NA> G ##> rs141778433 22 50300113 * | <NA> C ##> rs182170314 22 50300166 * | <NA> C ##> ... ... ... ... . ... ... ##> rs187302552 22 50999536 * | <NA> A ##> rs9628178 22 50999538 * | <NA> A ##> rs5770892 22 50999681 * | <NA> A ##> rs144055359 22 50999830 * | <NA> G ##> rs114526001 22 50999964 * | <NA> G ##> ALT QUAL FILTER ##> <DNAStringSetList> <numeric> <character> ##> rs7410291 G 100 PASS ##> rs147922003 T 100 PASS ##> rs114143073 A 100 PASS ##> rs141778433 T 100 PASS ##> rs182170314 T 100 PASS ##> ... ... ... ... ##> rs187302552 G 100 PASS ##> rs9628178 G 100 PASS ##> rs5770892 G 100 PASS ##> rs144055359 A 100 PASS ##> rs114526001 C 100 PASS ##> ------- ##> seqinfo: 1 sequence from an unspecified genome; no seqlengths fixed(vcf) ##> DataFrame with 10376 rows and 4 columns ##> REF ALT QUAL FILTER ##> <DNAStringSet> <DNAStringSetList> <numeric> <character> ##> 1 A G 100 PASS ##> 2 C T 100 PASS ##> 3 G A 100 PASS ##> 4 C T 100 PASS ##> 5 C T 100 PASS ##> ... ... ... ... ... ##> 10372 A G 100 PASS ##> 10373 A G 100 PASS ##> 10374 A G 100 PASS ##> 10375 G A 100 PASS ##> 10376 G C 100 PASS
続いてVCFの8列目、INFO列のデータを取得してみます。この列についても各バリアントごとのデータを格納していますが、 その内容については各VCFファイルで異なり、どのようなデータが入っているかはヘッダ行の “##INFO=”に記載されています。
これらの情報は以下の様に取得できます。
info(vcf) ##> DataFrame with 10376 rows and 22 columns ##> LDAF AVGPOST RSQ ERATE THETA CIEND ##> <numeric> <numeric> <numeric> <numeric> <numeric> <IntegerList> ##> rs7410291 0.3431 0.989 0.9856 0.002 5e-04 NA,NA ##> rs147922003 0.0091 0.9963 0.8398 5e-04 0.0011 NA,NA ##> rs114143073 0.0098 0.9891 0.5919 7e-04 8e-04 NA,NA ##> rs141778433 0.0062 0.995 0.6756 9e-04 3e-04 NA,NA ##> rs182170314 0.0041 0.9981 0.7909 7e-04 4e-04 NA,NA ##> ... ... ... ... ... ... ... ##> rs187302552 9e-04 0.9992 0.5571 3e-04 0.0026 NA,NA ##> rs9628178 0.0727 0.9997 0.9983 3e-04 0.0011 NA,NA ##> rs5770892 0.3678 0.9868 0.9784 0.0021 7e-04 NA,NA ##> rs144055359 0.0011 0.9987 0.5323 5e-04 4e-04 NA,NA ##> rs114526001 0.0543 0.9622 0.7595 0.0021 5e-04 NA,NA ##> CIPOS END HOMLEN HOMSEQ SVLEN ##> <IntegerList> <integer> <IntegerList> <CharacterList> <integer> ##> rs7410291 NA,NA NA NA NA NA ##> rs147922003 NA,NA NA NA NA NA ##> rs114143073 NA,NA NA NA NA NA ##> rs141778433 NA,NA NA NA NA NA ##> rs182170314 NA,NA NA NA NA NA ##> ... ... ... ... ... ... ##> rs187302552 NA,NA NA NA NA NA ##> rs9628178 NA,NA NA NA NA NA ##> rs5770892 NA,NA NA NA NA NA ##> rs144055359 NA,NA NA NA NA NA ##> rs114526001 NA,NA NA NA NA NA ##> SVTYPE AC AN AA AF AMR_AF ##> <character> <IntegerList> <integer> <character> <numeric> <numeric> ##> rs7410291 NA 751 2184 N 0.34 0.2 ##> rs147922003 NA 20 2184 c 0.01 0.01 ##> rs114143073 NA 20 2184 G 0.01 0.0028 ##> rs141778433 NA 12 2184 C 0.01 0.01 ##> rs182170314 NA 8 2184 C 0.0037 0.01 ##> ... ... ... ... ... ... ... ##> rs187302552 NA 1 2184 a 5e-04 NA ##> rs9628178 NA 158 2184 a 0.07 0.03 ##> rs5770892 NA 801 2184 a 0.37 0.32 ##> rs144055359 NA 1 2184 g 5e-04 NA ##> rs114526001 NA 113 2184 g 0.05 0.01 ##> ASN_AF AFR_AF EUR_AF VT SNPSOURCE ##> <numeric> <numeric> <numeric> <character> <CharacterList> ##> rs7410291 0.19 0.83 0.22 SNP LOWCOV ##> rs147922003 NA 0.02 0.01 SNP LOWCOV ##> rs114143073 0.02 0.01 0.01 SNP LOWCOV ##> rs141778433 NA 0.02 NA SNP LOWCOV ##> rs182170314 NA 0.01 NA SNP LOWCOV ##> ... ... ... ... ... ... ##> rs187302552 0.0017 NA NA SNP LOWCOV ##> rs9628178 0.01 0.17 0.08 SNP LOWCOV ##> rs5770892 0.38 0.59 0.23 SNP LOWCOV ##> rs144055359 NA NA 0.0013 SNP LOWCOV ##> rs114526001 0.01 0.16 0.04 SNP LOWCOV info(header(vcf)) ##> DataFrame with 22 rows and 3 columns ##> Number Type ##> <character> <character> ##> LDAF 1 Float ##> AVGPOST 1 Float ##> RSQ 1 Float ##> ERATE 1 Float ##> THETA 1 Float ##> ... ... ... ##> ASN_AF 1 Float ##> AFR_AF 1 Float ##> EUR_AF 1 Float ##> VT 1 String ##> SNPSOURCE . String ##> Description ##> <character> ##> LDAF MLE Allele Frequency Accounting for LD ##> AVGPOST Average posterior probability from MaCH/Thunder ##> RSQ Genotype imputation quality from MaCH/Thunder ##> ERATE Per-marker Mutation rate from MaCH/Thunder ##> THETA Per-marker Transition rate from MaCH/Thunder ##> ... ... ##> ASN_AF Allele Frequency for samples from ASN based on AC/AN ##> AFR_AF Allele Frequency for samples from AFR based on AC/AN ##> EUR_AF Allele Frequency for samples from EUR based on AC/AN ##> VT indicates what type of variant the line represents ##> SNPSOURCE indicates if a snp was called when analysing the low coverage or exome alignment data
エクセルなどでVCFファイルを開くとINFO列は1つの列として表示されますが、上記出力からもわかるように、Rを用いれば、INFO列を的確にパースして、複数の列に分解してくれます。
最後に最も重要な実際のバリアントコールのデータを取得してみます。これについては geno 関数で取得することができます。
VCFにもよりますが、一般的に、バリアントコールのデータ（10行目以降）にはコールされた遺伝子型の他に
Likelihood などが格納されています。
そのため、geno 関数によって返される値も複数の表形式データの“リスト”になります。
多くの場合、遺伝子型のデータはGTという名前が振ってあるのでリストの要素にアクセスする記号
$ を用いて以下のようにアクセスてきます。 また、今回のデータでは10列目以降に含まれるデータは GT, DS, GL
となっていますが、それぞれのデータが意味することは以下のようにしてみることができます。
geno(vcf) ##> List of length 3 ##> names(3): GT DS GL geno(header(vcf)) ##> DataFrame with 3 rows and 3 columns ##> Number Type Description ##> <character> <character> <character> ##> GT 1 String Genotype ##> DS 1 Float Genotype dosage from MaCH/Thunder ##> GL G Float Genotype Likelihoods 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のデータをフィルタリングしたい場合があると思います。その場合は subset 関数を用いることが一つの解決法になります。
以下に例を示します。
# SNVだけを取得する fixed(subset(vcf, isSNV(vcf))) ##> DataFrame with 9969 rows and 4 columns ##> REF ALT QUAL FILTER ##> <DNAStringSet> <DNAStringSetList> <numeric> <character> ##> 1 A G 100 PASS ##> 2 C T 100 PASS ##> 3 G A 100 PASS ##> 4 C T 100 PASS ##> 5 C T 100 PASS ##> ... ... ... ... ... ##> 9965 A G 100 PASS ##> 9966 A G 100 PASS ##> 9967 A G 100 PASS ##> 9968 G A 100 PASS ##> 9969 G C 100 PASS # InDelだけを取得する fixed(subset(vcf, isIndel(vcf))) ##> DataFrame with 406 rows and 4 columns ##> REF ALT QUAL FILTER ##> <DNAStringSet> <DNAStringSetList> <numeric> <character> ##> 1 C CT 425 PASS ##> 2 CA C 298 PASS ##> 3 G GA 497 PASS ##> 4 CA C 213 PASS ##> 5 T TC 242 PASS ##> ... ... ... ... ... ##> 402 G GT 488 PASS ##> 403 A AG 516 PASS ##> 404 T TAAAATAAATAAATAAATA 821 PASS ##> 405 A ATGT 1598 PASS ##> 406 A AT 378 PASS # REF列が"A"のバリアントだけを取得する fixed(subset(vcf, ref(vcf) == "A")) ##> DataFrame with 1459 rows and 4 columns ##> REF ALT QUAL FILTER ##> <DNAStringSet> <DNAStringSetList> <numeric> <character> ##> 1 A G 100 PASS ##> 2 A G 100 PASS ##> 3 A C 100 PASS ##> 4 A G 100 PASS ##> 5 A G 100 PASS ##> ... ... ... ... ... ##> 1455 A C 100 PASS ##> 1456 A AT 378 PASS ##> 1457 A G 100 PASS ##> 1458 A G 100 PASS ##> 1459 A G 100 PASS # QUAL列が2000以上のバリアントだけを取得する fixed(subset(vcf, qual(vcf) >= 2000)) ##> DataFrame with 6 rows and 4 columns ##> REF ALT QUAL FILTER ##> <DNAStringSet> <DNAStringSetList> <numeric> <character> ##> 1 A AAAACAATACCCAC 2137 PASS ##> 2 AGTGT A 2188 PASS ##> 3 TTGTG T 2641 PASS ##> 4 A AGT 2110 PASS ##> 5 ACAGT A 2487 PASS ##> 6 G GGAGGCA 2088 PASS
可視化
それでは最後にデータの可視化を少しだけしてみようと思います。 その前にデータをふるいにかけておこうと思います（このフィルタリングは特に意味ありません）。
vcf <- subset(vcf, info(vcf)$LDAF > 0.5) dim(vcf) ##> [1] 385 5
今回は遺伝子型をヒートマップ形式で表示してみたいと思います。まずはデータの取得と整形です。 表形式（Matrix型）のデータから
ggplot2 に対応した形式にするのには reshape2::melt が便利です。
# Genotype の取得 gt <- geno(vcf)$GT # 行および列がそれぞれ何を表しているかを指定しておく names(dimnames(gt)) <- c("marker", "sample") # 行列形式のデータを行名、列名、セルの値をそれぞれ新しい列に持つデータフレームに変換 gt <- reshape2::melt(gt, value.name = "genotype") # どんなデータ形式に変換されたか確認 head(gt) ##> marker sample genotype ##> 1 rs28523661 HG00096 0|0 ##> 2 rs28479081 HG00096 0|0 ##> 3 rs28695790 HG00096 1|1 ##> 4 rs28681372 HG00096 1|1 ##> 5 rs28465520 HG00096 1|1 ##> 6 rs4078443 HG00096 1|1
ではプロットしてみます。
ggplot(gt, aes(sample, marker, fill = genotype)) + geom_raster() + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
# A, B, H の形式でプロット code_genotype <- c("0|0" = "A", "0|1" = "H", "1|0" = "H", "1|1" = "B") gt <- dplyr::mutate(gt, genotype = code_genotype[genotype]) ggplot(gt, aes(sample, marker, fill = genotype)) + geom_raster() + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
おまけ
VCFの遺伝子型はREF型を 0、ALT型を 1~ の数字で表していますが、それを実際の塩基で表すには
genotypeCodesToNucleotides 関数を用います。
gt <- geno(genotypeCodesToNucleotides(vcf))$GT head(gt) ##> HG00096 HG00097 HG00099 HG00100 HG00101 ##> rs28523661 "T|T" "C|T" "T|C" "C|T" "T|C" ##> rs28479081 "T|T" "C|T" "T|C" "T|C" "C|T" ##> rs28695790 "T|T" "C|T" "T|T" "T|T" "T|T" ##> rs28681372 "A|A" "G|A" "G|G" "A|A" "A|A" ##> rs28465520 "G|G" "C|G" "C|C" "G|G" "G|G" ##> rs4078443 "C|C" "T|C" "T|T" "C|C" "C|C"
