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