ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/bam2fastx.cpp
Revision: 37
Committed: Mon Aug 22 16:17:34 2011 UTC (8 years, 9 months ago) by gpertea
File size: 3056 byte(s)
Log Message:
wip

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 static bool is_fastq=true; //default is fastq
10 static bool sam_input=false; //default is BAM
11
12 static char qseq[2048];
13
14 const char *short_options = "aq";
15
16 enum {
17 OPT_FASTA = 127,
18 OPT_FASTQ,
19 OPT_SAM
20 };
21
22 static struct option long_options[] = {
23 {"fasta", no_argument, 0, OPT_FASTA},
24 {"fastq", no_argument, 0, OPT_FASTQ},
25 {"sam", no_argument, 0, OPT_SAM},
26 {0, 0, 0, 0} // terminator
27 };
28
29 int parse_options(int argc, char** argv)
30 {
31 int option_index = 0;
32 int next_option;
33 do {
34 next_option = getopt_long(argc, argv, short_options, long_options, &option_index);
35 switch (next_option) {
36 case -1:
37 break;
38 case 'a':
39 case OPT_FASTA:
40 is_fastq = false;
41 break;
42 case 'q':
43 case OPT_FASTQ:
44 is_fastq = true;
45 break;
46 case 's':
47 case 't':
48 case OPT_SAM: //sam (text) input
49 sam_input = true;
50 break;
51 default:
52 return 1;
53 }
54 } while(next_option != -1);
55
56 return 0;
57 }
58
59
60 void showfastq(const bam1_t *b, samfile_t* fp) {
61 //uint32_t *cigar = bam1_cigar(b);
62 //const bam1_core_t *c = &b->core;
63 char *name = bam1_qname(b);
64 char *qual = (char*)bam1_qual(b);
65 char *s = (char*)bam1_seq(b);
66 int i;
67 for(i=0;i<(b->core.l_qseq);i++) {
68 int8_t v = bam1_seqi(s,i);
69 qseq[i] = bam_nt16_rev_table[v];
70 }
71 qseq[i] = 0;
72
73 printf("@%s\n%s\n",name, qseq);
74 //print SAM text format of this BAM record:
75 //char * samtext=bam_format1(fp_in->header, b);
76 //printf("%s\n",samtext);
77
78 //printf("s cigar: %s\n",cigar);
79
80 //printf("qual : ");
81 for(i=0;i<(b->core.l_qseq);i++) {
82 //printf(" %d",qual[n]);
83 qseq[i]=qual[i]+33;
84 }
85 qseq[i]=0;
86 printf("+\n%s\n",qseq);
87 }
88
89 void showfasta(const bam1_t *b, samfile_t* fp) {
90 char *name = bam1_qname(b);
91 char *s = (char*)bam1_seq(b);
92 int i;
93 for(i=0;i<(b->core.l_qseq);i++) {
94 int8_t v = bam1_seqi(s,i);
95 qseq[i] = bam_nt16_rev_table[v];
96 }
97 qseq[i] = 0;
98 printf(">%s\n%s\n",name, qseq);
99 }
100
101
102 #define USAGE "Usage: bam2fastx [--fasta|-a|--fastq|-q] [--sam|-s|-t] <in.bam>|<in.sam> \n"
103
104 int main(int argc, char *argv[])
105 {
106 samfile_t *fp;
107 char* fname=NULL;
108
109 if (parse_options(argc, argv) || optind>=argc) {
110 fprintf(stderr, USAGE);
111 return -1;
112 }
113 fname=argv[optind++];
114
115 if (fname==NULL || fname[0]==0) {
116 fprintf(stderr, USAGE);
117 return 1;
118 }
119 if (sam_input)
120 fp = samopen(fname, "r", 0);
121 else
122 fp = samopen(fname, "rb", 0);
123 if (fp == 0) {
124 fprintf(stderr, "Error: bam2fastx failed to open BAM file %s\n", fname);
125 return 1;
126 }
127
128 bam1_t *b = bam_init1();
129 if (is_fastq) {
130 while (samread(fp, b) >= 0) showfastq(b, fp);
131 }
132 else {
133 while (samread(fp, b) >= 0) showfasta(b, fp);
134 }
135
136 bam_destroy1(b);
137
138 samclose(fp);
139 return 0;
140 }