ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/bam_merge.cpp
(Generate patch)
# Line 2 | Line 2
2   #include <cstdio>
3   #include <cstring>
4  
5 < #include "bam/bam.h"
6 < #include "bam/sam.h"
5 > #include <bam/bam.h>
6 > #include <bam/sam.h>
7 >
8 > #include "common.h"
9   #include "GBase.h"
10   #include "GList.hh"
11  
# Line 11 | Line 13
13  
14   #define ERR_BAM_OPEN "Error: bam_merge failed to open BAM file %s\n"
15  
16 + void print_usage()
17 + {
18 +  fprintf(stderr, USAGE);
19 + }
20 +
21   samfile_t **srcfiles; //array of SAM file handles
22   samfile_t *fw; //output SAM handle
23  
# Line 55 | Line 62
62  
63   int main(int argc, char *argv[])
64   {
65 +  int parse_ret = parse_options(argc, argv, print_usage);
66 +  if (parse_ret)
67 +    return parse_ret;
68 +
69      char* outfname=NULL;
70 <    if (argc<4) {
71 <       fprintf(stderr, USAGE);
70 >    if (argc-optind<3) {
71 >       print_usage();
72         if (argc>1)
73           fprintf(stderr, "Error: only %d arguments given.\n", argc-1);
74         return -1;
75         }
76 <    outfname=argv[1];
77 <    int num_src=argc-2;
76 >    outfname=argv[optind];
77 >    int num_src=argc-optind-1;
78      GMALLOC(srcfiles, (num_src*sizeof(samfile_t*)));
79 <    for (int i=2;i<argc;i++) {
80 <       int fno=i-2;
79 >    for (int i=optind+1;i<argc;i++) {
80 >       int fno=i-optind-1;
81         samfile_t* fp=samopen(argv[i], "rb", 0);
82         if (fp==0) {
83                 fprintf(stderr, ERR_BAM_OPEN, argv[i]);
# Line 84 | Line 95
95      fw=samopen(outfname, "wb", srcfiles[lines[0]->fileno]->header);
96      if (fw==NULL)
97          GError("Error creating output file %s\n", outfname);
98 +
99 +    FILE* findex = NULL;
100 +    if (!index_outfile.empty()) {
101 +      findex = fopen(index_outfile.c_str(), "w");
102 +      if (findex == NULL)
103 +        err_die("Error: cannot create file %s\n", index_outfile.c_str());
104 +    }
105 +
106      int last;
107 +    uint64_t count = 0, last_id = 0;
108      while ((last=lines.Count()-1)>=0) {
109          CBamLine* from=lines[last]; //should have the smallest read_id
110          if (from->fileno<0 || from->b==NULL)
111                    GError("Invalid processTopLine() call with uninitialized value.");
112 <        samwrite(fw, from->b);
112 >        if (findex != NULL)
113 >          {
114 >            ++count;
115 >            if (count >= 1000 && from->read_id != last_id)
116 >              {
117 >                int64_t offset = bam_tell(fw->x.bam);
118 >                int block_offset = offset & 0xFFFF;
119 >                
120 >                // daehwan - this is a bug in bgzf.c in samtools
121 >                // I'll report this bug with a temporary solution, soon.
122 >                if (block_offset <= 0xF000)
123 >                  {
124 >                    fprintf(findex, "%lu\t%ld\n", from->read_id, offset);
125 >                    count = 0;
126 >                  }
127 >              }
128 >            last_id = from->read_id;
129 >          }
130 >        samwrite(fw, from->b);
131 >
132          if (samread(srcfiles[from->fileno], from->b)>0) {
133             from->b_init();
134             //adjust the position in the sorted lines list
# Line 117 | Line 156
156      for (int i=0;i<num_src;i++)
157          samclose(srcfiles[i]);
158      GFREE(srcfiles);
159 +
160 +    if (findex != NULL)
161 +      fclose(findex);
162 +    
163      return 0;
164   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines