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