ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/bam2fastx.cpp
Revision: 97
Committed: Thu Oct 6 19:45:36 2011 UTC (8 years, 7 months ago) by gpertea
File size: 3421 byte(s)
Log Message:
by default bam2fastx shows unmapped only

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