gclib/bamread/bam2fastx.cpp

#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <getopt.h>

#include "bam/bam.h"
#include "bam/sam.h"

bool is_fastq=true; //default is fastq
bool sam_input=false; //default is BAM
bool all_reads=false;
bool mapped_only=false;

#define USAGE "Usage: bam2fastx [--fasta|-a|--fastq|-q] [--sam|-s|-t] \n\
        [-M|--mapped-only|-A|--all] <in.bam>|<in.sam> \n"

char qseq[2048];

const char *short_options = "aqstMU";

enum {
   OPT_FASTA = 127,
   OPT_FASTQ,
   OPT_SAM,
   OPT_MAPPED_ONLY,
   OPT_ALL
   };
   
struct option long_options[] = {
  {"fasta", no_argument, 0, OPT_FASTA},
  {"fastq", no_argument, 0, OPT_FASTQ},
  {"sam", no_argument, 0, OPT_SAM},
  {"mapped-only", no_argument, 0, OPT_MAPPED_ONLY},
  {"all", no_argument, 0, OPT_ALL},
  {0, 0, 0, 0} // terminator
  };

int parse_options(int argc, char** argv)
{
  int option_index = 0;
  int next_option;
  do {
     next_option = getopt_long(argc, argv, short_options, long_options, &option_index);
     switch (next_option) {
       case -1:    
         break;
       case 'a':
       case OPT_FASTA:
         is_fastq = false;
         break;
       case 'q':
       case OPT_FASTQ:
         is_fastq = true;
         break;
       case 's':
       case 't':
       case OPT_SAM: //sam (text) input
         sam_input = true;
         break;
       case 'M':
       case OPT_MAPPED_ONLY:
         mapped_only = true;
         break;
       case 'A':
       case OPT_ALL:
         all_reads = true;
         break;
       default:
         return 1;
       }
   } while(next_option != -1);
  if (all_reads && mapped_only) {
    fprintf(stderr, "Error: incompatible options !\n");
    exit(2);
    }
  return 0;
}

#define bam_unmapped(b) (((b)->core.flag & BAM_FUNMAP) != 0)

void showfastq(const bam1_t *b, samfile_t* fp) {
  char *name  = bam1_qname(b);
  char *qual  = (char*)bam1_qual(b);
  char *s    = (char*)bam1_seq(b);
  int i;
  bool ismapped=((b->core.flag & BAM_FUNMAP) == 0);
  if (ismapped && !all_reads) return;
  if (mapped_only && !ismapped) return;
   
  for(i=0;i<(b->core.l_qseq);i++) {
    int8_t v = bam1_seqi(s,i);
    qseq[i] = bam_nt16_rev_table[v];
    }
  qseq[i] = 0;
  
  printf("@%s\n%s\n",name, qseq);
  for(i=0;i<(b->core.l_qseq);i++) {
    qseq[i]=qual[i]+33;
    }
  qseq[i]=0;
  printf("+\n%s\n",qseq);
}

void showfasta(const bam1_t *b, samfile_t* fp) {
  char *name  = bam1_qname(b);
  char *s    = (char*)bam1_seq(b);
  int i;
  for(i=0;i<(b->core.l_qseq);i++) {
    int8_t v = bam1_seqi(s,i);
    qseq[i] = bam_nt16_rev_table[v];
  }
  qseq[i] = 0;
  printf(">%s\n%s\n",name, qseq);
}


int main(int argc, char *argv[])
{
    samfile_t *fp;
    char* fname=NULL;

    if (parse_options(argc, argv) || optind>=argc) {
       fprintf(stderr, USAGE);
       return -1;
       }
    fname=argv[optind++];

    if (fname==NULL || fname[0]==0) {
        fprintf(stderr, USAGE);
        return 1;
        }
    if (sam_input)
        fp = samopen(fname, "r", 0);
      else 
        fp = samopen(fname, "rb", 0);
    if (fp == 0) {
        fprintf(stderr, "Error: bam2fastx failed to open BAM file %s\n", fname);
        return 1;
        }
    
    bam1_t *b = bam_init1();
    if (is_fastq) {
        while (samread(fp, b) >= 0) showfastq(b, fp);
        }
      else {
        while (samread(fp, b) >= 0) showfasta(b, fp);
        }
    
    bam_destroy1(b);
    
    samclose(fp);
    return 0;
}
Revision:	123
Committed:	Mon Dec 5 22:21:21 2011 UTC (7 years, 8 months ago) by gpertea
File size:	3419 byte(s)
Log Message:	bamread files
Line	File contents
1	#include <cstdlib>
2	#include <cstdio>
3	#include <cstring>
4	#include <getopt.h>
5
6	#include "bam/bam.h"
7	#include "bam/sam.h"
8
9	bool is_fastq=true; //default is fastq
10	bool sam_input=false; //default is BAM
11	bool all_reads=false;
12	bool mapped_only=false;
13
14	#define USAGE "Usage: bam2fastx [--fasta\|-a\|--fastq\|-q] [--sam\|-s\|-t] \n\
15	[-M\|--mapped-only\|-A\|--all] <in.bam>\|<in.sam> \n"
16
17	char qseq[2048];
18
19	const char *short_options = "aqstMU";
20
21	enum {
22	OPT_FASTA = 127,
23	OPT_FASTQ,
24	OPT_SAM,
25	OPT_MAPPED_ONLY,
26	OPT_ALL
27	};
28
29	struct option long_options[] = {
30	{"fasta", no_argument, 0, OPT_FASTA},
31	{"fastq", no_argument, 0, OPT_FASTQ},
32	{"sam", no_argument, 0, OPT_SAM},
33	{"mapped-only", no_argument, 0, OPT_MAPPED_ONLY},
34	{"all", no_argument, 0, OPT_ALL},
35	{0, 0, 0, 0} // terminator
36	};
37
38	int parse_options(int argc, char** argv)
39	{
40	int option_index = 0;
41	int next_option;
42	do {
43	next_option = getopt_long(argc, argv, short_options, long_options, &option_index);
44	switch (next_option) {
45	case -1:
46	break;
47	case 'a':
48	case OPT_FASTA:
49	is_fastq = false;
50	break;
51	case 'q':
52	case OPT_FASTQ:
53	is_fastq = true;
54	break;
55	case 's':
56	case 't':
57	case OPT_SAM: //sam (text) input
58	sam_input = true;
59	break;
60	case 'M':
61	case OPT_MAPPED_ONLY:
62	mapped_only = true;
63	break;
64	case 'A':
65	case OPT_ALL:
66	all_reads = true;
67	break;
68	default:
69	return 1;
70	}
71	} while(next_option != -1);
72	if (all_reads && mapped_only) {
73	fprintf(stderr, "Error: incompatible options !\n");
74	exit(2);
75	}
76	return 0;
77	}
78
79	#define bam_unmapped(b) (((b)->core.flag & BAM_FUNMAP) != 0)
80
81	void showfastq(const bam1_t b, samfile_t fp) {
82	char *name = bam1_qname(b);
83	char qual = (char)bam1_qual(b);
84	char s = (char)bam1_seq(b);
85	int i;
86	bool ismapped=((b->core.flag & BAM_FUNMAP) == 0);
87	if (ismapped && !all_reads) return;
88	if (mapped_only && !ismapped) return;
89
90	for(i=0;i<(b->core.l_qseq);i++) {
91	int8_t v = bam1_seqi(s,i);
92	qseq[i] = bam_nt16_rev_table[v];
93	}
94	qseq[i] = 0;
95
96	printf("@%s\n%s\n",name, qseq);
97	for(i=0;i<(b->core.l_qseq);i++) {
98	qseq[i]=qual[i]+33;
99	}
100	qseq[i]=0;
101	printf("+\n%s\n",qseq);
102	}
103
104	void showfasta(const bam1_t b, samfile_t fp) {
105	char *name = bam1_qname(b);
106	char s = (char)bam1_seq(b);
107	int i;
108	for(i=0;i<(b->core.l_qseq);i++) {
109	int8_t v = bam1_seqi(s,i);
110	qseq[i] = bam_nt16_rev_table[v];
111	}
112	qseq[i] = 0;
113	printf(">%s\n%s\n",name, qseq);
114	}
115
116
117
118	int main(int argc, char *argv[])
119	{
120	samfile_t *fp;
121	char* fname=NULL;
122
123	if (parse_options(argc, argv) \|\| optind>=argc) {
124	fprintf(stderr, USAGE);
125	return -1;
126	}
127	fname=argv[optind++];
128
129	if (fname==NULL \|\| fname[0]==0) {
130	fprintf(stderr, USAGE);
131	return 1;
132	}
133	if (sam_input)
134	fp = samopen(fname, "r", 0);
135	else
136	fp = samopen(fname, "rb", 0);
137	if (fp == 0) {
138	fprintf(stderr, "Error: bam2fastx failed to open BAM file %s\n", fname);
139	return 1;
140	}
141
142	bam1_t *b = bam_init1();
143	if (is_fastq) {
144	while (samread(fp, b) >= 0) showfastq(b, fp);
145	}
146	else {
147	while (samread(fp, b) >= 0) showfasta(b, fp);
148	}
149
150	bam_destroy1(b);
151
152	samclose(fp);
153	return 0;
154	}