Skip to content

Instantly share code, notes, and snippets.

@lindenb
Created December 10, 2010 10:25
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save lindenb/736059 to your computer and use it in GitHub Desktop.
Save lindenb/736059 to your computer and use it in GitHub Desktop.
/*
gcc -m64 -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 <assert.h>
#include "bam.h"
#include "sam.h"
int main(int argc, char *argv[]) {
samfile_t *fp_in = NULL;
bam1_t *b=NULL;
if(argc!=2) return EXIT_FAILURE;
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;
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("N"); break;
case BAM_CSOFT_CLIP: printf("S");break;
case BAM_CHARD_CLIP: printf("R");break;
case BAM_CPAD: printf("P");break;
default:printf("?");assert(0);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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment