目录
相信经常使用R的同学们对tidy格式的数据并不陌生。接近标准格式的数据框，非常便于操作。其实vcf数据也可以通过vcfR转变成tidy格式的数据。这次我们会继续使用vcfR自带测试文件vcfR_test来教学。
library(vcfR)
data("vcfR_test")
vcfR_test
***** Object of Class vcfR *****
3 samples
1 CHROMs
5 variants
Object size: 0 Mb
0 percent missing data
***** ***** *****
函数vcfR2tidy()会将这个数据变成tibble形式的tidy数据。在此之前我们可以通过vcf_field_names()函数来查看这个vcf里包含着哪些类型的数据。比方说查看一下FORMAT，结果显示FORMAT里有四种类型GT,GQ,DP,HQ，各自包含几个数据，分别代表什么意思等等。
vcf_field_names(vcfR_test, tag = "FORMAT")
> vcf_field_names(vcfR_test, tag = "FORMAT")
# A tibble: 4 x 5
Tag ID Number Type Description
<chr> <chr> <chr> <chr> <chr>
1 FORMAT GT 1 String Genotype
2 FORMAT GQ 1 Integer Genotype Quality
3 FORMAT DP 1 Integer Read Depth
4 FORMAT HQ 2 Integer Haplotype Quality
提取GT,DP并转变数据。形成一个list。
> Z <- vcfR2tidy(vcfR_test, format_fields = c("GT", "DP"))
Extracting gt element GT
Extracting gt element DP
查看一下刚才转变出来的数据Z里面有点什么。Z里有fix, gt, meta 三组数据。
> names(Z)
[1] "fix" "gt" "meta"
再分别看一下吧。
Z$meta
> Z$meta
# A tibble: 8 x 5
Tag ID Number Type Description
<chr> <chr> <chr> <chr> <chr>
1 INFO NS 1 Integer Number of Samples With Data
2 INFO DP 1 Integer Total Depth
3 INFO AF A Float Allele Frequency
4 INFO AA 1 String Ancestral Allele
5 INFO DB 0 Flag dbSNP membership, build 129
6 INFO H2 0 Flag HapMap2 membership
7 FORMAT gt_GT 1 String Genotype
8 FORMAT gt_DP 1 Integer Read Depth
> Z$gt
# A tibble: 15 x 6
ChromKey POS Indiv gt_GT gt_DP gt_GT_alleles
<int> <int> <chr> <chr> <int> <chr>
1 1 14370 NA00001 0|0 1 G|G
2 1 17330 NA00001 0|0 3 T|T
3 1 1110696 NA00001 1|2 6 G|T
4 1 1230237 NA00001 0|0 7 T|T
5 1 1234567 NA00001 0/1 4 GTC/G
6 1 14370 NA00002 1|0 8 A|G
7 1 17330 NA00002 0|1 5 T|A
8 1 1110696 NA00002 2|1 0 T|G
9 1 1230237 NA00002 0|0 4 T|T
10 1 1234567 NA00002 0/2 2 GTC/GTCT
11 1 14370 NA00003 1/1 5 A/A
12 1 17330 NA00003 0/0 3 T/T
13 1 1110696 NA00003 2/2 4 T/T
14 1 1230237 NA00003 0/0 2 T/T
15 1 1234567 NA00003 1/1 3 G/G
> Z$fix
# A tibble: 5 x 14
ChromKey CHROM POS ID REF ALT QUAL FILTER NS DP AF AA DB
<int> <chr> <int> <chr> <chr> <chr> <dbl> <chr> <int> <int> <chr> <chr> <lgl>
1 1 20 1.44e4 rs60… G A 29 PASS 3 14 0.5 NA TRUE
2 1 20 1.73e4 NA T A 3 q10 3 11 0.017 NA FALSE
3 1 20 1.11e6 rs60… A G,T 67 PASS 2 10 0.33… T TRUE
4 1 20 1.23e6 NA T NA 47 PASS 3 13 NA T FALSE
5 1 20 1.23e6 micr… GTC G,GT… 50 PASS 3 9 NA G FALSE
# … with 1 more variable: H2 <lgl>
这些数据看上去应该很眼熟了吧，可以直接用tidyverse包来操作。至于tidyverse怎么用可以参照我的文集。
