Блог

Я работаю с c для обработки файлов VCF, для этого я использую библиотеку vcf из htslib (https://github.com/samtools/htslib/blob/develop/htslib/vcf.h ). Я знаю, что могут быть библиотеки получше, но я также работаю с другими форматами файлов, для которых htslib также имеет библиотеки, поэтому я хотел бы придерживаться htslib.

Я нашел несколько примеров кода для открытия чтения в файле и создания правильной структуры. заголовок и использование некоторой информации из файла VCF здесь: https://gist.github.com/gatoravi/cad922bdf2b625a91126 и http://wresch.github.io/2014/11/18/process-vcf-file-with-htslib.html

Но если мы будем придерживаться первого примера, я «расшифровал» код следующим образом с моими комментариями к коду:

 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, чтобы в моих комментариях было ясно, каковы некоторые функциональные возможности, т.е. «rid» — это хромосома. Насколько я могу прочитать для библиотеки vcf, все имена «rid» или «nfmt» предварительно определены. Запустив этот код, я могу печатать несколько вещей, таких как имена хромосом, положение среди прочего. Но у меня есть несколько проблем:

Мой файл VCF имеет общую структуру в ФОРМАТЕ #CHROM POS ID REF ALT ИНФОРМАЦИЯ о КАЧЕСТВЕННОМ ФИЛЬТРЕ, с небольшим примером нескольких строк, показывающих только первые 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 amp; ALT = G, для строки 5: REF = T amp; ALT = TATGTTATG.

Может ли кто-нибудь помочь мне понять, что именно мне нужно сделать, чтобы извлечь эти два поля? Я не вижу в описании библиотеки, как использовать «test_record->» для их извлечения?

Я надеюсь, что мой вопрос имеет некоторый смысл. Спасибо за ваше время и помощь.

Я знаю, что немного поздно, так как вы опубликовали это 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;
 

 int bcf_unpack(bcf1_t *b, int which);
 

 typedef struct bcf_dec_t {
    int m_fmt, m_info, m_id, m_als, m_allele, m_f<
    int n_f<
    int *f<
    char *id, *als;     // ID and REF ALT block (-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;
 

 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] == '') {
            res.push_back(str);
            str.clear();
            j  ;
        } else {
            str  = als[i];
        }
        i  ;
    }
    return res;
}
 

 // 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...