How do I get the strand of a read when using htslib?

Question

[Click downvote]

If I want to read a bam like this

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <htslib/sam.h>

int main(int argc, char *argv[]){

    samFile *fp_in = hts_open(argv[1],"r"); //open bam file
    bam_hdr_t *bamHdr = sam_hdr_read(fp_in); //read header
    bam1_t *aln = bam_init1(); //initialize an alignment

    char *chrom = argv[2];
    int locus = atoi(argv[3]);
    int comp ;

    printf("%s\t%d\n", chrom, locus);

    //header parse
    //uint32_t *tar = bamHdr->text ;
    //uint32_t *tarlen = bamHdr->target_len ;

    //printf("%d\n",tar);

    while(sam_read1(fp_in,bamHdr,aln) > 0){

        int32_t pos = aln->core.pos +1; //left most position of alignment in zero based coordianate (+1)
        char *chr = bamHdr->target_name[aln->core.tid] ; //contig name (chromosome)
        uint32_t len = aln->core.l_qseq; //length of the read.

        uint8_t *q = bam_get_seq(aln); //quality string
        uint32_t q2 = aln->core.qual ; //mapping quality


        char *qseq = (char *)malloc(len);

        for(int i=0; i< len ; i++){
            qseq[i] = seq_nt16_str[bam_seqi(q,i)]; //gets nucleotide id and converts them into IUPAC id.
        }

        //printf("%s\t%d\t%d\t%s\t%s\t%d\n",chr,pos,len,qseq,q,q2);

         if(strcmp(chrom, chr) == 0){

             if(locus > pos+len){
                 printf("%s\t%d\t%d\t%s\t%s\t%d\n",chr,pos,len,qseq,q,q2);
             }
         }
    }

    bam_destroy1(aln);
    sam_close(fp_in);

    return 0;
}

How do I get the strand info? In sam files it seems to be represented as 16 for reverse, 0 for positive.

Answer 1

[Pierre Lindenbaum]

in sam.h:

/*! @abstract the read is mapped to the reverse strand */
#define BAM_FREVERSE      16

usage:

(aln->core.flag & BAM_FREVERSE)

[Click downvote]

So if the above is true, the alignment is reverse. Thanks :)

Answer 2

[Alex Reynolds]

I think you could call the following for each read:

bool is_reverse = bam_is_rev(aln);

See: https://github.com/samtools/htslib/blob/6e86e386bed5c80c5222bd2e8eb2fd8c99234909/htslib/sam.h#L200-L205

It's a convenience macro for the call described in Pierre's answer. You might need to add:

#include <stdbool.h>

to add support for the bool type.