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 "Transcript.h"
  18 #include "Transcripts.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 groupF[STRLEN], tiF[STRLEN];
  38 char imdName[STRLEN];
  39 char datF[STRLEN], cntF[STRLEN];
  40
  41 GroupInfo gi;
  42 Transcripts transcripts;
  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, transcripts, 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                                 hits.updateRI();
 100                                 nHits += hits.getNHits();
 101                                 nMulti += hits.calcNumGeneMultiReads(gi);
 102                                 nIsoMulti += hits.calcNumIsoformMultiReads();
 103                                 hits.write(hit_out);
 104
 105                                 iter = counter.find(hits.getNHits());
 106                                 if (iter != counter.end()) {
 107                                         iter->second++;
 108                                 }
 109                                 else {
 110                                         counter[hits.getNHits()] = 1;
 111                                 }
 112                         }
 113
 114                         hits.clear();
 115                         record_val = val;
 116                         record_read = read; // no pointer, thus safe
 117                 }
 118
 119                 if (val == 1 || val == 5) {
 120                         hits.push_back(hit);
 121                 }
 122
 123                 ++cnt;
 124                 if (verbose && (cnt % 1000000 == 0)) { printf("Parsed %lld entries\n", cnt); }
 125         }
 126
 127         if (record_val >= 0) {
 128                 record_read.write(n_os, cat[record_val]);
 129                 ++N[record_val];
 130                 hits.updateRI();
 131                 nHits += hits.getNHits();
 132                 nMulti += hits.calcNumGeneMultiReads(gi);
 133                 nIsoMulti += hits.calcNumIsoformMultiReads();
 134                 hits.write(hit_out);
 135
 136                 iter = counter.find(hits.getNHits());
 137                 if (iter != counter.end()) {
 138                         iter->second++;
 139                 }
 140                 else {
 141                         counter[hits.getNHits()] = 1;
 142                 }
 143         }
 144
 145         nUnique = N[1] - nMulti;
 146 }
 147
 148 void release() {
 149         for (int i = 0; i < 3; i++) {
 150                 for (int j = 0; j < n_os; j++) {
 151                         ((ofstream*)cat[i][j])->close();
 152                         delete cat[i][j];
 153                 }
 154                 if (N[i] > 0) continue;
 155                 for (int j = 0; j < n_os; j++) {
 156                         remove(readOutFs[i][j]); //delete if the file is empty
 157                 }
 158         }
 159         delete parser;
 160 }
 161
 162 int main(int argc, char* argv[]) {
 163         bool quiet = false;
 164
 165         if (argc < 6) {
 166                 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");
 167                 exit(-1);
 168         }
 169
 170         strcpy(fn_list, "");
 171         read_type = 0;
 172         if (argc > 6) {
 173                 for (int i = 6; i < argc; i++) {
 174                         if (!strcmp(argv[i], "-t")) {
 175                                 read_type = atoi(argv[i + 1]);
 176                         }
 177                         if (!strcmp(argv[i], "-l")) {
 178                                 strcpy(fn_list, argv[i + 1]);
 179                         }
 180                         if (!strcmp(argv[i], "-tag")) {
 181                                 SamParser::setReadTypeTag(argv[i + 1]);
 182                         }
 183                         if (!strcmp(argv[i], "-q")) { quiet = true; }
 184                 }
 185         }
 186
 187         verbose = !quiet;
 188
 189         sprintf(groupF, "%s.grp", argv[1]);
 190         gi.load(groupF);
 191         sprintf(tiF, "%s.ti", argv[1]);
 192         transcripts.readFrom(tiF);
 193
 194         sprintf(imdName, "%s.temp/%s", argv[2], argv[3]);
 195         sprintf(datF, "%s.dat", imdName);
 196         sprintf(cntF, "%s.stat/%s.cnt", argv[2], argv[3]);
 197
 198         init(imdName, argv[4][0], argv[5]);
 199
 200         hit_out.open(datF);
 201
 202         string firstLine(59, ' ');
 203         firstLine.append(1, '\n');              //May be dangerous!
 204         hit_out<<firstLine;
 205
 206         switch(read_type) {
 207         case 0 : parseIt<SingleRead, SingleHit>(parser); break;
 208         case 1 : parseIt<SingleReadQ, SingleHit>(parser); break;
 209         case 2 : parseIt<PairedEndRead, PairedEndHit>(parser); break;
 210         case 3 : parseIt<PairedEndReadQ, PairedEndHit>(parser); break;
 211         }
 212
 213         hit_out.seekp(0, ios_base::beg);
 214         hit_out<<N[1]<<" "<<nHits<<" "<<read_type;
 215
 216         hit_out.close();
 217
 218         //cntF for statistics of alignments file
 219         ofstream fout(cntF);
 220         fout<<N[0]<<" "<<N[1]<<" "<<N[2]<<" "<<(N[0] + N[1] + N[2])<<endl;
 221         fout<<nUnique<<" "<<nMulti<<" "<<nIsoMulti<<endl;
 222         fout<<nHits<<" "<<read_type<<endl;
 223         for (iter = counter.begin(); iter != counter.end(); iter++) {
 224                 fout<<iter->first<<'\t'<<iter->second<<endl;
 225         }
 226         fout.close();
 227
 228         release();
 229
 230         if (verbose) { printf("Done!\n"); }
 231
 232         return 0;
 233 }