ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/bam2fastx.cpp
Revision: 29
Committed: Tue Aug 2 21:24:54 2011 UTC (8 years, 9 months ago) by gpertea
File size: 2751 byte(s)
Log Message:
adding tophat source work

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