将 Htslib 用于提取替代等位基因信息的 VCF 文件

Question

我正在使用 c++ 处理 VCF 文件，为此我使用来自 htslib (https://github.com/samtools/htslib/blob/develop/htslib/vcf.h) 的 vcf 库。我知道可能有一些更好的库，但我也在使用 htslib 也有库的其他文件格式，所以我想坚持使用 htslib。

我找到了一些代码示例，可以在文件中打开读取并创建正确的结构。标头并在此处使用 VCF 文件中的一些信息：https://gist.github.com/gatoravi/cad922bdf2b625a91126 和 http://wresch.github.io/2014/11/18/process-vcf-file-with-htslib.html

但是，如果我们坚持第一个示例，我已经将代码“解码”为以下代码，并将我的 cmets 解码为代码：

int main(int argc,char **argv){
  std::cerr << "Usage:subset.vcf " << std::endl;
  
  // htslib internally represents VCF as bcf1_t data structures

  htsFile *test_vcf = NULL;

  // creates header
  bcf_hdr_t *test_header = NULL;
  
  // initialize and allocate bcf1_t object
  bcf1_t *test_record = bcf_init();

  test_vcf = vcf_open("subset.vcf", "r");

  // returning a bcf_hdr_t struct 
  test_header = bcf_hdr_read(test_vcf);
  if(test_header == NULL) {throw std::runtime_error("Unable to read header.");}
  
  while(bcf_read(test_vcf, test_header, test_record) == 0){
    // std::cout << "pos " << test_record->pos << std::endl; //column 2 in VCF with the coordinate looks like its -1
    // std::cout << "length " << test_record->rlen << std::endl; // I assume its the length of the ALT
    // std::cout << "chrom " << test_record->rid; (-1) format or bcf_hdr_id2name(test_header, test_record->rid)
    // std::cout << "qual " << test_record->qual; //column 6
    // std::cout << "allele " << test_record->n_allele << std::endl; // number of alleles
    // std::cout << "info " << test_record->n_info << std::endl; // I dont know what this is
    // std::cout << "nfmt " << test_record->n_fmt << std::endl;
    // "sample " << test_record->n_sample // i dont know what this is
    std::cout << "chr" << bcf_hdr_id2name(test_header, test_record->rid) << ":" <<test_record->pos+1 << std::endl;

    std::cout << "------------------" << std::endl;
  }
  bcf_hdr_destroy(test_header);
  bcf_destroy(test_record); 
  bcf_close(test_vcf);
  return 0;
}

在上面的这段代码中，我在 while 循环中注释掉了多个 std::cout，以便让我的 cmets 清楚地了解其中的一些功能是什么——即“消除”是染色体。据我所知，vcf 库的名称“rid”或“nfmt”都是预定义的。运行此代码，我可以打印多个内容，例如染色体名称、位置等。但我有几个问题：

我的 VCF 文件具有 #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 的一般结构，其中几行仅显示前 6 列的小示例：

14  19058352    rs144287685 A   G   100
14  19066089    rs374285188 C   A,T 100
14  19075627    .   G   A,T 100
14  19075912    .   A   C,T 100
14  19237205    .   T   TATGTTATG   100

我的问题是在使用库时我希望打印出参考（第 4 列）和替代（第 5 列），所以对于第 1 行：REF = A & ALT = G，对于第 5 行：REF = T & ALT = TATGTTATG。

谁能帮助我准确理解提取这两个字段需要做什么？我在库描述中看不到如何使用“test_record->”来提取这些？

我希望我的问题有点道理。感谢您的时间和帮助。

Answer 1

我知道自从您 2 个月前发布此内容以来有点晚了，但也许这可以帮助其他人。我最近一直在努力使用 htslib，但我设法获得了 ALT 和 REF 值。

ALT 和 REF 的值存储在 bcf1_t 结构中，在名为 d 的字段中：

typedef struct bcf1_t {
    hts_pos_t pos;
    hts_pos_t rlen;
    int32_t rid;
    float qual;
    uint32_t n_info:16, n_allele:16;
    uint32_t n_fmt:8, n_sample:24;
    kstring_t shared, indiv;
    bcf_dec_t d;   //<----- HERE
    int max_unpack;
    int unpacked;
    int unpack_size[3];
    int errcode;
} bcf1_t;

当你初始化一个bcf1_t对象时，该字段默认是不填充的，所以首先你必须调用函数bcf_unpack。该函数的第一个参数是指向记录的指针，第二个参数取决于您希望它解包的值。在您的情况下，对于 ALT 和 REF，我认为第一个参数应该是 test_record，第二个参数应该是 BCF_UN_STR。在 htslib 的源代码中，所有可用的值都被注释掉了。

int bcf_unpack(bcf1_t *b, int which);

现在您可以查看 d 字段。该字段的类型是另一个称为 bcf_dec_t 的结构。在这里你必须看看字段als。

typedef struct bcf_dec_t {
    int m_fmt, m_info, m_id, m_als, m_allele, m_flt;
    int n_flt;
    int *flt;
    char *id, *als;     // ID and REF+ALT block (\0-separated)
    char **allele;
    bcf_info_t *info;
    bcf_fmt_t *fmt;
    bcf_variant_t *var;called
    int n_var, var_type;
    int shared_dirty;
    int indiv_dirty;
} bcf_dec_t;

正如文档中所说，als 包含由“\0”分隔的值 REF 和 ALT。因此，如果您的值是：REF = T 和 ALT = TATGTTATG，als 包含以下字符数组：“T\0TATGTTATG\0”。

您可以解析该字符数组以获取 REF 和 ALT。我使用我编写的以 als 作为输入的函数来执行此操作，并返回一个包含分隔的 ALT 和 REF 的向量。我知道这可能不是最优化的功能，并且必须有一种使用 htslib 的方法，但它确实有效：

std::vector<std::string> extractAltRef(char *als) {
    std::vector<std::string> res;
    std::string str = "";
    int i = 0;
    int j = 0;
    while(j != 2) {
        if(als[i] == '\0') {
            res.push_back(str);
            str.clear();
            j++;
        } else {
            str += als[i];
        }
        i++;
    }
    return res;
}

所以在你的代码中，为了获取 REF 和 ALT 值，如果你使用我的函数，你应该做这样的事情（未经测试）：

// Pointer initializations...
bcf_unpack(test_record, BCF_UN_STR);
while(bcf_read(test_vcf, test_header, test_record) == 0){
    std::vector<std::string> altRef = extractAltRef(test_record->d.als);
    std::cout << "REF " << altRef[0] << std::endl;
    std::cout << "ALT " << altRef[1] << std::endl;
}
// Free memory...

如果你不使用我的代码，你将不得不找到一种分离 REF 和 ALT 的方法。

希望对你有帮助，

阿尔贝托。