Как я могу подсчитать количество особей в популяциях, перечисленных по порядку, из файла vcf

Я хотел бы получить количество особей в каждой популяции в порядке чтения популяций из файла 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...», который идентифицирует это.

Я хотел бы иметь возможность извлечь что-то вроде: Chalifor_2003* 35 Chaifour 2015* 45

При этом «35» и «45» обозначают количество особей в каждой популяции (хотя эти числа выдуманы). Меня совершенно не волнует формат вывода, мне нужны только числа, и важно, чтобы популяции были перечислены в том порядке, в котором они будут считываться в файл.

Любые предложения по способам получения этой информации будут высоко оценены.


r bioinformatics vcf-variant-call-format
person Ella Bowles    schedule 02.05.2018    source источник
comment
вы смотрели на пакет vcfR? У него есть read.vcfR, после чего вы можете просто извлечь префикс в новый столбец и использовать group_by() и count() из пакета dplyr.   -  person    schedule 02.05.2018

Ответы (1)


arrow_upward
0
arrow_downward

Используя пакет 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
person GordonShumway    schedule 02.05.2018
comment
@EllaBowles: если это рабочее решение, вам следует рассмотреть возможность предоставления ему флажок. - person AkselA; 25.05.2018