问题描述
我正在使用 htslib 提取包含在 C++ 中的 VCF 文件中的所有信息。
目前,多亏了VCF specification和文件vcf.h中的文档,我已经成功提取了header(Meta-information Lines)中的所有元数据信息,以及其中包含的大部分信息文件正文的每一行(数据行)。
但是,我不知道如何提取基因型信息(样本列)。
我正在使用 1000G 项目中的示例文件。这是文件两行的示例,它显示了格式字段和两个样本(该文件每行有 1000 多个样本,我想提取所有样本的数据):
FORMAT HG00096 HG00097
GT:DS:GL 0|0:0.050:-0.48,-0.48,-0.48 0|0:0.050:-0.24,-0.40,-1.49
GT:DS:GL 0|0:0.000:-0.10,-0.69,-4.70 0|0:0.000:-0.05,-0.94,-5.00
我知道这是一项繁重的任务,需要一些计算时间。我已经提取了每列的名称(HG00096,HG00077 ...),但我不知道如何将每个样本的信息提取为完整字符串(例如,“0|0:0.050:-0.48,- 0.48,-0.48"),作为键值对的集合(数组、映射、向量...)(例如,[("GT","0|0"),("DS","0.050"),("GL","-0.48,-0.48")),或简单地作为值数组(例如,["0|0","0.050",-0.48" ]. 我想对每个样本都这样做。
我一直在阅读 vcf.h 文件中的文档,我认为函数 bcf_get_genotypes(hdr,line,dst,ndst) 可能适用于此,但我不确定如何使用它用于将值提取为字符串。另外,我认为此信息可能存储在“bcf_fmt_t”的“p”指针内,但我不确定,它只包含一组 uint8_t 值,我不知道字符串(或char数组)可以按照我想要的方式提取。
typedef struct bcf_fmt_t {
int id;
int n,size,type;
uint8_t *p;
uint32_t p_len;
uint32_t p_off:31,p_free:1;
} bcf_fmt_t;
有没有我正在尝试的方法?
解决方法
我终于明白了。有一些函数可以执行此操作,具体取决于格式 ID 标头中指定的类型:这些函数位于 htslib 中的 vcf.h 文件中:
#define bcf_get_format_int32(hdr,line,tag,dst,ndst)
bcf_get_format_values(hdr,(void**)(dst),ndst,BCF_HT_INT)
#define bcf_get_format_float(hdr,BCF_HT_REAL)
#define bcf_get_format_char(hdr,ndst)
bcf_get_format_values(hdr,BCF_HT_STR)
#define bcf_get_genotypes(hdr,ndst)
bcf_get_format_values(hdr,"GT",BCF_HT_INT)
HTSLIB_EXPORT
int bcf_get_format_string(const bcf_hdr_t *hdr,bcf1_t *line,const char *tag,char ***dst,int *ndst);
HTSLIB_EXPORT
int bcf_get_format_values(const bcf_hdr_t *hdr,void **dst,int *ndst,int type);