2 * Assume any read should have a name other than ""
14 #include "GroupInfo.h"
16 #include "SingleRead.h"
17 #include "SingleReadQ.h"
18 #include "PairedEndRead.h"
19 #include "PairedEndReadQ.h"
20 #include "SingleHit.h"
21 #include "PairedEndHit.h"
23 #include "HitContainer.h"
24 #include "SamParser.h"
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;
34 char datF[STRLEN], cntF[STRLEN];
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];
45 void init(const char* imdName, char alignFType, const char* alignF) {
47 sprintf(datF, "%s.dat", imdName);
48 sprintf(cntF, "%s.cnt", imdName);
51 if (strcmp(fn_list, "")) aux = fn_list;
52 parser = new SamParser(alignFType, alignF, aux);
54 memset(cat, 0, sizeof(cat));
55 memset(readOutFs, 0, sizeof(readOutFs));
59 for (int i = 0; i < 3; i++) {
60 genReadFileNames(imdName, i, read_type, n_os, readOutFs[i]);
62 assert(tmp_n_os < 0 || tmp_n_os == n_os); tmp_n_os = n_os;
64 for (int j = 0; j < n_os; j++)
65 cat[i][j] = new ofstream(readOutFs[i][j]);
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
74 ReadType read, record_read;
76 HitContainer<HitType> hits;
79 nUnique = nMulti = nIsoMulti = 0;
80 memset(N, 0, sizeof(N));
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]);
92 nHits += hits.getNHits();
93 nMulti += hits.calcNumGeneMultiReads(gi);
94 nIsoMulti += hits.calcNumIsoformMultiReads();
100 record_read = read; // no pointer, thus safe
103 if (val == 1 || val == 5) {
108 if (verbose && (cnt % 1000000 == 0)) { printf("Parsed %lld entries\n", cnt); }
111 if (record_val >= 0) {
112 record_read.write(n_os, cat[record_val]);
115 nHits += hits.getNHits();
116 nMulti += hits.calcNumGeneMultiReads(gi);
117 nIsoMulti += hits.calcNumIsoformMultiReads();
121 nUnique = N[1] - nMulti;
125 for (int i = 0; i < 3; i++) {
126 for (int j = 0; j < n_os; j++) {
127 ((ofstream*)cat[i][j])->close();
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
138 int main(int argc, char* argv[]) {
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");
149 for (int i = 5; i < argc; i++) {
150 if (!strcmp(argv[i], "-t")) {
151 read_type = atoi(argv[i + 1]);
153 if (!strcmp(argv[i], "-l")) {
154 strcpy(fn_list, argv[i + 1]);
156 if (!strcmp(argv[i], "-tag")) {
157 SamParser::setReadTypeTag(argv[i + 1]);
159 if (!strcmp(argv[i], "-q")) { quiet = true; }
165 init(argv[2], argv[3][0], argv[4]);
167 sprintf(groupF, "%s.grp", argv[1]);
172 string firstLine(59, ' ');
173 firstLine.append(1, '\n'); //May be dangerous!
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;
183 hit_out.seekp(0, ios_base::beg);
184 hit_out<<N[1]<<" "<<nHits<<" "<<read_type;
188 //cntF for statistics of alignments file
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;
197 if (verbose) { printf("Done!\n"); }