--- /dev/null
+#include "alg/color.h"
+#include <math.h>
+
+Color::Color()
+{
+ set(0.0, 0.0, 0.0, 0.0);
+}
+
+Color::Color(const Color &c)
+{
+ set(c.colors[RedChannel],
+ c.colors[GreenChannel],
+ c.colors[BlueChannel],
+ c.colors[AlphaChannel]);
+}
+
+Color::Color(float red, float green, float blue, float alpha)
+{
+ set(red, green, blue, alpha);
+}
+
+void Color::set(float r, float g, float b, float a)
+{
+ colors[RedChannel] = r;
+ colors[GreenChannel] = g;
+ colors[BlueChannel] = b;
+ colors[AlphaChannel] = a;
+}
+
+const float *const Color::get() const
+{
+ return colors;
+}
+
+void Color::setR(float r)
+{
+ colors[RedChannel] = r;
+}
+
+float Color::r() const
+{
+ return colors[RedChannel];
+}
+
+void Color::setG(float g)
+{
+ colors[GreenChannel] = g;
+}
+
+float Color::g() const
+{
+ return colors[GreenChannel];
+}
+
+void Color::setB(float b)
+{
+ colors[BlueChannel] = b;
+}
+
+float Color::b() const
+{
+ return colors[BlueChannel];
+}
+
+void Color::setA(float a)
+{
+ colors[AlphaChannel] = a;
+}
+
+float Color::a() const
+{
+ return colors[AlphaChannel];
+}
+
+static bool close(float l, float r, float t=1e-6)
+{
+ float difference = fabsf(l-r);
+ l = fabsf(l);
+ r = fabsf(r);
+
+ return (difference <= (t*l) or difference <= (t*r));
+}
+bool operator==(const Color &x, const Color &y)
+{
+ return (close(x.colors[Color::RedChannel] ,y.colors[Color::RedChannel]) and
+ close(x.colors[Color::GreenChannel],y.colors[Color::GreenChannel]) and
+ close(x.colors[Color::BlueChannel] ,y.colors[Color::BlueChannel]) and
+ close(x.colors[Color::AlphaChannel],y.colors[Color::AlphaChannel]));
+}
+
+std::ostream &operator<<(std::ostream &out, const Color &c)
+{
+ out << "Color(" << c.r() << ", "
+ << c.g() << ", "
+ << c.b() << ", "
+ << c.a() << ")";
+}
--- /dev/null
+#ifndef _GLCOLOR_H_
+#define _GLCOLOR_H_
+
+#include <ostream>
+
+//! convienece class for handling opengl colors
+class Color
+{
+public:
+ Color();
+ Color(const Color &);
+ //! initialize with red, green, blue, alpha
+ Color(float r, float g, float b, float a=0.0);
+
+ //! set all channels simultaneously
+ void set(float r, float g, float b, float a=0.0);
+ //! return an array of 4 colors (stored by class)
+ const float* const get() const;
+ void setR(float);
+ float r() const;
+ void setG(float);
+ float g() const;
+ void setB(float);
+ float b() const;
+ void setA(float);
+ float a() const;
+
+ enum color_channels { RedChannel, GreenChannel, BlueChannel, AlphaChannel };
+ friend bool operator==(const Color&, const Color&);
+ friend std::ostream& operator<<(std::ostream&, const Color&);
+
+protected:
+ float colors[4];
+};
+
+#endif
+
//cout << all_matches->size() << "\n";
}
}
- cout << "windows in flp = " << all_matches->size() << endl;
+ //cout << "windows in flp = " << all_matches->size() << endl;
data_file.close();
}
} // end of loop thru the current offset sequence
} // end of loop thru the start positions of seq2 sequence
}
-
-/*
- cout << "fee\n";
- cout << "fie\n";
- cout << "foe\n";
- cout << "fum\n";
-*/
seq_x(0.0),
seq_y(0.0),
seq_z(1.0),
- seq_width(s.size()),
- seq_height(gl_track_height)
+ seq_height(gl_track_height),
+ drawColor(0.0, 0.0, 0.0),
+ char_pix_per_world_unit(5.0)
{
- assert (seq.size() == seq_width);
}
GlSequence::GlSequence(const GlSequence &s)
seq_x(s.seq_x),
seq_y(s.seq_y),
seq_z(s.seq_z),
- seq_width(s.seq_width),
- seq_height(s.seq_height)
+ seq_height(s.seq_height),
+ drawColor(s.drawColor),
+ char_pix_per_world_unit(s.char_pix_per_world_unit)
{
- assert (seq.size() == seq_width);
}
GlSequence &GlSequence::operator=(const GlSequence & s)
seq_x = s.seq_x;
seq_y = s.seq_y;
seq_z = s.seq_z;
- seq_width = s.seq_width;
seq_height = s.seq_height;
+ drawColor = s.drawColor;
+ assert(char_pix_per_world_unit == s.char_pix_per_world_unit);
}
return *this;
}
return seq_y;
}
-void GlSequence::setWidth(GLfloat value)
+GLfloat GlSequence::length() const
{
- seq_width = value;
+ return seq.size();
}
-GLfloat GlSequence::width() const
+Sequence::size_type GlSequence::leftbase(GLfloat left) const
{
- return seq_width;
+ assert (seq_x == 0);
+ left = ceil(left);
+ if (left < seq_x)
+ return 0;
+ else
+ return (Sequence::size_type)left;
+}
+
+Sequence::size_type GlSequence::rightbase(GLfloat right) const
+{
+ assert (seq_x == 0);
+ right = floor(right);
+ if (right > seq.size())
+ return seq.size();
+ else
+ return (Sequence::size_type)right;
+}
+
+Sequence::const_iterator GlSequence::sequence_begin() const
+{
+ return seq.begin();
+}
+
+Sequence::const_iterator GlSequence::sequence_end() const
+{
+ return seq.end();
+}
+
+Sequence::const_iterator
+GlSequence::sequence_begin(GLfloat left, GLfloat right) const
+{
+ // the following code will be wrong when sequences can be slid around
+ // so make sure we break.
+ assert (seq_x == 0);
+
+ if ( leftbase(left) > seq.size() or left > right )
+ return seq.end();
+ else
+ return seq.begin() + leftbase(left);
+}
+
+Sequence::const_iterator
+GlSequence::sequence_end(GLfloat left, GLfloat right) const
+{
+ // the following code will be wrong when sequences can be slid around
+ // so make sure we break.
+ assert (seq_x == 0);
+
+ if ( rightbase(right) > seq.size() or left > right )
+ return seq.end();
+ else
+ return seq.begin() + rightbase(right);
+}
+
+
+//! set default track draw color
+void GlSequence::setColor(Color &c)
+{
+ drawColor = c;
+}
+
+//! get default track draw color
+Color GlSequence::color()
+{
+ return drawColor;
+}
+
+
+int GlSequence::get_viewport_pixel_width()
+{
+ int viewport[4];
+ glGetIntegerv(GL_VIEWPORT, viewport);
+ return viewport[3]; // grab the viewport width
+}
+
+bool GlSequence::is_sequence_renderable(GLfloat left,
+ GLfloat right,
+ int viewport_width) const
+{
+ // if called with default argument, go get the viewable width
+ if (viewport_width == -1) {
+ viewport_width = get_viewport_pixel_width();
+ }
+ GLfloat world_width = right - left;
+ GLfloat pixels_needed = (char_pix_per_world_unit * world_width);
+
+ // if the number of pixels taken up by rendering the characters
+ // that'd show up in the current ortho width is less than the window
+ // width we can actually draw something
+ return pixels_needed < viewport_width;
}
void GlSequence::draw(GLfloat left, GLfloat right) const
+{
+ if ( not is_sequence_renderable(left, right) ) {
+ draw_bar(left, right);
+ } else {
+ draw_sequence(left, right);
+ }
+ draw_annotations(left, right);
+}
+
+void GlSequence::draw_bar(GLfloat left, GLfloat right) const
{
glLineWidth(seq_height);
- glColor3f(0.0, 0.0, 0.0);
+ glColor3fv(drawColor.get());
// draw main sequence track
glBegin(GL_LINES);
- glVertex3f(seq_x, seq_y, seq_z);
- glVertex3f(seq_x+seq_width, seq_y, seq_z);
- //clog << "drawing " << seq_x << " " << seq_y << " " << seq_width
+ glVertex3f(seq_x, seq_y, seq_z);
+ glVertex3f(seq_x+seq.size(), seq_y, seq_z);
+ //clog << "drawing " << seq_x << " " << seq_y << " " << seq.size()
// << std::endl;
glEnd();
+}
+
+void GlSequence::draw_annotations(GLfloat left, GLfloat right) const
+{
// draw annotations
GLfloat annotation_z = seq_z + 1.0;
std::list<annot> annots = seq.annotations();
glVertex3f(annot_itor->end , seq_y, annotation_z);
glEnd();
}
- draw_sequence(left, right);
}
const int PT = 1;
void GlSequence::draw_sequence(GLfloat left, GLfloat right) const
{
- int viewport[4];
- glGetIntegerv(GL_VIEWPORT, viewport);
- int port_width = viewport[3]; // grab the viewport width
- GLfloat world_width = right - left;
- int left_base = (left > 0) ? (int)ceil(left) : 0;
- int right_base = (int)floor(right);
- int render_count = (right_base < seq.size()) ? right_base : seq.size();
- render_count -= left_base;
-
- glLineWidth(0.01);
- GLfloat char_width = 1.1;
- GLfloat pixels_needed = char_width * world_width;
-
- cout << "seq width needed " << pixels_needed
- << " port width " << port_width
- << " count " << render_count
- << " left " << left_base
- << " right " << right_base
- << " size " << seq.size() << endl;
- // if the number of pixels taken up by rendering the characters
- // that'd show up in the current ortho width is less than the window
- // width we can actually draw something
- if (pixels_needed < port_width and render_count > 0) {
- cout << "rendering: ";
- string bases = seq.subseq(left_base, render_count);
- glColor3f(0.1, 0.1, 0.1);
- for (string::size_type basepair = 0; basepair != bases.size(); ++basepair)
- {
- glPushMatrix();
- glTranslatef( left_base + basepair, seq_y+20, 1.0 );
- glScalef(0.1, 1.0, 1.0);
- switch (bases[basepair]) {
- case 'A': case 'a':
- drawLetter(Adata);
+ glLineWidth(0.1);
+ glColor3fv(drawColor.get());
+
+ Sequence::const_iterator seq_itor = sequence_begin(left, right);
+ Sequence::const_iterator seq_end = sequence_end(left, right);
+ Sequence::size_type basepair = 0;
+
+ assert(seq_end - seq_itor >= 0);
+ while(seq_itor != seq_end)
+ {
+ assert ( basepair < seq.size() );
+ glPushMatrix();
+ glTranslatef( leftbase(left) + basepair, seq_y, 1.0 );
+ glScalef(0.1, 1.0, 1.0);
+ switch (*seq_itor) {
+ case 'A': case 'a':
+ drawLetter(Adata);
break;
- case 'T': case 't':
- drawLetter(Tdata);
+ case 'T': case 't':
+ drawLetter(Tdata);
break;
- case 'G': case 'g':
- drawLetter(Gdata);
+ case 'G': case 'g':
+ drawLetter(Gdata);
break;
- case 'C': case 'c':
- drawLetter(Cdata);
+ case 'C': case 'c':
+ drawLetter(Cdata);
break;
- case 'N': case 'n':
- drawLetter(Ndata);
+ case 'N': case 'n':
+ drawLetter(Ndata);
break;
- default:
- drawLetter(Xdata);
+ default:
+ drawLetter(Xdata);
break;
- }
- cout << bases[basepair];
- glPopMatrix();
}
- cout << endl;
+ glPopMatrix();
+ ++seq_itor;
+ ++basepair;
}
}
+
+bool operator==(const GlSequence &left, const GlSequence &right)
+{
+ return ( (left.seq_x == right.seq_x) and
+ (left.seq_y == right.seq_y) and
+ (left.seq_z == right.seq_z) and
+ (left.seq_height == right.seq_height) and
+ (left.drawColor == right.drawColor));
+}
+
#define _GL_SEQUENCE_H_
#include "alg/sequence.hh"
+#include "alg/color.h"
#include <GL/gl.h>
-
//! Manage rendering a mussa sequence track
/*! The idea is this will keep track of the location of where the sequence
* is being rendered, and handle displaying annotations on that track
GlSequence &operator=(const GlSequence &s);
//! draw a track
- /*! left and right are the current edges of the viewable world
+ /*! left and right are the edges of the current viewport
*/
void draw(GLfloat left, GLfloat right) const;
const Sequence& sequence() const;
+ //! set our starting x (horizontal) coordinate
void setX(GLfloat);
+ //! get our starting x (horizontal) coordinate
GLfloat x() const;
+ //! set our current y (vertical) position
void setY(GLfloat);
+ //! get our current y (vertical) position
GLfloat y() const;
- void setWidth(GLfloat);
- GLfloat width() const;
+ //! how long is our sequence track? (computed from the sequence)
+ GLfloat length() const;
+
+ //! return the left (lowest) base index that is fully visible
+ Sequence::size_type GlSequence::leftbase(GLfloat left) const;
+ //! return one past rightmost (highest ) base index that is fully visible
+ /*! done mostly so all the iterator logic continues to work correctly.
+ */
+ Sequence::size_type GlSequence::rightbase(GLfloat right) const;
+
+ //! return iterator to the start of the stored sequence
+ Sequence::const_iterator sequence_begin() const;
+ //! return iterator to the end of the stored sequence
+ Sequence::const_iterator sequence_end() const;
+ //! provide an iterator to the sequence starting at world coordinate left
+ Sequence::const_iterator sequence_begin(GLfloat left, GLfloat right) const;
+ //! provide an iterator to the sequence ending at world coordinate right
+ Sequence::const_iterator sequence_end(GLfloat left, GLfloat right) const;
+
+ //! set track color
+ void setColor(Color &);
+ Color color();
+
+ //! are we close enough that it would make sense to view the base pairs?
+ /*! though we don't actually check to see if there's sequence in our
+ * view, just that there's enough pixels to render something if
+ * there were.
+ * \param[in] left the left edge of the viewable region in world coordinates
+ * \param[in] right the right edge of the viewable region in world
+ * coordinates
+ * \param[in] pixel_width allow setting the current viewport pixel width
+ */
+ bool is_sequence_renderable(GLfloat left,
+ GLfloat right,
+ int pixel_width=-1) const;
+
+ friend bool operator==(const GlSequence &left, const GlSequence &right);
protected:
const Sequence& seq;
GLfloat seq_x;
GLfloat seq_y;
GLfloat seq_z;
- GLfloat seq_width;
GLfloat seq_height;
+ Color drawColor;
+ const GLfloat char_pix_per_world_unit;
+ //! Return the pixel width of the opengl viewport.
+ static int get_viewport_pixel_width();
+ void draw_bar(GLfloat, GLfloat) const;
+ void draw_annotations(GLfloat, GLfloat) const;
//! render a sequence (if we have enough space
/*! left and right are the current edges of the viewable world
*/
};
const float gl_track_height = 10.0;
-
#endif
CURDIR := $(BASEDIR)alg/
-SOURCES.cxx := conserved_path.cxx \
+SOURCES.cxx := color.cxx \
+ conserved_path.cxx \
flp.cxx \
flp_seqcomp.cxx \
mussa_class.cxx \
else if (param == "SEQUENCE")
{
seq_files.push_back(file_path_base + value);
- cout << "seq_file_name " << seq_files.back() << endl;
+ //cout << "seq_file_name " << seq_files.back() << endl;
fasta_index = 1;
annot_file = "";
sub_seq_start = 0;
sub_seq_start = atoi(value.c_str());
else if (param == "SEQ_END")
{
- cout << "hey! " << atoi(value.c_str()) << endl;
sub_seq_end = atoi(value.c_str());
}
//ignore empty lines or that start with '#'
else if ((param == "") || (param == "#")) {}
else
{
- cout << "Illegal/misplaced mussa parameter in file\n";
- cout << param << "\n";
+ clog << "Illegal/misplaced mussa parameter in file\n";
+ clog << param << "\n";
}
if (!did_seq)
para_file.close();
soft_thres = threshold;
- cout << "nway mupa: analysis_name = " << analysis_name
- << " window = " << window
- << " threshold = " << threshold << endl;
+ //cout << "nway mupa: analysis_name = " << analysis_name
+ // << " window = " << window
+ // << " threshold = " << threshold << endl;
}
// no file was loaded, signal error
else
throw mussa_analysis_error("you need to have at least 2 sequences to "
"do an analysis.");
}
- cout << "nway ana: seq_num = " << the_seqs.size() << endl;
+ //cout << "nway ana: seq_num = " << the_seqs.size() << endl;
t2 = time(NULL);
seqloadtime = difftime(t2, t1);
totaltime = difftime(end, begin);
- cout << "seqload\tseqcomp\tnway\tsave\ttotal\n";
- //setup\t
- //cout << setuptime << "\t";
- cout << seqloadtime << "\t";
- cout << seqcomptime << "\t";
- cout << nwaytime << "\t";
- cout << savetime << "\t";
- cout << totaltime << "\n";
+ //cout << "seqload\tseqcomp\tnway\tsave\ttotal\n";
+ //cout << seqloadtime << "\t";
+ //cout << seqcomptime << "\t";
+ //cout << nwaytime << "\t";
+ //cout << savetime << "\t";
+ //cout << totaltime << "\n";
}
aSeq.load_annot(*annot_files_i, *seq_starts_i, *seq_ends_i);
the_seqs.push_back(aSeq);
- cout << aSeq.get_header() << endl;
+ //cout << aSeq.get_header() << endl;
//cout << aSeq.get_seq() << endl;
aSeq.clear();
++seq_files_i; // advance all the iterators
for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++)
for(vector<Sequence>::size_type i2 = i+1; i2 < the_seqs.size(); i2++)
{
- cout << "seqcomping: " << i << " v. " << i2 << endl;
+ //cout << "seqcomping: " << i << " v. " << i2 << endl;
all_comps[i][i2].setup(window, threshold);
all_comps[i][i2].seqcomp(the_seqs[i].get_seq(), the_seqs[i2].get_seq(), false);
all_comps[i][i2].seqcomp(the_seqs[i].get_seq(),the_seqs[i2].rev_comp(),true);
//create_dir_cmd = "mkdir " + file_path_base + save_name; //osX
dir_create_status = system( (const char*) create_dir_cmd.c_str());
- cout << "action: " << dir_create_status << endl;
+ //cout << "action: " << dir_create_status << endl;
// save sequence and annots to a special mussa file
vector<FLPs> empty_FLP_vector;
FLPs dummy_comp;
- cout << "ana_file name " << ana_file << endl;
+ //cout << "ana_file name " << ana_file << endl;
ana_path = ana_file;
parsing_path = true;
end_index = ana_path.size()-1;
++start_index;
}
analysis_name = ana_path.substr(start_index, end_index-start_index+1);
- cout << " ana_name " << analysis_name << endl;
+ //cout << " ana_name " << analysis_name << endl;
file_path_base = ana_path.substr(0, start_index) + analysis_name
+ "/" + analysis_name;
a_file_path = file_path_base + ".muway";
- cout << " loading museq: " << a_file_path << endl;
+ //cout << " loading museq: " << a_file_path << endl;
the_paths.load(a_file_path);
int seq_num = the_paths.sequence_count();
for (i = 1; i <= seq_num; i++)
{
tmp_seq.clear();
- cout << "mussa_class: loading museq frag... " << a_file_path << endl;
+ //cout << "mussa_class: loading museq frag... " << a_file_path << endl;
tmp_seq.load_museq(a_file_path, i);
the_seqs.push_back(tmp_seq);
}
{
append_info.str("");
append_info << "_sp_" << i << "v" << i2;
- cout << append_info.str() << endl;
+ //cout << append_info.str() << endl;
a_file_path = file_path_base + append_info.str() + ".flp";
all_comps[i][i2].load(a_file_path);
- cout << "real size = " << all_comps[i][i2].size() << endl;
+ //cout << "real size = " << all_comps[i][i2].size() << endl;
}
}
}
int sp_depth = all_matches.size() - 1; // this roves up & down species list
while ( (not_advanced) && (sp_depth != -1) )
{
- //cout << "foe\n";
//cout << sp_depth << ".." << *(sp_itor_begin[sp_depth]) << endl;
(sp_itor_begin[sp_depth])++; // advance iter at current depth
//cout << sp_depth << ".." << *(sp_itor_begin[sp_depth]) << endl;
// if we reached the end of the list, reset iter & go up one depth
if (sp_itor_begin[sp_depth] == sp_itor_end[sp_depth])
{
- //cout << "fum\n";
sp_itor_begin[sp_depth] = all_matches[sp_depth].begin();
sp_depth--;
//cout << "depth = " << sp_depth << endl;
pathz.clear();
window_num = all_comparisons[0][1].size();
- cout << window_num << endl; // just a sanity check
// loop thru all windows in first species
for (win_i = 0; win_i < window_num; win_i++)
reset_species_iterators(all_matches, sp_itor_begin, sp_itor_end);
still_paths = true;
- //cout << "fie\n";
while (still_paths)
{
// add path that each species iterator is pointing to
vector<list<int>::iterator> sp_itor_end(all_comparisons.size());
- cout << "trans: softhres = " << soft_thres;
- cout << ", window = " << win_size << ", ";
+ //cout << "trans: softhres = " << soft_thres;
+ //cout << ", window = " << win_size << ", ";
pathz.clear();
window_num = all_comparisons[0][1].size();
- cout << "window number = " << window_num << endl; // just a sanity check
- cout << "trans: comparison size= " << all_comparisons.size() << endl;
+ //cout << "window number = " << window_num << endl; // just a sanity check
+ //cout << "trans: comparison size= " << all_comparisons.size() << endl;
// loop thru all windows in first species
for (win_i = 0; win_i != window_num; win_i++)
{
}
}
}
- clog << "all_cmp=" << all_comparisons.size();
- if (pathz.begin() != pathz.end())
- clog << " path_size=" << pathz.begin()->size();
- else
- clog << " path empty";
- clog << endl;
+ //clog << "all_cmp=" << all_comparisons.size();
+ //if (pathz.begin() != pathz.end())
+ // clog << " path_size=" << pathz.begin()->size();
+ //else
+ // clog << " path empty";
+ //clog << endl;
assert ( (pathz.begin() == pathz.end())
|| (pathz.begin()->size() == 0)
|| (all_comparisons.size() == pathz.begin()->size()));
win_size = w;
pathz.clear();
- cout << "nway: thres = " << threshold
- << ", soft threo = " << soft_thres << endl;
+ //cout << "nway: thres = " << threshold
+ // << ", soft threo = " << soft_thres << endl;
}
void
refined_pathz.clear();
- cout << "path number is: " << pathz.size() << endl;
+ //cout << "path number is: " << pathz.size() << endl;
pathz_i = pathz.begin();
// only try to extend when pathz isn't empty.
if (extending)
{
- //cout << "Raaaawwwrrr! I am extending\n";
win_ext_len++;
}
else
}
}
}
- cout << "r_path number is: " << refined_pathz.size() << endl;
+ //cout << "r_path number is: " << refined_pathz.size() << endl;
}
threshold = atoi (data.c_str());
file_data_line = file_data_line.substr(space_split_i+1);
- cout << "seq_num=" << species_num << " win=" << win_size;
- cout << " thres=" << threshold << endl;
+ //cout << "seq_num=" << species_num << " win=" << win_size;
+ //cout << " thres=" << threshold << endl;
// clear out the current data
refined_pathz.clear();
++new_nodes_i;
}
}
-// cout << " ****I have the power...\n";
-
/* use this if I ever get the friggin seqcomp match lists to sort...
if (binary_search(trans_check_nodes.begin(), trans_check_nodes.end(),
*new_nodes_i))
pathz.clear();
window_num = all_comparisons[0][1].size();
- cout << window_num << endl;
// loop thru all windows in first species
for (win_i = 0; win_i < window_num; win_i++)
{
new_nodes_end = new_nodes.end();
//if (new_nodes_i != new_nodes_end)
//cout << "* species 0 node: " << win_i << endl;
- // cout << "foookin hell\n";
path.push_back(0);
while(new_nodes_i != new_nodes_end)
{
{
}
-Sequence::Sequence(string seq):sequence(filter_sequence(seq))
+Sequence::Sequence(string seq)
{
+ set_filtered_sequence(seq);
+}
+
+Sequence &Sequence::operator=(const Sequence& s)
+{
+ if (this != &s) {
+ sequence = s.sequence;
+ header = s.header;
+ species = s.species;
+ annots = s.annots;
+ }
+ return *this;
+}
+
+Sequence &Sequence::operator=(const std::string& s)
+{
+ set_filtered_sequence(s);
+ return *this;
}
//! load a fasta file into a sequence
data_file.open(file_path.c_str(), ios::in);
if (!data_file)
- {
+ {
throw mussa_load_error("Sequence File: " + file_path + " not found");
}
// if file opened okay, read it
end_index = sequence_raw.size();
// sequence filtering for upcasing agctn and convert non AGCTN to N
- sequence = filter_sequence(sequence_raw, start_index, end_index-start_index);
+ set_filtered_sequence(sequence_raw, start_index, end_index-start_index);
}
}
-string Sequence::filter_sequence(const string &old_seq,
- string::size_type start,
- string::size_type count) const
+void Sequence::set_filtered_sequence(const string &old_seq,
+ string::size_type start,
+ string::size_type count)
{
char conversionTable[257];
- string new_seq;
if ( count == 0)
count = old_seq.size() - start;
- new_seq.reserve(count);
+ sequence.clear();
+ sequence.reserve(count);
// Make a conversion table
// finally, the actual conversion loop
for(string::size_type seq_index = 0; seq_index < count; seq_index++)
{
- new_seq += conversionTable[ (int)old_seq[seq_index+start]];
+ sequence += conversionTable[ (int)old_seq[seq_index+start]];
}
- return new_seq;
}
// this doesn't work properly under gcc 3.x ... it can't recognize toupper
an_annot.end = atoi (annot_value.c_str());
file_data_line = file_data_line.substr(space_split_i+1);
- cout << "seq, annots: " << an_annot.start << ", " << an_annot.end
- << endl;
+ //cout << "seq, annots: " << an_annot.start << ", " << an_annot.end
+ // << endl;
// get annot name
space_split_i = file_data_line.find(" ");
}
}
+bool Sequence::empty() const
+{
+ return (size() == 0);
+}
+
const std::list<annot> Sequence::annotations() const
{
return annots;
return sequence.size();
}
+Sequence::iterator Sequence::begin()
+{
+ return sequence.begin();
+}
+
+Sequence::const_iterator Sequence::begin() const
+{
+ return sequence.begin();
+}
+
+Sequence::iterator Sequence::end()
+{
+ return sequence.end();
+}
+
+Sequence::const_iterator Sequence::end() const
+{
+ return sequence.end();
+}
+
const string&
Sequence::get_seq() const
void
Sequence::set_seq(const string& a_seq)
{
- sequence = filter_sequence(a_seq);
+ set_filtered_sequence(a_seq);
}
annot_value = file_data_line.substr(0,space_split_i);
an_annot.end = atoi (annot_value.c_str());
- if (space_split_i == string::npos) // no entry for type or name
- {
- cout << "seq, annots - no type or name\n";
- an_annot.type = "";
- an_annot.name = "";
- }
- else // else get annot type
- {
- file_data_line = file_data_line.substr(space_split_i+1);
- space_split_i = file_data_line.find(" ");
- annot_value = file_data_line.substr(0,space_split_i);
- an_annot.type = annot_value;
- if (space_split_i == string::npos) // no entry for name
- {
- cout << "seq, annots - no name\n";
- an_annot.name = "";
- }
- else // get annot name
- {
- file_data_line = file_data_line.substr(space_split_i+1);
- space_split_i = file_data_line.find(" ");
- annot_value = file_data_line.substr(0,space_split_i);
- an_annot.type = annot_value;
- }
- }
- annots.push_back(an_annot); // don't forget to actually add the annot
+ if (space_split_i == string::npos) // no entry for type or name
+ {
+ cout << "seq, annots - no type or name\n";
+ an_annot.type = "";
+ an_annot.name = "";
+ }
+ else // else get annot type
+ {
+ file_data_line = file_data_line.substr(space_split_i+1);
+ space_split_i = file_data_line.find(" ");
+ annot_value = file_data_line.substr(0,space_split_i);
+ an_annot.type = annot_value;
+ if (space_split_i == string::npos) // no entry for name
+ {
+ cout << "seq, annots - no name\n";
+ an_annot.name = "";
+ }
+ else // get annot name
+ {
+ file_data_line = file_data_line.substr(space_split_i+1);
+ space_split_i = file_data_line.find(" ");
+ annot_value = file_data_line.substr(0,space_split_i);
+ an_annot.type = annot_value;
+ }
+ }
+ annots.push_back(an_annot); // don't forget to actually add the annot
}
- cout << "seq, annots: " << an_annot.start << ", " << an_annot.end
- << "-->" << an_annot.type << "::" << an_annot.name << endl;
+ //cout << "seq, annots: " << an_annot.start << ", " << an_annot.end
+ // << "-->" << an_annot.type << "::" << an_annot.name << endl;
}
}
load_file.close();
rev_comp += conversionTable[table_i];
}
- cout << "seq: " << a_motif << endl;
- cout << "rc: " << rev_comp << endl;
+ //cout << "seq: " << a_motif << endl;
+ //cout << "rc: " << rev_comp << endl;
return rev_comp;
}
valid_motif += 'B';
}
- cout << "valid_motif is: " << valid_motif << endl;
+ //cout << "valid_motif is: " << valid_motif << endl;
return valid_motif;
}
motif_match_starts.clear();
- cout << "motif is: " << a_motif << endl;
+ //cout << "motif is: " << a_motif << endl;
a_motif = motif_validate(a_motif);
//cout << "motif is: " << a_motif << endl;
if (a_motif != "")
{
- cout << "Sequence: none blank motif\n";
+ //cout << "Sequence: none blank motif\n";
motif_scan(a_motif, &motif_match_starts);
a_motif_rc = rc_motif(a_motif);
seq_i = 0;
while (seq_i < sequence.length())
{
- cout << seq_c[seq_i];
- //if ((seq_i > 10885) && (seq_i < 10917))
+ //cout << seq_c[seq_i];
//cout << seq_c[seq_i] << "?" << a_motif[motif_i] << ":" << motif_i << " ";
// this is pretty much a straight translation of Nora's python code
// to match iupac letter codes
cout << endl;
}
-/*
- // get annot start index
- space_split_i = file_data_line.find(" ");
- annot_value = file_data_line.substr(0,space_split_i);
- an_annot.name = annot_value;
- file_data_line = file_data_line.substr(space_split_i+1);
- // get annot start index
- space_split_i = file_data_line.find(" ");
- annot_value = file_data_line.substr(0,space_split_i);
- an_annot.type = annot_value;
-*/
void motif_scan(std::string a_motif, std::vector<int> * motif_match_starts);
std::string rc_motif(std::string a_motif);
std::string motif_validate(std::string a_motif);
-
public:
+ typedef std::string::iterator iterator;
+ typedef std::string::const_iterator const_iterator;
+ typedef std::string::size_type size_type;
+
Sequence();
Sequence(std::string seq);
//! assignment to constant sequences
- //Sequence &operator=(const Sequence&);
+ Sequence &operator=(const Sequence&);
+ Sequence &operator=(const std::string &);
+
+ //! set sequence to a (sub)string containing nothing but AGCTN
+ void set_filtered_sequence(const std::string& seq,
+ std::string::size_type start=0,
+ std::string::size_type count=0);
+
//! load sequence AGCT from fasta file
- //! throws mussa_load_error if something goes wrong
+ /*! \throws mussa_load_error
+ */
void load_fasta(const std::string file_path, int seq_num=1,
int start_index=0, int end_index=0);
//! load sequence annotations
- //! throws mussa_load_error if something goes wrong
+ /*! \throws mussa_load_error
+ */
void load_annot(const std::string file_path, int start_index, int end_index);
const std::list<annot> annotations() const;
+
+ // simple access functions
void set_seq(const std::string& a_seq);
const std::string& get_seq() const;
- const char * c_seq() const;
+ const char *c_seq() const;
std::string subseq(int start, int end) const;
std::string rev_comp() const;
- //! the number of base pairs in this sequence
- std::string::size_type size() const;
+
+ // string-like functions
+ iterator begin();
+ const_iterator begin() const;
+ iterator end();
+ const_iterator end() const;
+ bool empty() const;
//! alias for size, (included as this is much like a string)
std::string::size_type length() const;
+ //! the number of base pairs in this sequence
+ std::string::size_type size() const;
+ void clear();
+
const std::string& get_header() const;
- //std::string sp_name() const;
std::vector<int> find_motif(std::string a_motif);
- void clear();
void save(std::fstream &save_file);
void load_museq(std::string load_file_path, int seq_num);
- // FIXME: filter sequence should be a class not instance function
- //! convert a (sub)string to being just AGCTN mapping
- //! mapping lower case to upper case, and any invalid char to N
- std::string filter_sequence(const std::string& seq,
- std::string::size_type start=0,
- std::string::size_type count=0) const;
};
#endif
SOURCES.cxx := test_conserved_path.cxx \
test_flp.cxx \
test_glsequence.cxx \
+ test_color.cxx \
test_main.cxx \
test_mussa.cxx \
test_nway.cxx \
--- /dev/null
+#include <boost/test/auto_unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+#include <string>
+
+#include "alg/color.h"
+
+using namespace std;
+
+BOOST_AUTO_TEST_CASE ( color )
+{
+ const float tolerance = 1e-6;
+ const float zero = 0.0;
+ const float dotone = 0.1;
+ const float dottwo = 0.2;
+ const float dotthree = 0.3;
+ const float dotfour = 0.4;
+ Color c1;
+ BOOST_CHECK_CLOSE (c1.r(), zero, tolerance );
+ BOOST_CHECK_CLOSE (c1.g(), zero, tolerance );
+ BOOST_CHECK_CLOSE (c1.b(), zero, tolerance );
+ BOOST_CHECK_CLOSE (c1.a(), zero, tolerance );
+
+ Color c2(dotone, dottwo, dotthree, dotfour);
+ BOOST_CHECK_CLOSE (c2.r(), dotone, tolerance );
+ BOOST_CHECK_CLOSE (c2.g(), dottwo, tolerance );
+ BOOST_CHECK_CLOSE (c2.b(), dotthree, tolerance );
+ BOOST_CHECK_CLOSE (c2.a(), dotfour, tolerance );
+
+ Color c3(c2);
+ BOOST_CHECK_CLOSE (c3.r(), dotone, tolerance );
+ BOOST_CHECK_CLOSE (c3.g(), dottwo, tolerance );
+ BOOST_CHECK_CLOSE (c3.b(), dotthree, tolerance );
+ BOOST_CHECK_CLOSE (c3.a(), dotfour, tolerance );
+ BOOST_CHECK_EQUAL ( c2, c3 );
+
+ const float *colors = c3.get();
+ BOOST_CHECK_CLOSE (colors[Color::RedChannel ], dotone, tolerance );
+ BOOST_CHECK_CLOSE (colors[Color::GreenChannel ], dottwo, tolerance );
+ BOOST_CHECK_CLOSE (colors[Color::BlueChannel ], dotthree, tolerance );
+ BOOST_CHECK_CLOSE (colors[Color::AlphaChannel ], dotfour, tolerance );
+}
+
+
GlSequence glseq0(seq0);
BOOST_CHECK (glseq0.sequence().get_seq() == s0);
- // width of a sequence should be number of base pairs (aka chars)
- BOOST_CHECK (glseq0.width() == s0.size());
GlSequence glseq1(seq1);
GlSequence glseq_copy0(glseq0);
BOOST_CHECK( glseq0.sequence().get_seq() == s1 );
}
+
+BOOST_AUTO_TEST_CASE( glsequence_color )
+{
+ Color black(0.0, 0.0, 0.0, 0.0);
+ Color c(0.1, 0.2, 0.3, 0.4);
+ Sequence seq("AAGGCCTT");
+ GlSequence s(seq);
+
+ BOOST_CHECK_EQUAL(s.color(), black );
+ s.setColor( c );
+ BOOST_CHECK_EQUAL( s.color(), c );
+}
+
+BOOST_AUTO_TEST_CASE( glsequence_renderable )
+{
+ Sequence seq("AAGGCCTT");
+ GlSequence s(seq);
+
+ // way more base pairs than viewport pixel width
+ BOOST_CHECK_EQUAL(s.is_sequence_renderable( 0, 1000, 500), false );
+ // way fewer basepairs than viewport pixel width
+ BOOST_CHECK_EQUAL(s.is_sequence_renderable( 0, 10, 500), true);
+}
+
+BOOST_AUTO_TEST_CASE( glsequence_sequence )
+{
+ string seq_string("AAGGCCTTNNAAGGCCTTNNAAGGCCTTNN");
+ string::size_type seqlen = seq_string.size();
+ Sequence seq(seq_string);
+ GlSequence glseq(seq);
+
+ BOOST_CHECK( glseq.sequence_begin(0, 50) == seq.begin() );
+ // always make sure we return seq.end() regardless of how much extra
+ // is asked for
+ BOOST_CHECK( glseq.sequence_end(0, seqlen+10) == seq.end() );
+ // do we get the right end pointer?
+ BOOST_CHECK( glseq.sequence_end(0, 5) == seq.begin()+5 );
+
+ // when we request far too much sequence what do we get?
+ BOOST_CHECK( glseq.sequence_begin(seqlen+10, seqlen+20) == seq.end() );
+ BOOST_CHECK( glseq.sequence_end(seqlen+10, seqlen+20) == seq.end() );
+
+ // we cant ask for reversed sequences with sequence_begin/end
+ BOOST_CHECK( glseq.sequence_begin(10, 5) == seq.end() );
+ BOOST_CHECK( glseq.sequence_end(10, 5) == seq.end() );
+
+ Sequence::const_iterator seq_itor;
+
+ // if we as for an empty segment? start and end should equal
+ seq_itor = glseq.sequence_begin(10, 10);
+ BOOST_CHECK( seq_itor == glseq.sequence_end(10, 10) );
+
+ // reuse seq_itor downhere
+ string::const_iterator str_itor;
+ for(str_itor = seq.begin(),
+ seq_itor = glseq.sequence_begin();
+ str_itor != seq.end() and
+ seq_itor != glseq.sequence_end();
+ ++str_itor, ++seq_itor)
+ {
+ BOOST_CHECK_EQUAL( *str_itor, *seq_itor );
+ }
+}
+
+// make sure the computation of the leftmost and rightmost base is correct
+BOOST_AUTO_TEST_CASE( glsequence_leftright_base )
+{
+ std::string seq_string = "AAGGCCTT";
+ Sequence seq(seq_string);
+ GlSequence glseq(seq);
+
+ BOOST_CHECK_EQUAL( glseq.leftbase( -50.0 ), 0 );
+ BOOST_CHECK_EQUAL( glseq.leftbase( 0.5 ), 1 );
+ BOOST_CHECK_EQUAL( glseq.rightbase( 1000.0 ), seq_string.size() );
+ BOOST_CHECK_EQUAL( glseq.rightbase( seq_string.size()-0.5),
+ seq_string.size()-1);
+}
string s0("ATATGCGC");
string s1("GGGGGGGC");
Sequence seq1(s1);
- cout << "seq1 rev = " << seq1.rev_comp() << endl;
Mussa analysis;
analysis.add_a_seq(s0);
BOOST_CHECK_EQUAL(s3.get_seq(), "ANNNG");
BOOST_CHECK_EQUAL(s3.subseq(0,2), "AN");
- BOOST_CHECK_EQUAL(s3.filter_sequence("AAGGCCTT", 0, 2), "AA");
- BOOST_CHECK_EQUAL(s3.filter_sequence("AAGGCCTT", 2, 2), "GG");
- BOOST_CHECK_EQUAL(s3.filter_sequence("AAGGCCTT", 4), "CCTT");
+ s3.set_filtered_sequence("AAGGCCTT", 0, 2);
+ BOOST_CHECK_EQUAL(s3.get_seq(), "AA");
+ s3.set_filtered_sequence("AAGGCCTT", 2, 2);
+ BOOST_CHECK_EQUAL( s3.get_seq(), "GG");
+ s3.set_filtered_sequence("AAGGCCTT", 4);
+ BOOST_CHECK_EQUAL( s3.get_seq(), "CCTT");
s3.clear();
BOOST_CHECK_EQUAL(s3.get_seq(), "");
"muscle creatine kinase gene (CKMM), "
"5' flank");
}
+
+BOOST_AUTO_TEST_CASE ( sequence_empty )
+{
+ Sequence s;
+ BOOST_CHECK_EQUAL( s.empty(), true );
+ s = "AAAGGG";
+ BOOST_CHECK_EQUAL( s.empty(), false );
+}
+
+BOOST_AUTO_TEST_CASE ( sequence_iterators )
+{
+ std::string seq_string = "AAGGCCTTNNTATA";
+ Sequence s(seq_string);
+ const Sequence cs(s);
+ std::string::size_type count = 0;
+
+ std::string::iterator str_itor;
+ Sequence::iterator s_itor;
+ Sequence::const_iterator cs_itor;
+
+ for( str_itor = seq_string.begin(),
+ s_itor = s.begin(),
+ cs_itor = cs.begin();
+ str_itor != seq_string.end() and
+ s_itor != s.end() and
+ cs_itor != cs.end();
+ ++str_itor, ++s_itor, ++cs_itor, ++count)
+ {
+ BOOST_CHECK_EQUAL ( *str_itor, *s_itor );
+ BOOST_CHECK_EQUAL ( *s_itor, *cs_itor );
+ BOOST_CHECK_EQUAL ( *cs_itor, *str_itor );
+ }
+ BOOST_CHECK_EQUAL( seq_string.size(), count );
+ BOOST_CHECK_EQUAL( s.size(), count );
+ BOOST_CHECK_EQUAL( cs.size(), count );
+}
qui/ThresholdWidget.h \
qui/ImageScaler.h \
qui/ImageSaveDialog.h \
+ alg/color.h \
alg/conserved_path.h \
alg/flp.hh \
alg/glsequence.h \
qui/ThresholdWidget.cxx \
qui/ImageScaler.cxx \
qui/ImageSaveDialog.cxx \
+ alg/color.cxx \
alg/conserved_path.cxx \
alg/flp.cxx \
alg/flp_seqcomp.cxx \
return_internal_reference<>())
.add_property("x", &GlSequence::x, &GlSequence::setX)
.add_property("y", &GlSequence::y, &GlSequence::setY)
- .add_property("width", &GlSequence::width, &GlSequence::setWidth)
+ .add_property("length", &GlSequence::length)
;
}
return viewport_center;
}
-static float max(float a, float b)
-{
- if ( a < b)
- return b;
- else
- return a;
-}
-
-static float min(float a, float b)
-{
- if ( a < b)
- return a;
- else
- return b;
-}
-
void PathScene::updateViewport(float center, int new_zoom)
{
float max_width = maxOrtho2d.width();
gl_seq = new GlSequence(*seq_itor);
gl_seq->setX(0);
gl_seq->setY(y);
- if (gl_seq->width() > max_base_pairs ) max_base_pairs = gl_seq->width();
+ if (gl_seq->length() > max_base_pairs ) max_base_pairs = gl_seq->length();
tracks.push_back( *gl_seq );
}
maxOrtho2d.setWidth(max_base_pairs + 100.0);
glColor3f(0.0, 0.0, 1.0);
} else {
// else be dim
- glColor3f(0.5, 0.5, 1.0);
+ glColor3f(0.7, 0.7, 1.0);
}
} else {
// both current and previous path have the same orientation
if (not selectedMode or selectedPaths[objid] == true) {
- glColor3f(1.0, 0.0, 0.2);
+ glColor3f(1.0, 0.0, 0.0);
} else {
- glColor3f(1.0, 0.5, 0.5);
+ glColor3f(1.0, 0.7, 0.7);
}
}
prevReversed = reversed;
projects / mussa.git / commitdiff