0

我想从 vcf 文件中获取每个种群中的个体数量，按照种群的读取顺序。我文件的字段如下所示

##fileformat=VCFv4.2                                                
##fileDate=20180425                                             
##source="Stacks v1.45"                                             
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">                                              
##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">                                               
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">                                                
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">                                             
##FORMAT=<ID=AD,Number=1,Type=Integer,Description="Allele Depth">                                               
##FORMAT=<ID=GL,Number=.,Type=Float,Description="Genotype Likelihood">                                              
##INFO=<ID=locori,Number=1,Type=Character,Description="Orientation the 
corresponding Stacks locus aligns in">                                              
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT   
CHALIFOUR_2003_ChHis-1  CHALIFOUR_2003_ChHis-13 CHALIFOUR_2003_ChHis-14  
CHALIFOUR_2003_ChHis-15
un  1027    13_65   C   T   .   PASS    NS=69;AF=0.188;locori=p GT:DP:AD     
0/1:16:9,7  0/0:39:39,0 0/0:17:17,0 0/0:39:39,0

在此处查看示例文件vcf 文件

例如，在我链接到的文件中，我有两个群体，Chalifour 2003 和 Chalifour 2015。个人有一个前缀“CHALIFOUR_2003 ...”来标识这一点。

我希望能够提取如下内容： Chalifour_2003* 35 Chalifour 2015* 45

用“35”和“45”表示每个种群中的个体数量（尽管这些数字是虚构的）。我根本不关心输出的格式，我只需要数字，重要的是要按照将人口读入文件的顺序列出人口。

任何有关尝试获取此信息的途径的建议将不胜感激。

4

1 回答 1

0

使用data.table包读取 vcf 文件，您可以执行以下操作：

library(data.table)
df <- fread("~/Downloads/ChaliNoOddsWithOuts.vcf")
samples <- colnames(df)[-c(1:9)]
table(gsub("(.*_.*)_.*","\\1", samples))

如果您不坚持使用，R那么这是bash可以完成工作的一种衬里

grep "#CHROM" file.vcf | tr "\t" "\n " | tail -n +10 | cut -f1,2 -d'_' | uniq -c
于 2018-05-02T16:39:34.347 回答