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(1); // no genome, just transcript set
  21 char groupF[STRLEN], tiF[STRLEN], refFastaF[STRLEN], chromListF[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, fout2;
  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         sprintf(chromListF, "%s.chrlist", refName);
  66
  67         transcripts.writeTo(tiF);
  68         if (verbose) { printf("Transcript Information File is generated!\n"); }
  69
  70         fout.open(groupF);
  71         cur_gene_id = "";
  72         for (int i = 1; i <= M; i++) {
  73                 const Transcript& transcript = transcripts.getTranscriptAt(i);
  74                 if (cur_gene_id != transcript.getGeneID()) {
  75                         fout<<i<<endl;
  76                         cur_gene_id = transcript.getGeneID();
  77                 }
  78         }
  79         fout<<M + 1<<endl;
  80         fout.close();
  81         if (verbose) { printf("Group File is generated!\n"); }
  82
  83         fout2.open(chromListF);
  84         fout.open(refFastaF);
  85         for (int i = 1; i <= M; i++) {
  86                 name = transcripts.getTranscriptAt(i).getTranscriptID();
  87                 iter = name2seq.find(name);
  88                 if (iter == name2seq.end()) {
  89                         fprintf(stderr, "Cannot recognize transcript ID %s!\n", name.c_str());
  90                         exit(-1);
  91                 }
  92                 fout<<">"<<name<<endl;
  93                 fout<<iter->second<<endl;
  94
  95                 fout2<<name<<'\t'<<iter->second.length()<<endl;
  96         }
  97         fout.close();
  98         fout2.close();
  99
 100         if (verbose) {
 101           printf("Chromosome List File is generated!\n");
 102           printf("Extracted Sequences File is generated!\n");
 103         }
 104 }
 105
 106 int main(int argc, char* argv[]) {
 107   if (argc < 5 || ((hasMappingFile = atoi(argv[3])) && argc < 6)) {
 108                 printf("Usage: synthesisRef refName quiet hasMappingFile [mappingFile] reference_file_1 [reference_file_2 ...]\n");
 109                 exit(-1);
 110         }
 111
 112         verbose = !atoi(argv[2]);
 113
 114         if (hasMappingFile) { loadMappingInfo(argv[4]); }
 115
 116         int start = hasMappingFile ? 5 : 4;
 117
 118         ifstream fin;
 119         string line, gseq;
 120         string seqname, gene_id;
 121         void* pt;
 122
 123         vector<Interval> vec;
 124
 125         M = 0;
 126         name2seq.clear();
 127         for (int i = start; i < argc; i++) {
 128                 fin.open(argv[i]);
 129                 if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", argv[i]); exit(-1); }
 130                 pt = getline(fin, line);
 131                 while (pt != 0 && line[0] == '>') {
 132                         istringstream strin(line.substr(1));
 133                         strin>>seqname;
 134
 135                         gseq = "";
 136                         while((pt = getline(fin, line)) && line[0] != '>') {
 137                             gseq += line;
 138                         }
 139
 140                         int len = gseq.length();
 141                         assert(len > 0);
 142                         for (int j = 0; j < len; j++) gseq[j] = check(gseq[j]);
 143
 144                         name2seq[seqname] = gseq;
 145
 146                         if (hasMappingFile) {
 147                               mi_iter = mi_table.find(seqname);
 148                               if (mi_iter == mi_table.end()) {
 149                                   fprintf(stderr, "Mapping Info is not correct, cannot find %s's gene_id!\n", seqname.c_str());
 150                                   exit(-1);
 151                               }
 152                               //assert(iter != table.end());
 153                               gene_id = mi_iter->second;
 154                         }
 155                         else gene_id = seqname;
 156
 157                         vec.clear();
 158                         vec.push_back(Interval(1, len));
 159                         transcripts.add(Transcript(seqname, gene_id, seqname, '+', vec, ""));
 160                         ++M;
 161
 162                         if (verbose && M % 1000000 == 0) { printf("%d sequences are processed!\n", M); }
 163                 }
 164                 fin.close();
 165         }
 166
 167         if (M < 1) {
 168                 fprintf(stderr, "Number of transcripts in the reference is less than 1!\n");
 169                 exit(-1);
 170         }
 171
 172         assert(M == transcripts.getM());
 173         transcripts.sort();
 174
 175         writeResults(argv[1]);
 176
 177         return 0;
 178 }