目录
提取矩阵数据
一般的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数据读取自由，关键是可以直接写入内存进行下一步的统计分析和数据可视化，个人感觉是很有效的提高了生产力。值得花时间学习一下这个工具。
