1
$\begingroup$

I would like to print the alternative allele from a vcf file which i load in with htslib. The description of some functions: https://github.com/samtools/htslib/blob/develop/htslib/vcf.h

I have used bcf1_t objects by creating a header and reading the information.

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("Testfile.vcf", "r");

  // returning a bcf_hdr_t struct 
  test_header = bcf_hdr_read(test_vcf);
  
  while(bcf_read(test_vcf, test_header, test_record) == 0){
    std::cout << "chr" << bcf_hdr_id2name(test_header, test_record->rid) << ":" <<test_record->pos+1 << std::endl;
  }
  bcf_hdr_destroy(test_header);
  bcf_destroy(test_record); 
  bcf_close(test_vcf);

But im not able to extract the reference or alternative allele information using this structure bcf1_t.

As far as i can read its perhaps possible with the structure and function:

bcf_dec_t
bcf_unpack()

But i am not able to get it to work since all the other functions seems limited to the bcf1_t struct.

Can anyone help with extracting the REF and ALT sequence? Thanks for your time

$\endgroup$
1
$\begingroup$

see the code in https://github.com/samtools/htslib/blob/develop/vcf.c#L3104

 bcf_unpack((bcf1_t*)v, BCF_UN_ALL);// if you need genotypes too, might use another constant to avoid decoding all the structure
 (...)
 kputc('\t', s); // ALT
 if (v->n_allele > 1) {
        for (i = 1; i < v->n_allele; ++i) {
            if (i > 1) kputc(',', s);
            kputs(v->d.allele[i], s);
        }
    }
$\endgroup$
4
  • $\begingroup$ Thanks for reply, i'll try that!. $\endgroup$ Nov 24 '20 at 12:36
  • $\begingroup$ Okay, so as far as i understand based on this function, s is a kstring_t* which you can create using mem "kstring_t s = (kstring_t)&test_header->mem;", and then all the information is added to the pointer "s". But how are we then able to print or access the information in "s" ? I can see its a pointer declared with a reference, but then i should be able to print the values with *s. Again thanks for your time. $\endgroup$ Nov 24 '20 at 13:20
  • $\begingroup$ kstring_t is just the implementation of a dynamic string. you could just use ` cout << v->d.allele[i] << endl; ` $\endgroup$
    – Pierre
    Nov 24 '20 at 14:51
  • $\begingroup$ okay, that was my initial try, i just looked into it, and the problem is my n_allele for some reason is 0. So ill play around with it some more. Thanks again. $\endgroup$ Nov 24 '20 at 14:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.