sequence.cc

   1 //                        ----------------------------------------
   2 //                           ---------- sequence.cc -----------
   3 //                        ----------------------------------------
   4
   5 #include "sequence.hh"
   6
   7
   8 Sequence::Sequence()
   9 {
  10 }
  11
  12
  13 void
  14 Sequence::load_fasta(string file_path, int seq_num)
  15 {
  16   fstream data_file;
  17   string file_data_line;
  18   int header_counter = 0;
  19   bool read_seq = true;
  20
  21
  22   data_file.open(file_path.c_str(), ios::in);
  23
  24
  25   // search for the header of the fasta sequence we want
  26   while ( (!data_file.eof()) && (header_counter < seq_num) )
  27   {
  28     getline(data_file,file_data_line);
  29     if (file_data_line.substr(0,1) == ">")
  30       header_counter++;
  31   }
  32
  33   header = file_data_line.substr(1);
  34
  35   while ( !data_file.eof() && read_seq )
  36   {
  37     getline(data_file,file_data_line);
  38     if (file_data_line.substr(0,1) == ">")
  39       read_seq = false;
  40     else sequence += file_data_line;
  41   }
  42
  43   data_file.close();
  44
  45   // need more sequence filtering for non AGCTN and convert to N
  46   transform(sequence.begin(), sequence.end(), sequence.begin(), toupper);
  47 }
  48
  49
  50 void
  51 Sequence::load_annot(string file_path)
  52 {
  53   fstream data_file;
  54   string file_data_line;
  55   annot an_annot;
  56
  57   data_file.open(file_path.c_str(), ios::in);
  58
  59   data_file.close();
  60
  61 }
  62
  63
  64 int
  65 Sequence::len()
  66 {
  67   return sequence.length();
  68 }
  69
  70
  71 string
  72 Sequence::seq()
  73 {
  74   return sequence;
  75 }
  76
  77 const char *
  78 Sequence::c_seq()
  79 {
  80   return sequence.c_str();
  81 }
  82
  83 string
  84 Sequence::rev_comp()
  85 {
  86   string rev_comp;
  87   char conversionTable[257];
  88   int seq_i, table_i, len;
  89
  90   len = sequence.length();
  91   rev_comp.reserve(len);
  92   // make a conversion table
  93   // init all parts of conversion table to '~' character
  94   // '~' I doubt will ever appear in a sequence file (jeez, I hope)
  95   // and may the fleas of 1000 camels infest the genitals of any biologist (and
  96   // seven generations of their progeny) who decides to make it mean
  97   // something special!!!
  98   // PS - double the curse for any smartass non-biologist who tries it as well
  99   for(table_i=0; table_i < 256; table_i++)
 100   {
 101     conversionTable[table_i] = '~';
 102   }
 103   // add end of string character for printing out table for testing purposes
 104   conversionTable[256] = '\0';
 105
 106   // add in the characters for the bases we want to convert
 107   conversionTable['A'] = 'T';
 108   conversionTable['T'] = 'A';
 109   conversionTable['G'] = 'C';
 110   conversionTable['C'] = 'G';
 111   conversionTable['N'] = 'N';
 112
 113   // finally, the actual conversion loop
 114   for(seq_i = len - 1; seq_i >= 0; seq_i--)
 115   {
 116     table_i = (int) sequence[seq_i];
 117     rev_comp += conversionTable[table_i];
 118   }
 119
 120   return rev_comp;
 121 }
 122
 123
 124 string
 125 Sequence::hdr()
 126 {
 127   return header;
 128 }
 129
 130
 131 void
 132 Sequence::set_seq(string a_seq)
 133 {
 134   sequence = a_seq;
 135 }
 136
 137
 138 /*
 139 string
 140 Sequence::species()
 141 {
 142   return species;
 143 }
 144 */
 145
 146 void
 147 Sequence::clear()
 148 {
 149   sequence = "";
 150   length = 0;
 151   header = "";
 152   species = "";
 153   species_num = -4;
 154   annots.clear();
 155 }