ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/bamread/bamread.cpp
Revision: 131
Committed: Thu Dec 8 02:53:00 2011 UTC (8 years, 3 months ago) by gpertea
File size: 3246 byte(s)
Log Message:
working

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    
20     #define bam_unmapped(b) (((b)->core.flag & BAM_FUNMAP) != 0)
21    
22 gpertea 128 void showfastq(GBamRecord& rec, FILE* fout) {
23     if (rec.isMapped() && !all_reads) return;
24     if (mapped_only && rec.isUnmapped()) return;
25     char* qseq=rec.sequence();
26     fprintf(fout, "@%s\n%s\n", rec.name(), qseq);
27     GFREE(qseq);
28     qseq=rec.qualities();
29 gpertea 123 fprintf(fout, "+\n%s\n",qseq);
30 gpertea 128 GFREE(qseq);
31 gpertea 123 }
32    
33 gpertea 128 void showfasta(GBamRecord& rec, FILE* fout) {
34     if (rec.isMapped() && !all_reads) return;
35     if (mapped_only && rec.isUnmapped()) return;
36     char* qseq=rec.sequence();
37     fprintf(fout, ">%s\n%s\n", rec.name(), qseq);
38     GFREE(qseq);
39 gpertea 123 }
40    
41 gpertea 128 void showgff(GBamRecord& rec, FILE* fout) {
42     if (rec.isUnmapped()) return;
43     char tstrand=rec.spliceStrand();
44     fprintf(fout, "%s\tbam\tmRNA\t%d\t%d\t.\t%c\t.\tID=%s\n", rec.refName(),
45     rec.start, rec.end, tstrand, rec.name());
46     for (int i=0;i<rec.exons.Count();i++) {
47     fprintf(fout, "%s\tbam\texon\t%d\t%d\t.\t%c\t.\tParent=%s\n", rec.refName(),
48     rec.exons[i].start, rec.exons[i].end, tstrand, rec.name());
49     }
50 gpertea 123 }
51    
52     int main(int argc, char *argv[])
53     {
54    
55     GArgs args(argc, argv, "fasta;fastq;gff;all;help;mapped-only;sam;hsSAMGaqo:");
56     args.printError(USAGE, true);
57     if (args.getOpt('h') || args.getOpt("help") || args.startNonOpt()==0) {
58     GMessage(USAGE);
59     return 1;
60     }
61    
62     sam_input=(args.getOpt('S') || args.getOpt('s') || args.getOpt("sam"));
63     all_reads=(args.getOpt('A') || args.getOpt("all"));
64     mapped_only=(args.getOpt('M') || args.getOpt("mapped-only"));
65     if ((all_reads && mapped_only)) {
66     GError("Error: incompatible options !\n");
67     }
68    
69     if (args.getOpt('a') || args.getOpt("fasta"))
70     out_type=outFASTA;
71     else if (args.getOpt('G') || args.getOpt("gff")) {
72     out_type=outGFF;
73     mapped_only=true;
74     }
75     char* fname=args.nextNonOpt();
76     if (fname==NULL || fname[0]==0) {
77     GMessage(USAGE);
78     return 1;
79     }
80 gpertea 128 /*
81 gpertea 123 if (sam_input)
82     fp = samopen(fname, "r", 0);
83     else
84     fp = samopen(fname, "rb", 0);
85 gpertea 128 */
86     GBamReader bamreader(fname);
87 gpertea 123 FILE* fout=stdout;
88     char* outfname=args.getOpt('o');
89     if (outfname) {
90     fout=fopen(outfname, "w");
91     if (fout==NULL) {
92     fprintf(stderr, "Error creating output file %s\n", outfname);
93     return 2;
94     }
95     }
96 gpertea 128 GBamRecord *aln=NULL;
97 gpertea 123 if (out_type==outFASTA) {
98 gpertea 128 while ((aln=bamreader.next())!=NULL) {
99     showfasta(*aln, fout);
100     delete aln;
101     }
102 gpertea 123 }
103     else if (out_type==outGFF) {
104 gpertea 128 while ((aln=bamreader.next())!=NULL) {
105     showgff(*aln, fout);
106     delete aln;
107     }
108 gpertea 123 }
109     else {
110 gpertea 128 while ((aln=bamreader.next())!=NULL) {
111     showfastq(*aln, fout);
112     delete aln;
113     }
114 gpertea 123 }
115 gpertea 131 if (fout!=stdout) fclose(fout);
116 gpertea 123 return 0;
117     }