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 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     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     }