Transcripts.h

   1 /*
   2  * transcripts are numbered from 1. 0 is reserved for noise isoform
   3  */
   4 #ifndef TRANSCRIPTS_H_
   5 #define TRANSCRIPTS_H_
   6
   7 #include<cstdio>
   8 #include<cstdlib>
   9 #include<cassert>
  10 #include<fstream>
  11 #include<vector>
  12 #include<algorithm>
  13 #include<map>
  14 #include<string>
  15
  16 #include "my_assert.h"
  17 #include "Transcript.h"
  18
  19
  20 class Transcripts {
  21 public:
  22         Transcripts(int type = 0) {
  23                 M = 0; this->type = type;
  24                 transcripts.clear();
  25                 transcripts.push_back(Transcript());
  26
  27                 e2i.clear(); i2e.clear();
  28         }
  29
  30         int getM() { return M; }
  31         int getType() { return type; }
  32
  33         const Transcript& getTranscriptAt(int pos) {
  34                 assert(pos > 0 && pos <= M);
  35                 return transcripts[pos];
  36         }
  37
  38         void add(const Transcript& transcript) {
  39                 transcripts.push_back(transcript);
  40                 ++M;
  41         }
  42
  43         void sort() {
  44                 std::sort(transcripts.begin(), transcripts.end());
  45         }
  46
  47         void readFrom(const char*);
  48         void writeTo(const char*);
  49
  50         //Eid: external sid
  51         int getInternalSid(int eid) {
  52                 assert(eid > 0 && eid <= M);
  53                 return e2i[eid];
  54         }
  55
  56         const Transcript& getTranscriptViaEid(int eid) {
  57                 return transcripts[getInternalSid(eid)];
  58         }
  59
  60         void buildMappings(int, char**);
  61
  62 private:
  63         int M, type; // type 0 from genome , 1 standalone transcriptome
  64         std::vector<Transcript> transcripts;
  65
  66         std::vector<int> e2i, i2e; // external sid to internal sid, internal sid to external sid
  67 };
  68
  69 void Transcripts::readFrom(const char* inpF) {
  70         std::string line;
  71         std::ifstream fin(inpF);
  72
  73         if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
  74
  75         fin>>M>>type;
  76         getline(fin, line);
  77         transcripts.assign(M + 1, Transcript());
  78         for (int i = 1; i <= M; i++) {
  79                 transcripts[i].read(fin);
  80         }
  81         fin.close();
  82 }
  83
  84 void Transcripts::writeTo(const char* outF) {
  85         std::ofstream fout(outF);
  86         fout<<M<<" "<<type<<std::endl;
  87         for (int i = 1; i <= M; i++) {
  88                 transcripts[i].write(fout);
  89         }
  90         fout.close();
  91 }
  92
  93 void Transcripts::buildMappings(int n_targets, char** target_name) {
  94         std::map<std::string, int> dict;
  95         std::map<std::string, int>::iterator iter;
  96
  97         general_assert(n_targets == M, "Number of transcripts does not match! Please check if the reads are aligned to a transcript set (instead of a genome)!");
  98
  99         dict.clear();
 100         for (int i = 1; i <= M; i++) {
 101                 const std::string& tid = transcripts[i].getTranscriptID();
 102                 iter = dict.find(tid);
 103                 assert(iter == dict.end());
 104                 dict[tid] = i;
 105         }
 106
 107         e2i.assign(M + 1, 0);
 108         i2e.assign(M + 1, 0);
 109         for (int i = 0; i < n_targets; i++) {
 110                 iter = dict.find(std::string(target_name[i]));
 111                 general_assert(iter != dict.end(), "RSEM can not recognize transcript " + cstrtos(target_name[i]) + "!");
 112                 general_assert(iter->second > 0, "Reference sequence name " + cstrtos(target_name[i]) + " is duplicated!");
 113                 e2i[i + 1] = iter->second;
 114                 i2e[iter->second] = i + 1;
 115                 iter->second = -1;
 116         }
 117 }
 118
 119 #endif /* TRANSCRIPTS_H_ */