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 User Rev File contents
1 gpertea 29 #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 gpertea 37 static bool sam_input=false; //default is BAM
11    
12 gpertea 29 static char qseq[2048];
13    
14     const char *short_options = "aq";
15    
16     enum {
17     OPT_FASTA = 127,
18 gpertea 37 OPT_FASTQ,
19     OPT_SAM
20 gpertea 29 };
21    
22     static struct option long_options[] = {
23     {"fasta", no_argument, 0, OPT_FASTA},
24     {"fastq", no_argument, 0, OPT_FASTQ},
25 gpertea 37 {"sam", no_argument, 0, OPT_SAM},
26 gpertea 29 {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 gpertea 37 case 's':
47     case 't':
48     case OPT_SAM: //sam (text) input
49     sam_input = true;
50     break;
51 gpertea 29 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 gpertea 37 #define USAGE "Usage: bam2fastx [--fasta|-a|--fastq|-q] [--sam|-s|-t] <in.bam>|<in.sam> \n"
103 gpertea 29
104     int main(int argc, char *argv[])
105     {
106     samfile_t *fp;
107     char* fname=NULL;
108 gpertea 37
109 gpertea 29 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 gpertea 37 if (sam_input)
120     fp = samopen(fname, "r", 0);
121     else
122     fp = samopen(fname, "rb", 0);
123     if (fp == 0) {
124 gpertea 29 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     }