User:Lindenb/Notebook/UMR915/20101209
From OpenWetWare
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; }