ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/bamread/bamread.cpp
Revision: 123
Committed: Mon Dec 5 22:21:21 2011 UTC (8 years, 3 months ago) by gpertea
File size: 3385 byte(s)
Log Message:
bamread files

Line User Rev File contents
1 gpertea 123 #include "GBam.h"
2     #include "GArgs.h"
3     #include "GStr.h"
4    
5     #define USAGE "Usage: bamread [--fasta|-a|--fastq|-q|-G|--gff] [--sam|-S] \n\
6     [-M|--mapped-only|-A|--all] [-o <outfile>] <in.bam>|<in.sam> \n"
7    
8    
9     enum OutType {
10     outFASTQ,
11     outFASTA,
12     outGFF
13     };
14    
15     OutType out_type=outFASTQ;
16     bool sam_input=false; //default is BAM
17     bool all_reads=false;
18     bool mapped_only=false;
19     char qseq[2048];
20    
21     #define bam_unmapped(b) (((b)->core.flag & BAM_FUNMAP) != 0)
22    
23     void showfastq(const bam1_t *b, FILE* fout) {
24     char *name = bam1_qname(b);
25     char *qual = (char*)bam1_qual(b);
26     char *s = (char*)bam1_seq(b);
27     int i;
28     bool ismapped=((b->core.flag & BAM_FUNMAP) == 0);
29     if (ismapped && !all_reads) return;
30     if (mapped_only && !ismapped) return;
31    
32     for(i=0;i<(b->core.l_qseq);i++) {
33     int8_t v = bam1_seqi(s,i);
34     qseq[i] = bam_nt16_rev_table[v];
35     }
36     qseq[i] = 0;
37    
38     fprintf(fout, "@%s\n%s\n",name, qseq);
39     for(i=0;i<(b->core.l_qseq);i++) {
40     qseq[i]=qual[i]+33;
41     }
42     qseq[i]=0;
43     fprintf(fout, "+\n%s\n",qseq);
44     }
45    
46     void showfasta(const bam1_t *b, FILE* fout) {
47     char *name = bam1_qname(b);
48     char *s = (char*)bam1_seq(b);
49     int i;
50     for(i=0;i<(b->core.l_qseq);i++) {
51     int8_t v = bam1_seqi(s,i);
52     qseq[i] = bam_nt16_rev_table[v];
53     }
54     qseq[i] = 0;
55     fprintf(fout, ">%s\n%s\n",name, qseq);
56     }
57    
58    
59     void showgff(const bam1_t *b, FILE* fout) {
60     char *name = bam1_qname(b);
61     char *s = (char*)bam1_seq(b);
62     int i;
63     for(i=0;i<(b->core.l_qseq);i++) {
64     int8_t v = bam1_seqi(s,i);
65     qseq[i] = bam_nt16_rev_table[v];
66     }
67     qseq[i] = 0;
68     fprintf(fout, ">%s\n%s\n",name, qseq);
69     }
70    
71    
72    
73     int main(int argc, char *argv[])
74     {
75    
76     GArgs args(argc, argv, "fasta;fastq;gff;all;help;mapped-only;sam;hsSAMGaqo:");
77     args.printError(USAGE, true);
78     if (args.getOpt('h') || args.getOpt("help") || args.startNonOpt()==0) {
79     GMessage(USAGE);
80     return 1;
81     }
82     samfile_t *fp;
83    
84     sam_input=(args.getOpt('S') || args.getOpt('s') || args.getOpt("sam"));
85     all_reads=(args.getOpt('A') || args.getOpt("all"));
86     mapped_only=(args.getOpt('M') || args.getOpt("mapped-only"));
87     if ((all_reads && mapped_only)) {
88     GError("Error: incompatible options !\n");
89     }
90    
91     if (args.getOpt('a') || args.getOpt("fasta"))
92     out_type=outFASTA;
93     else if (args.getOpt('G') || args.getOpt("gff")) {
94     out_type=outGFF;
95     mapped_only=true;
96     }
97     char* fname=args.nextNonOpt();
98     if (fname==NULL || fname[0]==0) {
99     GMessage(USAGE);
100     return 1;
101     }
102     if (sam_input)
103     fp = samopen(fname, "r", 0);
104     else
105     fp = samopen(fname, "rb", 0);
106     if (fp == 0) {
107     fprintf(stderr, "Error: bamread failed to open BAM file %s\n", fname);
108     return 1;
109     }
110     FILE* fout=stdout;
111     char* outfname=args.getOpt('o');
112     if (outfname) {
113     fout=fopen(outfname, "w");
114     if (fout==NULL) {
115     fprintf(stderr, "Error creating output file %s\n", outfname);
116     return 2;
117     }
118     }
119     bam1_t *b = bam_init1();
120     if (out_type==outFASTA) {
121     while (samread(fp, b) >= 0) showfasta(b, fout);
122     }
123     else if (out_type==outGFF) {
124     while (samread(fp, b) >= 0) showgff(b, fout); //showgff here
125     }
126     else {
127     while (samread(fp, b) >= 0) showfastq(b, fout);
128     }
129    
130     bam_destroy1(b);
131    
132     samclose(fp);
133     return 0;
134     }