Transcript.h

   1 #ifndef TRANSCRIPT_H_
   2 #define TRANSCRIPT_H_
   3
   4 #include<cstdio>
   5 #include<cstdlib>
   6 #include<cassert>
   7 #include<string>
   8 #include<vector>
   9 #include<fstream>
  10
  11 #include "utils.h"
  12
  13 struct Interval {
  14         int start, end;
  15
  16         Interval(int start, int end) {
  17                 this->start = start;
  18                 this->end = end;
  19         }
  20 };
  21
  22 class Transcript {
  23 public:
  24         Transcript() {
  25                 length = 0;
  26                 structure.clear();
  27                 strand = 0;
  28                 seqname = gene_id = transcript_id = "";
  29                 left = "";
  30         }
  31
  32         Transcript(const std::string& transcript_id, const std::string& gene_id, const std::string& seqname,
  33                         const char& strand, const std::vector<Interval>& structure, const std::string& left) {
  34                 this->structure = structure;
  35                 this->strand = strand;
  36                 this->seqname = seqname;
  37                 this->gene_id = gene_id;
  38                 this->transcript_id = transcript_id;
  39
  40                 //eliminate prefix spaces in string variable "left"
  41                 int pos = 0;
  42                 int len = left.length();
  43                 while (pos < len && left[pos] == ' ') ++pos;
  44                 this->left = left.substr(pos);
  45
  46                 length = 0;
  47                 int s = structure.size();
  48                 for (int i = 0; i < s; i++) length += structure[i].end + 1 - structure[i].start;
  49         }
  50
  51         bool operator< (const Transcript& o) const {
  52           return gene_id < o.gene_id || (gene_id == o.gene_id && transcript_id < o.transcript_id);
  53         }
  54
  55         const std::string& getTranscriptID() const { return transcript_id; }
  56
  57         const std::string& getGeneID() const { return gene_id; }
  58
  59         const std::string& getSeqName() const { return seqname; }
  60
  61         char getStrand() const { return strand; }
  62
  63         const std::string& getLeft() const { return left; }
  64
  65         int getLength() const { return length; }
  66
  67         const std::vector<Interval>& getStructure() const { return structure; }
  68
  69         void extractSeq (const std::string&, std::string&) const;
  70
  71         void read(std::ifstream&);
  72         void write(std::ofstream&);
  73
  74 private:
  75         int length; // transcript length
  76         std::vector<Interval> structure; // transcript structure , coordinate starts from 1
  77         char strand;
  78         std::string seqname, gene_id, transcript_id; // follow GTF definition
  79         std::string left;
  80 };
  81
  82 //gseq : genomic sequence
  83 void Transcript::extractSeq(const std::string& gseq, std::string& seq) const {
  84         seq = "";
  85         int s = structure.size();
  86         size_t glen = gseq.length();
  87
  88         if (structure[0].start < 1 || (size_t)structure[s - 1].end > glen) {
  89                 fprintf(stderr, "Transcript %s is out of chromosome %s's boundary!\n", transcript_id.c_str(), seqname.c_str());
  90                 exit(-1);
  91         }
  92
  93         switch(strand) {
  94         case '+':
  95                 for (int i = 0; i < s; i++) {
  96                         seq += gseq.substr(structure[i].start - 1, structure[i].end - structure[i].start + 1); // gseq starts from 0!
  97                 }
  98                 break;
  99         case '-':
 100                 for (int i = s - 1; i >= 0; i--) {
 101                         for (int j = structure[i].end; j >= structure[i].start; j--) {
 102                                 seq += getOpp(gseq[j - 1]);
 103                         }
 104                 }
 105                 break;
 106         default: assert(false);
 107         }
 108
 109         assert(seq.length() > 0);
 110 }
 111
 112 void Transcript::read(std::ifstream& fin) {
 113         int s;
 114         std::string tmp;
 115
 116         getline(fin, transcript_id);
 117         getline(fin, gene_id);
 118         getline(fin, seqname);
 119         fin>>tmp>>length;
 120         assert(tmp.length() == 1 && (tmp[0] == '+' || tmp[0] == '-'));
 121         strand = tmp[0];
 122         structure.clear();
 123         fin>>s;
 124         for (int i = 0; i < s; i++) {
 125                 int start, end;
 126                 fin>>start>>end;
 127                 structure.push_back(Interval(start, end));
 128         }
 129         getline(fin, tmp); //get the end of this line
 130         getline(fin, left);
 131 }
 132
 133 void Transcript::write(std::ofstream& fout) {
 134         int s = structure.size();
 135
 136         fout<<transcript_id<<std::endl;
 137         fout<<gene_id<<std::endl;
 138         fout<<seqname<<std::endl;
 139         fout<<strand<<" "<<length<<std::endl;
 140         fout<<s;
 141         for (int i = 0; i < s; i++) fout<<" "<<structure[i].start<<" "<<structure[i].end;
 142         fout<<std::endl;
 143         fout<<left<<std::endl;
 144 }
 145
 146 #endif /* TRANSCRIPT_H_ */