用R语言对vcf文件进行数据挖掘.3 从vcf文件里提取有用信息

目录

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

提取矩阵数据

一般的VCF文件都很大，用手动提取里面的信息肯定不大现实。用vcfR就可以轻松实现。
vcfR自带测试文件vcfR_test。就用这个文件来操作一下吧。

library(vcfR)
data("vcfR_test")
head(vcfR_test)
[1] "***** Object of class 'vcfR' *****"
[1] "***** Meta section *****"
[1] "##fileformat=VCFv4.3"
[1] "##fileDate=20090805"
[1] "##source=myImputationProgramV3.1"
[1] "##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta"
[1] "##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d [Truncated]"
[1] "##phasing=partial"
[1] "First 6 rows."
[1] 
[1] "***** Fixed section *****"
     CHROM POS       ID          REF   ALT      QUAL FILTER
[1,] "20"  "14370"   "rs6054257" "G"   "A"      "29" "PASS"
[2,] "20"  "17330"   NA          "T"   "A"      "3"  "q10" 
[3,] "20"  "1110696" "rs6040355" "A"   "G,T"    "67" "PASS"
[4,] "20"  "1230237" NA          "T"   NA       "47" "PASS"
[5,] "20"  "1234567" "microsat1" "GTC" "G,GTCT" "50" "PASS"
[1] 
[1] "***** Genotype section *****"
     FORMAT        NA00001          NA00002          NA00003       
[1,] "GT:GQ:DP:HQ" "0|0:48:1:51,51" "1|0:48:8:51,51" "1/1:43:5:.,."
[2,] "GT:GQ:DP:HQ" "0|0:49:3:58,50" "0|1:3:5:65,3"   "0/0:41:3"    
[3,] "GT:GQ:DP:HQ" "1|2:21:6:23,27" "2|1:2:0:18,2"   "2/2:35:4"    
[4,] "GT:GQ:DP:HQ" "0|0:54:7:56,60" "0|0:48:4:51,51" "0/0:61:2"    
[5,] "GT:GQ:DP"    "0/1:35:4"       "0/2:17:2"       "1/1:40:3"    
[1] 
[1] "Unique GT formats:"
[1] "GT:GQ:DP:HQ" "GT:GQ:DP"   
[1] 

在分区Genotype里，通过观察FORMAT列可以看到一共有四种类型的数据GT:GQ:DP:HQ，至于这四种类型的数据个各自代表什么意思大家可以查阅知乎百度谷歌。我们可以提取出我们想要的数据类型。比方说最重要的GT(genotype)。

gt <- extract.gt(vcfR_test)
gt
           NA00001 NA00002 NA00003
rs6054257  "0|0"   "1|0"   "1/1"  
20_17330   "0|0"   "0|1"   "0/0"  
rs6040355  "1|2"   "2|1"   "2/2"  
20_1230237 "0|0"   "0|0"   "0/0"  
microsat1  "0/1"   "0/2"   "1/1" 

同样，我们也可以提取例如DP(测序深度Read Depth)的数字矩阵。

gt <- extract.gt(vcfR_test, element = 'DP', as.numeric = TRUE)
gt
           NA00001 NA00002 NA00003
rs6054257        1       8       5
20_17330         3       5       3
rs6040355        6       0       4
20_1230237       7       4       2
microsat1        4       2       3

值的注意的是这里用到了参数as.numeric = TRUE使得数据自动转换成了数字。但是并不是对所有类型的数据都有效，比方说我们重复一下提取gt。

> gt <- extract.gt(vcfR_test, element = 'GT', as.numeric = TRUE)
> gt
           NA00001 NA00002 NA00003
rs6054257        0       1       1
20_17330         0       0       0
rs6040355        1       2       2
20_1230237       0       0       0
microsat1        0       0       1

在没有任何报错的情况下gt变成了一堆毫无意义的数字，很明显不合理，不要用这些经过错误转换的数据进行下一步分析，比方说喜闻乐见的主成分分析。

数据拆分

在一些类型的数据里可能会出现一个以上的结果，比方说上面的HQ数据。

> gt <- extract.gt(vcfR_test, element = 'HQ')
> gt
           NA00001 NA00002 NA00003
rs6054257  "51,51" "51,51" ".,."  
20_17330   "58,50" "65,3"  NA     
rs6040355  "23,27" "18,2"  NA     
20_1230237 "56,60" "51,51" NA     
microsat1  NA      NA      NA  

一般情况下我们只需要每一列的第一个数字

> myHQ1 <- masplit(gt[,1:2], sort = 0)
> myHQ1
           NA00001 NA00002
rs6054257       51      51
20_17330        58      65
rs6040355       23      18
20_1230237      56      51
microsat1       NA      NA

不需要samtools之类的软件我们也可以实现vcf数据读取自由，关键是可以直接写入内存进行下一步的统计分析和数据可视化，个人感觉是很有效的提高了生产力。值得花时间学习一下这个工具。