用R语言对vcf文件进行数据挖掘.4 tidy vcfR

目录

1. 前言
2. 方法简介
3. 从vcf文件里提取有用信息
4. tidy vcfR
5. vcf可视化1
6. vcf可视化2
7. 测序深度覆盖度
8. 窗口缩放
9. 如何单独分离染色体
10. 利用vcf信息判断物种染色体倍数
11. CNV分析

相信经常使用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怎么用可以参照我的文集。