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