Как я могу посчитать количество людей в группах населения, как указано в порядке, из файла VCF - PullRequest
Новая
       13

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

0 голосов
/ 02 мая 2018

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

Например, в файле, с которым я связан, у меня есть две группы пользователей: Chalifour 2003 и Chalifour 2015. У отдельных лиц есть префикс "CHALIFOUR_2003 ...", который определяет это.

Я бы хотел иметь возможность извлечь что-то вроде: Chalifour_2003 * 35 Чалифур 2015 * 45

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

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

С уважением, Элла

1 Ответ

0 голосов
/ 02 мая 2018

Используя пакет 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
...