ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/bamread/bamread.cpp
Revision: 134
Committed: Fri Dec 9 17:52:20 2011 UTC (8 years, 3 months ago) by gpertea
File size: 2993 byte(s)
Log Message:
Line User Rev File contents
1 gpertea 123 #include "GBam.h"
2     #include "GArgs.h"
3     #include "GStr.h"
4    
5 gpertea 134 #define USAGE "Usage: bamread [--fasta|-a|--fastq|-q|-G|--gff]\n\
6 gpertea 123 [-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 all_reads=false;
17     bool mapped_only=false;
18    
19     #define bam_unmapped(b) (((b)->core.flag & BAM_FUNMAP) != 0)
20    
21 gpertea 128 void showfastq(GBamRecord& rec, FILE* fout) {
22     if (rec.isMapped() && !all_reads) return;
23     if (mapped_only && rec.isUnmapped()) return;
24     char* qseq=rec.sequence();
25     fprintf(fout, "@%s\n%s\n", rec.name(), qseq);
26     GFREE(qseq);
27     qseq=rec.qualities();
28 gpertea 123 fprintf(fout, "+\n%s\n",qseq);
29 gpertea 128 GFREE(qseq);
30 gpertea 123 }
31    
32 gpertea 128 void showfasta(GBamRecord& rec, FILE* fout) {
33     if (rec.isMapped() && !all_reads) return;
34     if (mapped_only && rec.isUnmapped()) return;
35     char* qseq=rec.sequence();
36     fprintf(fout, ">%s\n%s\n", rec.name(), qseq);
37     GFREE(qseq);
38 gpertea 123 }
39    
40 gpertea 128 void showgff(GBamRecord& rec, FILE* fout) {
41     if (rec.isUnmapped()) return;
42     char tstrand=rec.spliceStrand();
43     fprintf(fout, "%s\tbam\tmRNA\t%d\t%d\t.\t%c\t.\tID=%s\n", rec.refName(),
44     rec.start, rec.end, tstrand, rec.name());
45     for (int i=0;i<rec.exons.Count();i++) {
46     fprintf(fout, "%s\tbam\texon\t%d\t%d\t.\t%c\t.\tParent=%s\n", rec.refName(),
47     rec.exons[i].start, rec.exons[i].end, tstrand, rec.name());
48     }
49 gpertea 123 }
50    
51     int main(int argc, char *argv[])
52     {
53    
54 gpertea 134 GArgs args(argc, argv, "fasta;fastq;gff;all;help;mapped-only;hAMGaqo:");
55 gpertea 123 args.printError(USAGE, true);
56     if (args.getOpt('h') || args.getOpt("help") || args.startNonOpt()==0) {
57     GMessage(USAGE);
58     return 1;
59     }
60    
61     all_reads=(args.getOpt('A') || args.getOpt("all"));
62     mapped_only=(args.getOpt('M') || args.getOpt("mapped-only"));
63     if ((all_reads && mapped_only)) {
64     GError("Error: incompatible options !\n");
65     }
66    
67     if (args.getOpt('a') || args.getOpt("fasta"))
68     out_type=outFASTA;
69     else if (args.getOpt('G') || args.getOpt("gff")) {
70     out_type=outGFF;
71     mapped_only=true;
72     }
73     char* fname=args.nextNonOpt();
74     if (fname==NULL || fname[0]==0) {
75     GMessage(USAGE);
76     return 1;
77     }
78 gpertea 128 GBamReader bamreader(fname);
79 gpertea 123 FILE* fout=stdout;
80     char* outfname=args.getOpt('o');
81     if (outfname) {
82     fout=fopen(outfname, "w");
83     if (fout==NULL) {
84     fprintf(stderr, "Error creating output file %s\n", outfname);
85     return 2;
86     }
87     }
88 gpertea 128 GBamRecord *aln=NULL;
89 gpertea 123 if (out_type==outFASTA) {
90 gpertea 128 while ((aln=bamreader.next())!=NULL) {
91     showfasta(*aln, fout);
92     delete aln;
93     }
94 gpertea 123 }
95     else if (out_type==outGFF) {
96 gpertea 128 while ((aln=bamreader.next())!=NULL) {
97     showgff(*aln, fout);
98     delete aln;
99     }
100 gpertea 123 }
101     else {
102 gpertea 128 while ((aln=bamreader.next())!=NULL) {
103     showfastq(*aln, fout);
104     delete aln;
105     }
106 gpertea 123 }
107 gpertea 131 if (fout!=stdout) fclose(fout);
108 gpertea 123 return 0;
109     }