parseIt.cpp

   1 /*
   2  * Assume any read should have a name other than ""
   3  */
   4 #include<cstdio>
   5 #include<cstring>
   6 #include<cstdlib>
   7 #include<cassert>
   8 #include<iostream>
   9 #include<fstream>
  10 #include<string>
  11
  12 #include "utils.h"
  13
  14 #include "GroupInfo.h"
  15
  16 #include "SingleRead.h"
  17 #include "SingleReadQ.h"
  18 #include "PairedEndRead.h"
  19 #include "PairedEndReadQ.h"
  20 #include "SingleHit.h"
  21 #include "PairedEndHit.h"
  22
  23 #include "HitContainer.h"
  24 #include "SamParser.h"
  25
  26 using namespace std;
  27
  28 int read_type; // 0 SingleRead, 1 SingleReadQ, 2 PairedEndRead, 3 PairedEndReadQ
  29 int N[3]; // note, N = N0 + N1 + N2 , but may not be equal to the total number of reads in data
  30 int nHits; // # of hits
  31 int nUnique, nMulti, nIsoMulti;
  32 char fn_list[STRLEN];
  33 char groupF[STRLEN];
  34 char datF[STRLEN], cntF[STRLEN];
  35
  36 GroupInfo gi;
  37
  38 SamParser *parser;
  39 ofstream hit_out;
  40
  41 int n_os; // number of ostreams
  42 ostream *cat[3][2]; // cat : category  1-dim 0 N0 1 N1 2 N2; 2-dim  0 mate1 1 mate2
  43 char readOutFs[3][2][STRLEN];
  44
  45 void init(const char* imdName, char alignFType, const char* alignF) {
  46
  47         sprintf(datF, "%s.dat", imdName);
  48         sprintf(cntF, "%s.cnt", imdName);
  49
  50         char* aux = 0;
  51         if (strcmp(fn_list, "")) aux = fn_list;
  52         parser = new SamParser(alignFType, alignF, aux);
  53
  54         memset(cat, 0, sizeof(cat));
  55         memset(readOutFs, 0, sizeof(readOutFs));
  56
  57         int tmp_n_os = -1;
  58
  59         for (int i = 0; i < 3; i++) {
  60                 genReadFileNames(imdName, i, read_type, n_os, readOutFs[i]);
  61
  62                 assert(tmp_n_os < 0 || tmp_n_os == n_os); tmp_n_os = n_os;
  63
  64                 for (int j = 0; j < n_os; j++)
  65                         cat[i][j] = new ofstream(readOutFs[i][j]);
  66         }
  67 }
  68
  69 //Do not allow duplicate for unalignable reads and supressed reads in SAM input
  70 template<class ReadType, class HitType>
  71 void parseIt(SamParser *parser) {
  72         // record_val & record_read are copies of val & read for record purpose
  73         int val, record_val;
  74         ReadType read, record_read;
  75         HitType hit;
  76         HitContainer<HitType> hits;
  77
  78         nHits = 0;
  79         nUnique = nMulti = nIsoMulti = 0;
  80         memset(N, 0, sizeof(N));
  81
  82         long long cnt = 0;
  83
  84         record_val = -2; //indicate no recorded read now
  85         while ((val = parser->parseNext(read, hit)) >= 0) {
  86                 if (val >= 0 && val <= 2) {
  87                         // flush out previous read's info if needed
  88                         if (record_val >= 0) {
  89                                 record_read.write(n_os, cat[record_val]);
  90                                 ++N[record_val];
  91                                 hits.updateRI();
  92                                 nHits += hits.getNHits();
  93                                 nMulti += hits.calcNumGeneMultiReads(gi);
  94                                 nIsoMulti += hits.calcNumIsoformMultiReads();
  95                                 hits.write(hit_out);
  96                         }
  97
  98                         hits.clear();
  99                         record_val = val;
 100                         record_read = read; // no pointer, thus safe
 101                 }
 102
 103                 if (val == 1 || val == 5) {
 104                         hits.push_back(hit);
 105                 }
 106
 107                 ++cnt;
 108                 if (verbose && (cnt % 1000000 == 0)) { printf("Parsed %lld entries\n", cnt); }
 109         }
 110
 111         if (record_val >= 0) {
 112                 record_read.write(n_os, cat[record_val]);
 113                 ++N[record_val];
 114                 hits.updateRI();
 115                 nHits += hits.getNHits();
 116                 nMulti += hits.calcNumGeneMultiReads(gi);
 117                 nIsoMulti += hits.calcNumIsoformMultiReads();
 118                 hits.write(hit_out);
 119         }
 120
 121         nUnique = N[1] - nMulti;
 122 }
 123
 124 void release() {
 125         for (int i = 0; i < 3; i++) {
 126                 for (int j = 0; j < n_os; j++) {
 127                         ((ofstream*)cat[i][j])->close();
 128                         delete cat[i][j];
 129                 }
 130                 if (N[i] > 0) continue;
 131                 for (int j = 0; j < n_os; j++) {
 132                         remove(readOutFs[i][j]); //delete if the file is empty
 133                 }
 134         }
 135         delete parser;
 136 }
 137
 138 int main(int argc, char* argv[]) {
 139         bool quiet = false;
 140
 141         if (argc < 5) {
 142                 printf("Usage : rsem-parse-alignments refName imdName alignFType('s' for sam, 'b' for bam) alignF [-t Type] [-l fn_list] [-tag tagName] [-q]\n");
 143                 exit(-1);
 144         }
 145
 146         strcpy(fn_list, "");
 147         read_type = 0;
 148         if (argc > 5) {
 149                 for (int i = 5; i < argc; i++) {
 150                         if (!strcmp(argv[i], "-t")) {
 151                                 read_type = atoi(argv[i + 1]);
 152                         }
 153                         if (!strcmp(argv[i], "-l")) {
 154                                 strcpy(fn_list, argv[i + 1]);
 155                         }
 156                         if (!strcmp(argv[i], "-tag")) {
 157                                 SamParser::setReadTypeTag(argv[i + 1]);
 158                         }
 159                         if (!strcmp(argv[i], "-q")) { quiet = true; }
 160                 }
 161         }
 162
 163         verbose = !quiet;
 164
 165         init(argv[2], argv[3][0], argv[4]);
 166
 167         sprintf(groupF, "%s.grp", argv[1]);
 168         gi.load(groupF);
 169
 170         hit_out.open(datF);
 171
 172         string firstLine(59, ' ');
 173         firstLine.append(1, '\n');              //May be dangerous!
 174         hit_out<<firstLine;
 175
 176         switch(read_type) {
 177         case 0 : parseIt<SingleRead, SingleHit>(parser); break;
 178         case 1 : parseIt<SingleReadQ, SingleHit>(parser); break;
 179         case 2 : parseIt<PairedEndRead, PairedEndHit>(parser); break;
 180         case 3 : parseIt<PairedEndReadQ, PairedEndHit>(parser); break;
 181         }
 182
 183         hit_out.seekp(0, ios_base::beg);
 184         hit_out<<N[1]<<" "<<nHits<<" "<<read_type;
 185
 186         hit_out.close();
 187
 188         //cntF for statistics of alignments file
 189         ofstream fout(cntF);
 190         fout<<N[0]<<" "<<N[1]<<" "<<N[2]<<" "<<(N[0] + N[1] + N[2])<<endl;
 191         fout<<nUnique<<" "<<nMulti<<" "<<nIsoMulti<<endl;
 192         fout<<nHits<<" "<<read_type<<endl;
 193         fout.close();
 194
 195         release();
 196
 197         if (verbose) { printf("Done!\n"); }
 198
 199         return 0;
 200 }