Используйте Htslib для файлов VCF, извлекая информацию об альтернативных аллелях

Question

Я работаю с 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 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 содержит значения REF и ALT, разделенные '\ 0'. Итак, если ваши значения: 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.

Надеюсь, это поможет,

Альберто.