目錄
提取矩陣數(shù)據(jù)
一般的VCF文件都很大歼跟，用手動提取里面的信息肯定不大現(xiàn)實断凶。用vcfR就可以輕松實現(xiàn)。
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]
在分區(qū)Genotype里，通過觀察FORMAT列可以看到一共有四種類型的數(shù)據(jù)GT:GQ:DP:HQ，至于這四種類型的數(shù)據(jù)個各自代表什么意思大家可以查閱知乎百度谷歌。我們可以提取出我們想要的數(shù)據(jù)類型荧止。比方說最重要的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)的數(shù)字矩陣跃巡。
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
值的注意的是這里用到了參數(shù)as.numeric = TRUE使得數(shù)據(jù)自動轉(zhuǎn)換成了數(shù)字。但是并不是對所有類型的數(shù)據(jù)都有效牧愁，比方說我們重復(fù)一下提取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變成了一堆毫無意義的數(shù)字，很明顯不合理猪半，不要用這些經(jīng)過錯誤轉(zhuǎn)換的數(shù)據(jù)進(jìn)行下一步分析兔朦，比方說喜聞樂見的主成分分析。
數(shù)據(jù)拆分
在一些類型的數(shù)據(jù)里可能會出現(xiàn)一個以上的結(jié)果磨确，比方說上面的HQ數(shù)據(jù)沽甥。
> 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
一般情況下我們只需要每一列的第一個數(shù)字
> 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之類的軟件我們也可以實現(xiàn)vcf數(shù)據(jù)讀取自由，關(guān)鍵是可以直接寫入內(nèi)存進(jìn)行下一步的統(tǒng)計分析和數(shù)據(jù)可視化乏奥，個人感覺是很有效的提高了生產(chǎn)力摆舟。值得花時間學(xué)習(xí)一下這個工具。
