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 READ_INT_TYPE N[3]; // note, N = N0 + N1 + N2 , but may not be equal to the total number of reads in data
  32 HIT_INT_TYPE nHits; // # of hits
  33 READ_INT_TYPE nUnique, nMulti, nIsoMulti;
  34 char fn_list[STRLEN];
  35 char groupF[STRLEN], tiF[STRLEN];
  36 char datF[STRLEN], cntF[STRLEN];
  37
  38 GroupInfo gi;
  39 Transcripts transcripts;
  40
  41 SamParser *parser;
  42 ofstream hit_out;
  43
  44 int n_os; // number of ostreams
  45 ostream *cat[3][2]; // cat : category  1-dim 0 N0 1 N1 2 N2; 2-dim  0 mate1 1 mate2
  46 char readOutFs[3][2][STRLEN];
  47
  48 map<int, READ_INT_TYPE> counter;
  49 map<int, READ_INT_TYPE>::iterator iter;
  50
  51 void init(const char* imdName, char alignFType, const char* alignF) {
  52
  53         char* aux = 0;
  54         if (strcmp(fn_list, "")) aux = fn_list;
  55         parser = new SamParser(alignFType, alignF, aux, transcripts, imdName);
  56
  57         memset(cat, 0, sizeof(cat));
  58         memset(readOutFs, 0, sizeof(readOutFs));
  59
  60         int tmp_n_os = -1;
  61
  62         for (int i = 0; i < 3; i++) {
  63                 genReadFileNames(imdName, i, read_type, n_os, readOutFs[i]);
  64
  65                 assert(tmp_n_os < 0 || tmp_n_os == n_os); tmp_n_os = n_os;
  66
  67                 for (int j = 0; j < n_os; j++)
  68                         cat[i][j] = new ofstream(readOutFs[i][j]);
  69         }
  70
  71         counter.clear();
  72 }
  73
  74 //Do not allow duplicate for unalignable reads and supressed reads in SAM input
  75 template<class ReadType, class HitType>
  76 void parseIt(SamParser *parser) {
  77         // record_val & record_read are copies of val & read for record purpose
  78         int val, record_val;
  79         ReadType read, record_read;
  80         HitType hit;
  81         HitContainer<HitType> hits;
  82
  83         nHits = 0;
  84         nUnique = nMulti = nIsoMulti = 0;
  85         memset(N, 0, sizeof(N));
  86
  87         READ_INT_TYPE cnt = 0;
  88
  89         record_val = -2; //indicate no recorded read now
  90         while ((val = parser->parseNext(read, hit)) >= 0) {
  91                 if (val >= 0 && val <= 2) {
  92                         // flush out previous read's info if needed
  93                         if (record_val >= 0) {
  94                                 record_read.write(n_os, cat[record_val]);
  95                                 ++N[record_val];
  96                         }
  97
  98                         general_assert(record_val == 1 || hits.getNHits() == 0, "Read " + record_read.getName() + " is both unalignable and alignable according to the input SAM/BAM file!");
  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)) { cout<< "Parsed "<< cnt<< " entries"<< endl; }
 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 imdName statName 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(groupF, "%s.grp", argv[1]);
 196         gi.load(groupF);
 197         sprintf(tiF, "%s.ti", argv[1]);
 198         transcripts.readFrom(tiF);
 199
 200         sprintf(datF, "%s.dat", argv[2]);
 201         sprintf(cntF, "%s.cnt", argv[3]);
 202
 203         init(argv[2], argv[4][0], argv[5]);
 204
 205         hit_out.open(datF);
 206
 207         string firstLine(99, ' ');
 208         firstLine.append(1, '\n');              //May be dangerous!
 209         hit_out<<firstLine;
 210
 211         switch(read_type) {
 212         case 0 : parseIt<SingleRead, SingleHit>(parser); break;
 213         case 1 : parseIt<SingleReadQ, SingleHit>(parser); break;
 214         case 2 : parseIt<PairedEndRead, PairedEndHit>(parser); break;
 215         case 3 : parseIt<PairedEndReadQ, PairedEndHit>(parser); break;
 216         }
 217
 218         hit_out.seekp(0, ios_base::beg);
 219         hit_out<<N[1]<<" "<<nHits<<" "<<read_type;
 220
 221         hit_out.close();
 222
 223         //cntF for statistics of alignments file
 224         ofstream fout(cntF);
 225         fout<<N[0]<<" "<<N[1]<<" "<<N[2]<<" "<<(N[0] + N[1] + N[2])<<endl;
 226         fout<<nUnique<<" "<<nMulti<<" "<<nIsoMulti<<endl;
 227         fout<<nHits<<" "<<read_type<<endl;
 228         fout<<"0\t"<<N[0]<<endl;
 229         for (iter = counter.begin(); iter != counter.end(); iter++) {
 230                 fout<<iter->first<<'\t'<<iter->second<<endl;
 231         }
 232         fout<<"Inf\t"<<N[2]<<endl;
 233         fout.close();
 234
 235         release();
 236
 237         if (verbose) { printf("Done!\n"); }
 238
 239         return 0;
 240 }