getUnique.cpp

   1 #include<cstdio>
   2 #include<cstring>
   3 #include<cstdlib>
   4 #include<cassert>
   5 #include<string>
   6 #include<vector>
   7
   8 #include <stdint.h>
   9 #include "sam/bam.h"
  10 #include "sam/sam.h"
  11
  12 #include "utils.h"
  13
  14 using namespace std;
  15
  16 string cqname;
  17 samfile_t *in, *out;
  18 bam1_t *b;
  19 vector<bam1_t*> arr;
  20 bool unaligned;
  21
  22 void output() {
  23         if (unaligned || arr.size() == 0) return;
  24         bool isPaired = (arr[0]->core.flag & 0x0001);
  25         if ((isPaired && arr.size() != 2) || (!isPaired && arr.size() != 1)) return;
  26         for (size_t i = 0; i < arr.size(); i++) samwrite(out, arr[i]);
  27 }
  28
  29 int main(int argc, char* argv[]) {
  30         if (argc != 3) {
  31                 printf("Usage: rsem-get-unique unsorted_transcript_bam_input bam_output\n");
  32                 exit(-1);
  33         }
  34
  35         in = samopen(argv[1], "rb", NULL);
  36         assert(in != 0);
  37         out = samopen(argv[2], "wb", in->header);
  38         assert(out != 0);
  39
  40         HIT_INT_TYPE cnt = 0;
  41
  42         cqname = "";
  43         arr.clear();
  44         b = bam_init1();
  45         unaligned = false;
  46
  47         while (samread(in, b) >= 0) {
  48                 if (cqname != bam1_qname(b)) {
  49                         output();
  50                         cqname = bam1_qname(b);
  51                         for (size_t i = 0; i < arr.size(); i++) bam_destroy1(arr[i]);
  52                         arr.clear();
  53                         unaligned = false;
  54                 }
  55
  56                 unaligned = unaligned || (b->core.flag & 0x0004);
  57                 arr.push_back(bam_dup1(b));
  58
  59                 ++cnt;
  60                 if (cnt % 1000000 == 0) { printf("."); fflush(stdout); }
  61         }
  62
  63         if (cnt >= 1000000) printf("\n");
  64
  65         output();
  66
  67         bam_destroy1(b);
  68         samclose(in);
  69         samclose(out);
  70
  71         printf("done!\n");
  72
  73         return 0;
  74 }