User:Lindenb/Notebook/UMR915/20101209

20101208        Top        20101210       

SAM C API

read a BAM file, echo the data:

/*
gcc -Wall -I package/sam/samtools-0.1.7a -L package/sam/samtools-0.1.7a jeter.c -lbam -lz
*/
#include <stdlib.h>
#include <stdio.h>
   
#include "bam.h"
#include "sam.h"
 
int main(int argc, char *argv[]) {
 
  samfile_t *fp_in = NULL;
  bam1_t *b=NULL;

  fp_in = samopen(argv[1], "rb", 0);

  if(NULL == fp_in)
	{
       fprintf(stderr,"Could not open file\n");
       return EXIT_FAILURE;
       }
 
  b = bam_init1();
  int pos=0;
  int lastpos=0;

  while(samread(fp_in, b) > 0)
    {
    lastpos = pos;
    pos = b->core.pos;
         
    if(pos != lastpos) {
      printf("tid  : %d\n",b->core.tid);
     if(b->core.tid>=0)
	{
      	printf("seqname: %s\n", fp_in->header->target_name[b->core.tid]);
	printf("seqlength: %d\n", fp_in->header->target_len[b->core.tid]);
	}
      printf("pos  : %d\n",b->core.pos);
      char *name  = bam1_qname(b);
      unsigned char *qual  = bam1_qual(b);
      uint32_t* cigar= bam1_cigar(b);
      int n=0;
      char *qseq = (char *) malloc(b->core.l_qseq+1);
      unsigned char *s   = bam1_seq(b);
      for(n=0;n<(b->core.l_qseq);n++) {
        int v = bam1_seqi(s,n);
        qseq[n] = bam_nt16_rev_table[v];
      }
      qseq[n] = 0;

      int cigar_leng=bam_cigar2qlen(&(b->core),cigar);
      printf("name : %s\n",name);
      printf("qseq : %s\n",qseq);
      printf("cigar-size: %d\n",cigar_leng);
       printf("cigar:");
      int k;
      for( k=0;k< b->core.n_cigar;++k)
	{
        int cop =cigar[k] & BAM_CIGAR_MASK; // operation
	int cl = cigar[k] >> BAM_CIGAR_SHIFT; // length
	switch(cop)
		{
		case BAM_CMATCH: printf("M");break;
		case BAM_CINS: printf("I");break;
		case BAM_CDEL: printf("D");break;
		case BAM_CREF_SKIP: printf("S"); break;//TODO
		case BAM_CSOFT_CLIP: printf("A");break;//TODO
		case BAM_CHARD_CLIP: printf("R");break;//TODO
		case BAM_CPAD: printf("P");break;
		default:printf("?");break;
		}
	printf("%d",cl);
	}
      printf("\n");

      printf("qual :");
      for(n=0;n<(b->core.l_qseq);n++) {
        printf(" %d",qual[n]);
      }
      printf("\n\n");
    }
    bam_destroy1(b);
    b = bam_init1();
  }
  bam_destroy1(b);

  samclose(fp_in);
  return 0;
}