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 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
18
19 char qseq[2048];
20
21 const char *short_options = "aqstMU";
22
23 enum {
24 OPT_FASTA = 127,
25 OPT_FASTQ,
26 OPT_SAM,
27 OPT_MAPPED_ONLY,
28 OPT_ALL
29 };
30
31 struct option long_options[] = {
32 {"fasta", no_argument, 0, OPT_FASTA},
33 {"fastq", no_argument, 0, OPT_FASTQ},
34 {"sam", no_argument, 0, OPT_SAM},
35 {"mapped-only", no_argument, 0, OPT_MAPPED_ONLY},
36 {"all", no_argument, 0, OPT_ALL},
37 {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 case 's':
58 case 't':
59 case OPT_SAM: //sam (text) input
60 sam_input = true;
61 break;
62 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 default:
71 return 1;
72 }
73 } while(next_option != -1);
74 if (all_reads && mapped_only) {
75 fprintf(stderr, "Error: incompatible options !\n");
76 exit(2);
77 }
78 return 0;
79 }
80
81 #define bam_unmapped(b) (((b)->core.flag & BAM_FUNMAP) != 0)
82
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 bool ismapped=((b->core.flag & BAM_FUNMAP) == 0);
89 if (ismapped && !all_reads) return;
90 if (mapped_only && !ismapped) return;
91
92 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 }
96 qseq[i] = 0;
97
98 printf("@%s\n%s\n",name, qseq);
99 for(i=0;i<(b->core.l_qseq);i++) {
100 qseq[i]=qual[i]+33;
101 }
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
125 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 if (sam_input)
136 fp = samopen(fname, "r", 0);
137 else
138 fp = samopen(fname, "rb", 0);
139 if (fp == 0) {
140 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 }