synthesisRef.cpp

   1 #include<cstdio>
   2 #include<cstring>
   3 #include<cstdlib>
   4 #include<cassert>
   5 #include<fstream>
   6 #include<sstream>
   7 #include<map>
   8 #include<vector>
   9
  10 #include "utils.h"
  11 #include "my_assert.h"
  12 #include "Transcript.h"
  13 #include "Transcripts.h"
  14
  15 using namespace std;
  16
  17 int M;
  18
  19 map<string, string> name2seq;
  20 map<string, string>::iterator iter;
  21
  22 Transcripts transcripts(1); // no genome, just transcript set
  23 char groupF[STRLEN], tiF[STRLEN], refFastaF[STRLEN], chromListF[STRLEN];
  24 char gtF[STRLEN], taF[STRLEN]; // group info between gene and transcript, transcript and allele
  25
  26 int hasMappingFile;
  27 char mappingFile[STRLEN];
  28
  29 map<string, string> mi_table; // mapping info table
  30 map<string, string>::iterator mi_iter; //mapping info table's iterator
  31
  32 map<string, string> mi_table2; // allele_id to transcript_id
  33 map<string, string>::iterator mi_iter2; // corresponding iterator
  34
  35 void loadMappingInfo(int file_type, char* mappingF) {
  36   ifstream fin(mappingF);
  37   string line, key, value, value2;
  38
  39   general_assert(fin.is_open(), "Cannot open " + cstrtos(mappingF) + "! It may not exist.");
  40
  41   switch(file_type) {
  42   case 1:
  43     mi_table.clear();
  44     while (getline(fin, line)) {
  45       line = cleanStr(line);
  46       if (line[0] == '#') continue;
  47       istringstream strin(line);
  48       strin>>value>>key;
  49       mi_table[key] = value;
  50     }
  51     break;
  52   case 2:
  53     mi_table.clear();
  54     mi_table2.clear();
  55     while (getline(fin, line)) {
  56       line = cleanStr(line);
  57       if (line[0] == '#') continue;
  58       istringstream strin(line);
  59       strin>> value>> value2>> key;
  60       mi_table[key] = value;
  61       mi_table2[key] = value2;
  62     }
  63     break;
  64   default: assert(false);
  65   }
  66
  67   fin.close();
  68 }
  69
  70 char check(char c) {
  71         if (!isalpha(c)) { fprintf(stderr, "Sequence contains unknown letter '%c'!\n", c); exit(-1); }
  72         //assert(isalpha(c));
  73         if (isupper(c) && c != 'A' && c != 'C' && c != 'G' && c != 'T') c = 'N';
  74         if (islower(c) && c != 'a' && c != 'c' && c != 'g' && c != 't') c = 'n';
  75         return c;
  76 }
  77
  78 void writeResults(int option, char* refName) {
  79         ofstream fout, fout2;
  80         string cur_gene_id, cur_transcript_id, name;
  81         vector<int> gi, gt, ta;
  82
  83         sprintf(tiF, "%s.ti", refName);
  84         transcripts.writeTo(tiF);
  85         if (verbose) { printf("Transcript Information File is generated!\n"); }
  86
  87         cur_gene_id = ""; gi.clear();
  88         if (option == 2) { cur_transcript_id = ""; gt.clear(); ta.clear(); }
  89         for (int i = 1; i <= M; i++) {
  90                 const Transcript& transcript = transcripts.getTranscriptAt(i);
  91                 if (cur_gene_id != transcript.getGeneID()) {
  92                   gi.push_back(i);
  93                   if (option == 2) gt.push_back((int)ta.size());
  94                   cur_gene_id = transcript.getGeneID();
  95                 }
  96                 if ((option == 2) && (cur_transcript_id != transcript.getTranscriptID())) {
  97                     ta.push_back(i);
  98                     cur_transcript_id = transcript.getTranscriptID();
  99                 }
 100         }
 101         gi.push_back(M + 1);
 102         if (option == 2) { gt.push_back((int)ta.size()); ta.push_back(M + 1); }
 103
 104         sprintf(groupF, "%s.grp", refName);
 105         fout.open(groupF);
 106         for (int i = 0; i < (int)gi.size(); i++) fout<< gi[i]<< endl;
 107         fout.close();
 108         if (verbose) { printf("Group File is generated!\n"); }
 109
 110         if (option == 2) {
 111           sprintf(gtF, "%s.gt", refName);
 112           fout.open(gtF);
 113           for (int i = 0; i < (int)gt.size(); i++) fout<< gt[i]<< endl;
 114           fout.close();
 115           sprintf(taF, "%s.ta", refName);
 116           fout.open(taF);
 117           for (int i = 0; i < (int)ta.size(); i++) fout<< ta[i]<< endl;
 118           fout.close();
 119           if (verbose) { printf("Allele-specific group files are generated!\n"); }
 120         }
 121
 122         sprintf(refFastaF, "%s.transcripts.fa", refName);
 123         sprintf(chromListF, "%s.chrlist", refName);
 124         fout2.open(chromListF);
 125         fout.open(refFastaF);
 126         for (int i = 1; i <= M; i++) {
 127                 name = transcripts.getTranscriptAt(i).getSeqName();
 128                 iter = name2seq.find(name);
 129                 general_assert(iter != name2seq.end(), "Cannot recognize sequence ID" + name + "!");
 130                 fout<<">"<<name<<endl;
 131                 fout<<iter->second<<endl;
 132
 133                 fout2<<name<<'\t'<<iter->second.length()<<endl;
 134         }
 135         fout.close();
 136         fout2.close();
 137
 138         if (verbose) {
 139           printf("Chromosome List File is generated!\n");
 140           printf("Extracted Sequences File is generated!\n");
 141         }
 142 }
 143
 144 int main(int argc, char* argv[]) {
 145   if (argc < 5 || ((hasMappingFile = atoi(argv[3])) && argc < 6)) {
 146                 printf("Usage: synthesisRef refName quiet hasMappingFile<0,no;1,yes;2,allele-specific> [mappingFile] reference_file_1 [reference_file_2 ...]\n");
 147                 exit(-1);
 148         }
 149
 150         verbose = !atoi(argv[2]);
 151
 152         if (hasMappingFile) { loadMappingInfo(hasMappingFile, argv[4]); }
 153
 154         // allele-specific
 155         if (hasMappingFile == 2) { transcripts.setType(2); }
 156
 157         int start = hasMappingFile ? 5 : 4;
 158
 159         ifstream fin;
 160         string line, gseq;
 161         string seqname, gene_id, transcript_id;
 162
 163         vector<Interval> vec;
 164
 165         M = 0;
 166         name2seq.clear();
 167         for (int i = start; i < argc; i++) {
 168                 fin.open(argv[i]);
 169                 general_assert(fin.is_open(), "Cannot open " + cstrtos(argv[i]) + "! It may not exist.");
 170                 getline(fin, line);
 171                 while ((fin) && (line[0] == '>')) {
 172                         istringstream strin(line.substr(1));
 173                         strin>>seqname;
 174
 175                         gseq = "";
 176                         while((getline(fin, line)) && (line[0] != '>')) {
 177                             gseq += line;
 178                         }
 179
 180                         int len = gseq.length();
 181                         assert(len > 0);
 182                         for (int j = 0; j < len; j++) gseq[j] = check(gseq[j]);
 183
 184                         name2seq[seqname] = gseq;
 185
 186                         transcript_id = seqname;
 187                         gene_id = seqname;
 188
 189                         if (hasMappingFile) {
 190                               mi_iter = mi_table.find(seqname);
 191                               general_assert(mi_iter != mi_table.end(), "Mapping Info is not correct, cannot find " + seqname + "'s gene_id!");
 192                               gene_id = mi_iter->second;
 193                               if (hasMappingFile == 2) {
 194                                 mi_iter2 = mi_table2.find(seqname);
 195                                 general_assert(mi_iter2 != mi_table2.end(), "Mapping Info is not correct, cannot find allele " + seqname + "'s transcript_id!");
 196                                 transcript_id = mi_iter2->second;
 197                               }
 198                         }
 199
 200                         vec.clear();
 201                         vec.push_back(Interval(1, len));
 202                         transcripts.add(Transcript(transcript_id, gene_id, seqname, '+', vec, ""));
 203                         ++M;
 204
 205                         if (verbose && M % 1000000 == 0) { printf("%d sequences are processed!\n", M); }
 206                 }
 207                 fin.close();
 208         }
 209
 210         if (M < 1) {
 211                 fprintf(stderr, "Number of transcripts in the reference is less than 1!\n");
 212                 exit(-1);
 213         }
 214
 215         assert(M == transcripts.getM());
 216         transcripts.sort();
 217
 218         writeResults(hasMappingFile, argv[1]);
 219
 220         return 0;
 221 }