1

I'm using the htslib library for reading SAM/BAM files, it works perfectly. I can also write the alignments back to a new SAM/BAM file.

For example, the following code prints the DNA sequence of an alignment:

    bam1_t *b = ...;
    int i;
    for (i = 0; i < b->core.l_qseq; ++i) {
        printf("%c", seq_nt16_str[bam_seqi(bam_get_seq(b),i)]);
    }

Question: How do I change the query sequence? Say, change the first letter to 'T'? bam_get_seq returns the sequence of a read, but there is no bam_set_seq function? Ideally, I'm looking for something like:

 bam_set_seq(b, 'TTTT') # My new DNA sequence

If I can figure out how to do the update, I know how to write the information to a new SAM/BAM file.

Timur Shtatland
  • 12,024
  • 2
  • 30
  • 47
ABCD
  • 7,914
  • 9
  • 54
  • 90
  • Why not just use Samtools and linux to create a text file of your alignment and then in perl/python manipulate the sequence as a string? – Seekheart Jun 10 '16 at 12:03
  • Samtools seems to be focused on collecting stats/viewing existing data in various ways, etc., but not on facilitating arbitrary changes to underlying data? – flyingfinger Jun 16 '16 at 20:05
  • @flyingfinger I don't think it's possible. – ABCD Jun 16 '16 at 22:12

0 Answers0