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 "utils.h"
  17 #include "my_assert.h"
  18 #include "Transcript.h"
  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         // used in shrinking the transcripts
  33         void setM(int M) { this->M = M; transcripts.resize(M + 1); }
  34
  35         void move(int from, int to) {
  36           assert(from >= to);
  37           if (from > to) transcripts[to] = transcripts[from];
  38         }
  39
  40         int getType() { return type; }
  41         void setType(int type) { this->type = type; }
  42
  43         bool isAlleleSpecific() { return type == 2; }
  44
  45         const Transcript& getTranscriptAt(int pos) {
  46                 assert(pos > 0 && pos <= M);
  47                 return transcripts[pos];
  48         }
  49
  50         void add(const Transcript& transcript) {
  51                 transcripts.push_back(transcript);
  52                 ++M;
  53         }
  54
  55         void sort() {
  56                 std::sort(transcripts.begin(), transcripts.end());
  57         }
  58
  59         void readFrom(const char*);
  60         void writeTo(const char*);
  61
  62         //Eid: external sid
  63         int getInternalSid(int eid) {
  64                 assert(eid > 0 && eid <= M);
  65                 return e2i[eid];
  66         }
  67
  68         const Transcript& getTranscriptViaEid(int eid) {
  69                 return transcripts[getInternalSid(eid)];
  70         }
  71
  72         void buildMappings(int, char**, const char* = NULL);
  73
  74 private:
  75         int M, type; // type 0 from genome, 1 standalone transcriptome, 2 allele-specific
  76         std::vector<Transcript> transcripts;
  77
  78         std::vector<int> e2i, i2e; // external sid to internal sid, internal sid to external sid
  79 };
  80
  81 void Transcripts::readFrom(const char* inpF) {
  82         std::string line;
  83         std::ifstream fin(inpF);
  84
  85         if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
  86
  87         fin>>M>>type;
  88         getline(fin, line);
  89         transcripts.assign(M + 1, Transcript());
  90         for (int i = 1; i <= M; i++) {
  91                 transcripts[i].read(fin);
  92         }
  93         fin.close();
  94 }
  95
  96 void Transcripts::writeTo(const char* outF) {
  97         std::ofstream fout(outF);
  98         fout<<M<<" "<<type<<std::endl;
  99         for (int i = 1; i <= M; i++) {
 100                 transcripts[i].write(fout);
 101         }
 102         fout.close();
 103 }
 104
 105 void Transcripts::buildMappings(int n_targets, char** target_name, const char* imdName) {
 106         std::map<std::string, int> dict;
 107         std::map<std::string, int>::iterator iter;
 108         std::vector<bool> appeared;
 109
 110         general_assert(n_targets > 0, "The SAM/BAM file declares less than one reference sequence!");
 111         general_assert(n_targets <=  M, "The SAM/BAM file declares more reference sequences (" + itos(n_targets) + ") than RSEM knows (" + itos(M) + ")!");
 112         if (n_targets < M) printf("Warning: The SAM/BAM file declares less reference sequences (%d) than RSEM knows (%d)!\n", n_targets, M);
 113
 114         dict.clear();
 115         for (int i = 1; i <= M; i++) {
 116                 const std::string& tid = isAlleleSpecific() ? transcripts[i].getSeqName() : transcripts[i].getTranscriptID();
 117                 iter = dict.find(tid);
 118                 general_assert(iter == dict.end(), "RSEM's indices might be corrupted, " + tid + " appears more than once!");
 119                 dict[tid] = i;
 120         }
 121
 122         e2i.assign(M + 1, 0);
 123         i2e.assign(M + 1, 0);
 124         appeared.assign(M + 1, false);
 125         for (int i = 0; i < n_targets; i++) {
 126                 iter = dict.find(std::string(target_name[i]));
 127                 general_assert(iter != dict.end(), "RSEM can not recognize reference sequence name " + cstrtos(target_name[i]) + "!");
 128                 general_assert(iter->second > 0, "Reference sequence name " + cstrtos(target_name[i]) + " appears more than once in the SAM/BAM file!");
 129                 e2i[i + 1] = iter->second;
 130                 i2e[iter->second] = i + 1;
 131                 iter->second = -1;
 132                 appeared[e2i[i + 1]] = true;
 133         }
 134
 135         if (imdName != NULL) {
 136           char omitF[STRLEN];
 137           sprintf(omitF, "%s.omit", imdName);
 138           FILE *fo = fopen(omitF, "w");
 139           for (int i = 1; i <= M; i++)
 140             if (!appeared[i]) fprintf(fo, "%d\n", i);
 141           fclose(fo);
 142         }
 143 }
 144
 145 #endif /* TRANSCRIPTS_H_ */