2 * transcripts are numbered from 1. 0 is reserved for noise isoform
16 #include "my_assert.h"
17 #include "Transcript.h"
22 Transcripts(int type = 0) {
23 M = 0; this->type = type;
25 transcripts.push_back(Transcript());
27 e2i.clear(); i2e.clear();
30 int getM() { return M; }
31 int getType() { return type; }
33 const Transcript& getTranscriptAt(int pos) {
34 assert(pos > 0 && pos <= M);
35 return transcripts[pos];
38 void add(const Transcript& transcript) {
39 transcripts.push_back(transcript);
44 std::sort(transcripts.begin(), transcripts.end());
47 void readFrom(const char*);
48 void writeTo(const char*);
51 int getInternalSid(int eid) {
52 assert(eid > 0 && eid <= M);
56 const Transcript& getTranscriptViaEid(int eid) {
57 return transcripts[getInternalSid(eid)];
60 void buildMappings(int, char**);
63 int M, type; // type 0 from genome , 1 standalone transcriptome
64 std::vector<Transcript> transcripts;
66 std::vector<int> e2i, i2e; // external sid to internal sid, internal sid to external sid
69 void Transcripts::readFrom(const char* inpF) {
71 std::ifstream fin(inpF);
73 if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
77 transcripts.assign(M + 1, Transcript());
78 for (int i = 1; i <= M; i++) {
79 transcripts[i].read(fin);
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);
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;
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)!");
100 for (int i = 1; i <= M; i++) {
101 const std::string& tid = transcripts[i].getTranscriptID();
102 iter = dict.find(tid);
103 general_assert(iter == dict.end(), tid + " appears more than once!");
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;
119 #endif /* TRANSCRIPTS_H_ */