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