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, 4 months ago) by gpertea
File size: 3246 byte(s)
Log Message:
working

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] [--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 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 fprintf(fout, "+\n%s\n",qseq);
30 GFREE(qseq);
31 }
32
33 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 }
40
41 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 }
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 /*
81 if (sam_input)
82 fp = samopen(fname, "r", 0);
83 else
84 fp = samopen(fname, "rb", 0);
85 */
86 GBamReader bamreader(fname);
87 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 GBamRecord *aln=NULL;
97 if (out_type==outFASTA) {
98 while ((aln=bamreader.next())!=NULL) {
99 showfasta(*aln, fout);
100 delete aln;
101 }
102 }
103 else if (out_type==outGFF) {
104 while ((aln=bamreader.next())!=NULL) {
105 showgff(*aln, fout);
106 delete aln;
107 }
108 }
109 else {
110 while ((aln=bamreader.next())!=NULL) {
111 showfastq(*aln, fout);
112 delete aln;
113 }
114 }
115 if (fout!=stdout) fclose(fout);
116 return 0;
117 }