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