0
$\begingroup$

I am creating a threaded c++ file where i generate in silico bam files, using header, DNA sequence and read information.

First i use bam_init1() to create the bam1_t structure just named "b".

Then i use bam_set1 to create the actual sequence entry in the bam file

bam_set1(b,read_id_length,READ_ID,flag,chr_idx,min_beg,mapq,n_cigar,cigar,-1,-1,0,strlen(DNAsequence),DNAsequence,quality_string,l_aux)

And finally after bam_set1 I use sam_write1 to save the bam sequence entry to the samfile.

samFile *SAMout = NULL;
sam_hdr_t *SAMHeader;
htsFormat *fmt_hts =(htsFormat*) calloc(1,sizeof(htsFormat));
sam_write1(SAMout,SAMHeader,b)

The issue I have is that I want to add some auxiliary tags, such as the NM or MD tag.

I just want to create an MD tag that is just equal to the sequence length. But when i use the bam_aux_update_str function I get the error message "[E::bam_aux_get] Corrupted aux data for read (null)" when compiling and then my bam file is truncated.

I have tried using this function before I use bam_set1

const char* MDtag = "\tMD:Z:";
bam_aux_update_str(b,"MD",strlen(DNAsequence),MDtag);

bam_set1(b,read_id_length,READ_ID,flag,chr_idx,min_beg,mapq,n_cigar,cigar,-1,-1,0,strlen(DNAsequence),DNAsequence,quality_string,l_aux);

The info regarding the bam_aux_update_str function can be found here: https://docs.rs/hts-sys/1.11.1-fix1/hts_sys/fn.bam_aux_update_str.html

My example is similar to other examples I have found online, they create other tags, but using the same bam_aux_update_str function: https://www.itcodet.com/cpp/cpp-bam_aux2a-function-examples.html

The only difference is this link doesn't use bam_set1 or sam_write since it's opening an existing sam file, but i need these functions to create in silico sequence reads in sam/bam format.

https://github.com/samtools/samtools/blob/develop/bam_ampliconclip.c It uses sam_write.

Anyone have any suggestions to solve this specific issue adding NM and MD tags to the bam file using this function or any other function.?

$\endgroup$

0

Your Answer

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

Browse other questions tagged or ask your own question.