synthesisRef.cpp

   1 #include<cstdio>
   2 #include<cstring>
   3 #include<cstdlib>
   4 #include<cassert>
   5 #include<fstream>
   6 #include<sstream>
   7 #include<map>
   8
   9 #include "utils.h"
  10 #include "Transcript.h"
  11 #include "Transcripts.h"
  12
  13 using namespace std;
  14
  15 int M;
  16
  17 map<string, string> name2seq;
  18 map<string, string>::iterator iter;
  19
  20 Transcripts transcripts;
  21 char groupF[STRLEN], tiF[STRLEN], refFastaF[STRLEN];
  22
  23 bool hasMappingFile;
  24 char mappingFile[STRLEN];
  25
  26 map<string, string> mi_table; // mapping info table
  27 map<string, string>::iterator mi_iter; //mapping info table's iterator
  28
  29 void loadMappingInfo(char* mappingF) {
  30   ifstream fin(mappingF);
  31   string line, key, value;
  32
  33   if (!fin.is_open()) {
  34           fprintf(stderr, "Cannot open %s! It may not exist.\n", mappingF);
  35           exit(-1);
  36   }
  37
  38   mi_table.clear();
  39   while (getline(fin, line)) {
  40     line = cleanStr(line);
  41     if (line[0] == '#') continue;
  42     istringstream strin(line);
  43     strin>>value>>key;
  44     mi_table[key] = value;
  45   }
  46
  47   fin.close();
  48 }
  49
  50 char check(char c) {
  51         if (!isalpha(c)) { fprintf(stderr, "Sequence contains unknown letter '%c'!\n", c); exit(-1); }
  52         //assert(isalpha(c));
  53         if (isupper(c) && c != 'A' && c != 'C' && c != 'G' && c != 'T') c = 'N';
  54         if (islower(c) && c != 'a' && c != 'c' && c != 'g' && c != 't') c = 'n';
  55         return c;
  56 }
  57
  58 void writeResults(char* refName) {
  59         ofstream fout;
  60         string cur_gene_id, name;
  61
  62         sprintf(groupF, "%s.grp", refName);
  63         sprintf(tiF, "%s.ti", refName);
  64         sprintf(refFastaF, "%s.transcripts.fa", refName);
  65
  66         transcripts.writeTo(tiF);
  67         if (verbose) { printf("Transcript Information File is generated!\n"); }
  68
  69         fout.open(groupF);
  70         cur_gene_id = "";
  71         for (int i = 1; i <= M; i++) {
  72                 const Transcript& transcript = transcripts.getTranscriptAt(i);
  73                 if (cur_gene_id != transcript.getGeneID()) {
  74                         fout<<i<<endl;
  75                         cur_gene_id = transcript.getGeneID();
  76                 }
  77         }
  78         fout<<M + 1<<endl;
  79         fout.close();
  80         if (verbose) { printf("Group File is generated!\n"); }
  81
  82         // We have to generate this .transcripts.fa, even reference_file is only one. Reason : polyA choice 2, need ">transcript_id"
  83         fout.open(refFastaF);
  84         for (int i = 1; i <= M; i++) {
  85                 name = transcripts.getTranscriptAt(i).getTranscriptID();
  86                 iter = name2seq.find(name);
  87                 if (iter == name2seq.end()) {
  88                         fprintf(stderr, "Cannot recognize transcript ID %s!\n", name.c_str());
  89                         exit(-1);
  90                 }
  91                 fout<<">"<<name<<endl;
  92                 fout<<iter->second<<endl;
  93         }
  94         fout.close();
  95         if (verbose) { printf("Extracted Sequences File is generated!\n"); }
  96 }
  97
  98 int main(int argc, char* argv[]) {
  99         if (argc < 5 || (hasMappingFile = atoi(argv[3])) && argc < 6) {
 100                 printf("Usage: synthesisRef refName quiet hasMappingFile [mappingFile] reference_file_1 [reference_file_2 ...]\n");
 101                 exit(-1);
 102         }
 103
 104         verbose = !atoi(argv[2]);
 105
 106         if (hasMappingFile) { loadMappingInfo(argv[4]); }
 107
 108         int start = hasMappingFile ? 5 : 4;
 109
 110         ifstream fin;
 111         string line, gseq;
 112         string seqname, gene_id;
 113         void* pt;
 114
 115         vector<Interval> vec;
 116
 117         M = 0;
 118         name2seq.clear();
 119         for (int i = start; i < argc; i++) {
 120                 fin.open(argv[i]);
 121                 if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", argv[i]); exit(-1); }
 122                 pt = getline(fin, line);
 123                 while (pt != 0 && line[0] == '>') {
 124                         istringstream strin(line.substr(1));
 125                         strin>>seqname;
 126
 127                         gseq = "";
 128                         while((pt = getline(fin, line)) && line[0] != '>') {
 129                             gseq += line;
 130                         }
 131
 132                         int len = gseq.length();
 133                         assert(len > 0);
 134                         for (int j = 0; j < len; j++) gseq[j] = check(gseq[j]);
 135
 136                         name2seq[seqname] = gseq;
 137
 138                         if (hasMappingFile) {
 139                               mi_iter = mi_table.find(seqname);
 140                               if (mi_iter == mi_table.end()) {
 141                                   fprintf(stderr, "Mapping Info is not correct, cannot find %s's gene_id!\n", seqname.c_str());
 142                                   exit(-1);
 143                               }
 144                               //assert(iter != table.end());
 145                               gene_id = mi_iter->second;
 146                         }
 147                         else gene_id = seqname;
 148
 149                         vec.clear();
 150                         vec.push_back(Interval(1, len));
 151                         transcripts.add(Transcript(seqname, gene_id, seqname, '+', vec, ""));
 152                         ++M;
 153
 154                         if (verbose && M % 1000000 == 0) { printf("%d sequences are processed!\n", M); }
 155                 }
 156                 fin.close();
 157         }
 158
 159         if (M < 1) {
 160                 fprintf(stderr, "Number of transcripts in the reference is less than 1!\n");
 161                 exit(-1);
 162         }
 163
 164         assert(M == transcripts.getM());
 165         transcripts.sort();
 166
 167         writeResults(argv[1]);
 168
 169         return 0;
 170 }