X-Git-Url: https://git.donarmstrong.com/?p=rsem.git;a=blobdiff_plain;f=Transcripts.h;h=db8e64a879d50bdf1e2e11e8df1037a1e62c5467;hp=2750e0c36394d3c3be5cb7d63dfb74f7ae5687a3;hb=92b24279a3ecc72946e7e7c23149ad0d181f373a;hpb=a95154919f950f86de9104b2b9dcf1f0c7e83387

diff --git a/Transcripts.h b/Transcripts.h
index 2750e0c..db8e64a 100644
--- a/Transcripts.h
+++ b/Transcripts.h
@@ -10,20 +10,37 @@
 #include<fstream>
 #include<vector>
 #include<algorithm>
+#include<map>
+#include<string>
 
+#include "utils.h"
+#include "my_assert.h"
 #include "Transcript.h"
 
-
 class Transcripts {
 public:
 	Transcripts(int type = 0) {
 		M = 0; this->type = type;
 		transcripts.clear();
 		transcripts.push_back(Transcript());
+
+		e2i.clear(); i2e.clear();
 	}
 
 	int getM() { return M; }
+
+	// used in shrinking the transcripts
+	void setM(int M) { this->M = M; transcripts.resize(M + 1); } 
+	
+	void move(int from, int to) {
+	  assert(from >= to);
+	  if (from > to) transcripts[to] = transcripts[from];
+	}
+	
 	int getType() { return type; }
+	void setType(int type) { this->type = type; }
+
+	bool isAlleleSpecific() { return type == 2; }
 
 	const Transcript& getTranscriptAt(int pos) {
 		assert(pos > 0 && pos <= M);
@@ -42,9 +59,23 @@ public:
 	void readFrom(const char*);
 	void writeTo(const char*);
 
+	//Eid: external sid
+	int getInternalSid(int eid) {
+		assert(eid > 0 && eid <= M);
+		return e2i[eid];
+	}
+
+	const Transcript& getTranscriptViaEid(int eid) {
+		return transcripts[getInternalSid(eid)];
+	}
+
+	void buildMappings(int, char**, const char* = NULL);
+
 private:
-	int M, type; // type 0 from genome , 1 standalone transcriptome
+	int M, type; // type 0 from genome, 1 standalone transcriptome, 2 allele-specific 
 	std::vector<Transcript> transcripts;
+
+	std::vector<int> e2i, i2e; // external sid to internal sid, internal sid to external sid
 };
 
 void Transcripts::readFrom(const char* inpF) {
@@ -55,8 +86,7 @@ void Transcripts::readFrom(const char* inpF) {
 
 	fin>>M>>type;
 	getline(fin, line);
-	transcripts.clear();
-	transcripts.resize(M + 1);
+	transcripts.assign(M + 1, Transcript());
 	for (int i = 1; i <= M; i++) {
 		transcripts[i].read(fin);
 	}
@@ -72,4 +102,44 @@ void Transcripts::writeTo(const char* outF) {
 	fout.close();
 }
 
+void Transcripts::buildMappings(int n_targets, char** target_name, const char* imdName) {
+	std::map<std::string, int> dict;
+	std::map<std::string, int>::iterator iter;
+	std::vector<bool> appeared;
+
+	general_assert(n_targets > 0, "The SAM/BAM file declares less than one reference sequence!");
+	general_assert(n_targets <=  M, "The SAM/BAM file declares more reference sequences (" + itos(n_targets) + ") than RSEM knows (" + itos(M) + ")!");
+	if (n_targets < M) printf("Warning: The SAM/BAM file declares less reference sequences (%d) than RSEM knows (%d)!\n", n_targets, M);
+
+	dict.clear();
+	for (int i = 1; i <= M; i++) {
+		const std::string& tid = isAlleleSpecific() ? transcripts[i].getSeqName() : transcripts[i].getTranscriptID();
+		iter = dict.find(tid);
+		general_assert(iter == dict.end(), "RSEM's indices might be corrupted, " + tid + " appears more than once!");
+		dict[tid] = i;
+	}
+
+	e2i.assign(M + 1, 0);
+	i2e.assign(M + 1, 0);
+	appeared.assign(M + 1, false);
+	for (int i = 0; i < n_targets; i++) {
+		iter = dict.find(std::string(target_name[i]));
+		general_assert(iter != dict.end(), "RSEM can not recognize reference sequence name " + cstrtos(target_name[i]) + "!");
+		general_assert(iter->second > 0, "Reference sequence name " + cstrtos(target_name[i]) + " appears more than once in the SAM/BAM file!");
+		e2i[i + 1] = iter->second;
+		i2e[iter->second] = i + 1;
+		iter->second = -1;
+		appeared[e2i[i + 1]] = true;
+	}
+
+	if (imdName != NULL) {
+	  char omitF[STRLEN];
+	  sprintf(omitF, "%s.omit", imdName);
+	  FILE *fo = fopen(omitF, "w");
+	  for (int i = 1; i <= M; i++) 
+	    if (!appeared[i]) fprintf(fo, "%d\n", i);
+	  fclose(fo);
+	}
+}
+
 #endif /* TRANSCRIPTS_H_ */