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 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 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 }