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
  32         int getType() { return type; }
  33         void setType(int type) { this->type = type; }
  34
  35         bool isAlleleSpecific() { return type == 2; }
  36
  37         const Transcript& getTranscriptAt(int pos) {
  38                 assert(pos > 0 && pos <= M);
  39                 return transcripts[pos];
  40         }
  41
  42         void add(const Transcript& transcript) {
  43                 transcripts.push_back(transcript);
  44                 ++M;
  45         }
  46
  47         void sort() {
  48                 std::sort(transcripts.begin(), transcripts.end());
  49         }
  50
  51         void readFrom(const char*);
  52         void writeTo(const char*);
  53
  54         //Eid: external sid
  55         int getInternalSid(int eid) {
  56                 assert(eid > 0 && eid <= M);
  57                 return e2i[eid];
  58         }
  59
  60         const Transcript& getTranscriptViaEid(int eid) {
  61                 return transcripts[getInternalSid(eid)];
  62         }
  63
  64         void buildMappings(int, char**);
  65
  66 private:
  67         int M, type; // type 0 from genome, 1 standalone transcriptome, 2 allele-specific
  68         std::vector<Transcript> transcripts;
  69
  70         std::vector<int> e2i, i2e; // external sid to internal sid, internal sid to external sid
  71 };
  72
  73 void Transcripts::readFrom(const char* inpF) {
  74         std::string line;
  75         std::ifstream fin(inpF);
  76
  77         if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
  78
  79         fin>>M>>type;
  80         getline(fin, line);
  81         transcripts.assign(M + 1, Transcript());
  82         for (int i = 1; i <= M; i++) {
  83                 transcripts[i].read(fin);
  84         }
  85         fin.close();
  86 }
  87
  88 void Transcripts::writeTo(const char* outF) {
  89         std::ofstream fout(outF);
  90         fout<<M<<" "<<type<<std::endl;
  91         for (int i = 1; i <= M; i++) {
  92                 transcripts[i].write(fout);
  93         }
  94         fout.close();
  95 }
  96
  97 void Transcripts::buildMappings(int n_targets, char** target_name) {
  98         std::map<std::string, int> dict;
  99         std::map<std::string, int>::iterator iter;
 100
 101         general_assert(n_targets == M, "Number of reference sequences does not match! Please check if the reads are aligned to a transcript set (instead of a genome)!");
 102
 103         dict.clear();
 104         for (int i = 1; i <= M; i++) {
 105                 const std::string& tid = isAlleleSpecific() ? transcripts[i].getSeqName() : transcripts[i].getTranscriptID();
 106                 iter = dict.find(tid);
 107                 general_assert(iter == dict.end(), tid + " appears more than once!");
 108                 dict[tid] = i;
 109         }
 110
 111         e2i.assign(M + 1, 0);
 112         i2e.assign(M + 1, 0);
 113         for (int i = 0; i < n_targets; i++) {
 114                 iter = dict.find(std::string(target_name[i]));
 115                 general_assert(iter != dict.end(), "RSEM can not recognize reference sequence name " + cstrtos(target_name[i]) + "!");
 116                 general_assert(iter->second > 0, "Reference sequence name " + cstrtos(target_name[i]) + " is duplicated!");
 117                 e2i[i + 1] = iter->second;
 118                 i2e[iter->second] = i + 1;
 119                 iter->second = -1;
 120         }
 121 }
 122
 123 #endif /* TRANSCRIPTS_H_ */