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 (7 years, 8 months ago) by gpertea
File size: 2993 byte(s)
Log Message:
Line File contents
1 #include "GBam.h"
2 #include "GArgs.h"
3 #include "GStr.h"
4
5 #define USAGE "Usage: bamread [--fasta|-a|--fastq|-q|-G|--gff]\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 all_reads=false;
17 bool mapped_only=false;
18
19 #define bam_unmapped(b) (((b)->core.flag & BAM_FUNMAP) != 0)
20
21 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 fprintf(fout, "+\n%s\n",qseq);
29 GFREE(qseq);
30 }
31
32 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 }
39
40 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 }
50
51 int main(int argc, char *argv[])
52 {
53
54 GArgs args(argc, argv, "fasta;fastq;gff;all;help;mapped-only;hAMGaqo:");
55 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 GBamReader bamreader(fname);
79 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 GBamRecord *aln=NULL;
89 if (out_type==outFASTA) {
90 while ((aln=bamreader.next())!=NULL) {
91 showfasta(*aln, fout);
92 delete aln;
93 }
94 }
95 else if (out_type==outGFF) {
96 while ((aln=bamreader.next())!=NULL) {
97 showgff(*aln, fout);
98 delete aln;
99 }
100 }
101 else {
102 while ((aln=bamreader.next())!=NULL) {
103 showfastq(*aln, fout);
104 delete aln;
105 }
106 }
107 if (fout!=stdout) fclose(fout);
108 return 0;
109 }