// ---------- sequence.cc -----------
// ----------------------------------------
#include <boost/filesystem/fstream.hpp>
+#include <boost/filesystem/operations.hpp>
namespace fs = boost::filesystem;
#include <boost/spirit/core.hpp>
Sequence::Sequence(alphabet_ref alphabet_)
- : parent(0),
+ : seq(new SeqSpan("")),
alphabet(alphabet_),
- seq_start(0),
- seq_count(0),
strand(UnknownStrand)
{
}
}
Sequence::Sequence(const char *seq, alphabet_ref alphabet_)
- : parent(0),
- alphabet(alphabet_),
- seq_start(0),
- seq_count(0),
+ : alphabet(alphabet_),
strand(UnknownStrand),
header(""),
species("")
}
Sequence::Sequence(const std::string& seq, alphabet_ref alphabet_)
- : parent(0),
- alphabet(alphabet_),
- seq_start(0),
- seq_count(0),
+ : alphabet(alphabet_),
strand(UnknownStrand),
header(""),
species("")
}
Sequence::Sequence(const Sequence& o)
- : parent(o.parent),
- seq(o.seq),
+ : seq(o.seq),
alphabet(o.alphabet),
- seq_start(o.seq_start),
- seq_count(o.seq_count),
strand(o.strand),
header(o.header),
species(o.species),
{
}
+Sequence::Sequence(const Sequence* o)
+ : seq(o->seq),
+ alphabet(o->alphabet),
+ strand(o->strand),
+ header(o->header),
+ species(o->species),
+ annots(o->annots),
+ motif_list(o->motif_list)
+{
+}
+
+Sequence::Sequence(const SeqSpanRef& seq_ref, alphabet_ref alphabet_)
+ : seq(seq_ref),
+ alphabet(alphabet),
+ strand(UnknownStrand),
+ header(""),
+ species("")
+{
+}
+
Sequence &Sequence::operator=(const Sequence& s)
{
if (this != &s) {
- parent = s.parent;
seq = s.seq;
alphabet = s.alphabet;
- seq_start = s.seq_start;
- seq_count = s.seq_count;
strand = s.strand;
header = s.header;
species = s.species;
errormsg << file_path.native_file_string()
<< " did not have any fasta sequences" << std::endl;
throw sequence_empty_file_error(errormsg.str());
+ } catch(sequence_invalid_load_error e) {
+ std::ostringstream msg;
+ msg << file_path.native_file_string();
+ msg << " " << e.what();
+ throw sequence_invalid_load_error(msg.str());
}
}
}
{
std::string file_data_line;
int header_counter = 0;
+ size_t line_counter = 0;
bool read_seq = true;
std::string rev_comp;
std::string sequence_raw;
while ( (!data_file.eof()) && (header_counter < seq_num) )
{
multiplatform_getline(data_file, file_data_line);
+ ++line_counter;
if (file_data_line.substr(0,1) == ">")
header_counter++;
}
while ( !data_file.eof() && read_seq ) {
multiplatform_getline(data_file,file_data_line);
+ ++line_counter;
if (file_data_line.substr(0,1) == ">")
read_seq = false;
else {
if(alpha.exists(*line_i)) {
sequence_raw += *line_i;
} else {
- throw sequence_invalid_load_error("Unrecognized characters in fasta sequence");
+ std::ostringstream msg;
+ msg << "Unrecognized characters in fasta sequence at line ";
+ msg << line_counter;
+ throw sequence_invalid_load_error(msg.str());
}
}
}
alphabet = alphabet_;
if ( count == npos)
count = in_seq.size() - start;
- boost::shared_ptr<seq_string> new_seq(new seq_string);
- new_seq->reserve(count);
+ std::string new_seq;
+ new_seq.reserve(count);
// finally, the actual conversion loop
const Alphabet& alpha_impl = get_alphabet(); // go get one of our actual alphabets
for(size_type i = 0; i != count; ++i, ++seq_i)
{
if (alpha_impl.exists(*seq_i)) {
- new_seq->append(1, toupper(*seq_i));
+ new_seq.append(1, toupper(*seq_i));
} else {
- new_seq->append(1, 'N');
+ new_seq.append(1, 'N');
}
}
- parent = 0;
- seq = new_seq;
- seq_start = 0;
- seq_count = count;
+ SeqSpanRef new_seq_ref(new SeqSpan(new_seq));
+ seq = new_seq_ref;
strand = strand_;
}
void
Sequence::load_annot(fs::path file_path, int start_index, int end_index)
{
+ if (not fs::exists(file_path)) {
+ throw mussa_load_error("Annotation File " + file_path.string() + " was not found");
+ }
+ if (fs::is_directory(file_path)) {
+ throw mussa_load_error(file_path.string() +
+ " is a directory, please provide a file for annotations."
+ );
+ }
fs::fstream data_stream(file_path, std::ios::in);
if (!data_stream)
{
- throw mussa_load_error("Sequence File: " + file_path.string() + " not found");
+ throw mussa_load_error("Error loading annotation file " + file_path.string());
}
// so i should probably be passing the parse function some iterators
data.push_back(c);
}
data_stream.close();
-
- parse_annot(data, start_index, end_index);
+
+ try {
+ parse_annot(data, start_index, end_index);
+ } catch(annotation_load_error e) {
+ std::ostringstream msg;
+ msg << file_path.native_file_string()
+ << " "
+ << e.what();
+ throw annotation_load_error(msg.str());
+ }
}
/* If this works, yikes, this is some brain hurting code.
std::string seq;
std::list<annot> parsed_annots;
std::list<Sequence> query_seqs;
- int parsed=1;
+ int parsed=0;
bool ok = spirit::parse(data.begin(), data.end(),
(
return annots;
}
-Sequence
-Sequence::subseq(int start, int count)
+void Sequence::copy_children(Sequence &new_seq, size_type start, size_type count) const
{
- if (!seq) {
- Sequence new_seq;
- return new_seq;
- }
-
- // there might be an off by one error with start+count > size()
- if ( count == npos || start+count > size()) {
- count = size()-start;
- }
- Sequence new_seq(*this);
- new_seq.parent = this;
- new_seq.seq_start = seq_start+start;
- new_seq.seq_count = count;
-
new_seq.motif_list = motif_list;
new_seq.annots.clear();
- // attempt to copy & reannotate position based annotations
- int end = start+count;
for(std::list<annot>::const_iterator annot_i = annots.begin();
annot_i != annots.end();
++annot_i)
{
- int annot_begin= annot_i->begin;
- int annot_end = annot_i->end;
+ size_type annot_begin= annot_i->begin;
+ size_type annot_end = annot_i->end;
- if (annot_begin < end) {
+ if (annot_begin < start+count) {
if (annot_begin >= start) {
annot_begin -= start;
} else {
annot_begin = 0;
}
- if (annot_end < end) {
+ if (annot_end < start+count) {
annot_end -= start;
} else {
annot_end = count;
new_seq.annots.push_back(new_annot);
}
}
+
+}
+Sequence
+Sequence::subseq(size_type start, size_type count) const
+{
+ if (!seq) {
+ Sequence new_seq;
+ return new_seq;
+ }
+ Sequence new_seq = *this;
+ new_seq.seq = seq->subseq(start, count);
+
+ copy_children(new_seq, start, count);
+
return new_seq;
}
std::string Sequence::get_sequence() const
{
- if (seq)
- return *seq;
- else
- return std::string();
+ return seq->sequence();
}
Sequence::const_reference Sequence::operator[](Sequence::size_type i) const
return at(i);
}
-Sequence::const_reference Sequence::at(Sequence::size_type i) const
-{
- if (!seq) throw std::out_of_range("empty sequence");
- return seq->at(i+seq_start);
-}
-
void
Sequence::clear()
{
- parent = 0;
seq.reset();
- seq_start = 0;
- seq_count = 0;
strand = UnknownStrand;
header.clear();
species.clear();
motif_list.clear();
}
-const char *Sequence::c_str() const
-{
- if (seq)
- return seq->c_str()+seq_start;
- else
- return 0;
-}
-
-Sequence::const_iterator Sequence::begin() const
-{
- if (seq and seq_count != 0)
- return seq->begin()+seq_start;
- else
- return Sequence::const_iterator(0);
-}
-
-Sequence::const_iterator Sequence::end() const
-{
- if (seq and seq_count != 0) {
- return seq->begin() + seq_start + seq_count;
- } else {
- return Sequence::const_iterator(0);
- }
-}
-
-Sequence::const_reverse_iterator Sequence::rbegin() const
-{
- if (seq and seq_count != 0)
- return seq->rbegin()+(seq->size()-(seq_start+seq_count));
- else
- return Sequence::const_reverse_iterator();
-}
-
-Sequence::const_reverse_iterator Sequence::rend() const
-{
- if (seq and seq_count != 0) {
- return rbegin() + seq_count;
- } else {
- return Sequence::const_reverse_iterator();
- }
-}
-
-bool Sequence::empty() const
-{
- return (seq_count == 0) ? true : false;
-}
-
-Sequence::size_type Sequence::start() const
-{
- if (parent)
- return seq_start - parent->start();
- else
- return seq_start;
-}
-
-Sequence::size_type Sequence::stop() const
-{
- return start() + seq_count;
-}
-
-Sequence::size_type Sequence::size() const
-{
- return seq_count;
-}
-
-Sequence::size_type Sequence::length() const
-{
- return size();
-}
-
void
Sequence::save(fs::fstream &save_file)
{
Sequence::value_type motif_char;
Sequence::value_type seq_char;
- while (seq_i < seq->length())
+ while (seq_i < size())
{
// this is pretty much a straight translation of Nora's python code
// to match iupac letter codes
std::ostream& operator<<(std::ostream& out, const Sequence& s)
{
- for(Sequence::const_iterator s_i = s.begin(); s_i != s.end(); ++s_i) {
- out << *s_i;
+ if (s.seq) {
+ for(Sequence::const_iterator s_i = s.begin(); s_i != s.end(); ++s_i) {
+ out << *s_i;
+ }
}
return out;
}
}
}
-bool operator==(const Sequence& x, const Sequence& y)
+template <typename Iter1, typename Iter2>
+static
+bool sequence_insensitive_equality(Iter1 abegin, Iter1 aend, Iter2 bbegin, Iter2 bend)
{
- if (x.empty() and y.empty()) {
- // if there's no sequence in either sequence structure, they're equal
+ Iter1 aseq_i = abegin;
+ Iter2 bseq_i = bbegin;
+ if (aend-abegin == bend-bbegin) {
+ // since the length of the two sequences is equal, we only need to
+ // test one.
+ for(; aseq_i != aend; ++aseq_i, ++bseq_i) {
+ if (toupper(*aseq_i) != toupper(*bseq_i)) {
+ return false;
+ }
+ }
return true;
- } else if (x.empty() or y.empty()) {
- // if we fail the first test, and we discover one is empty,
- // we know they can't be equal. (and we need to do this
- // to prevent dereferencing an empty pointer)
- return false;
- } else if (x.seq_count != y.seq_count) {
- // if they're of different lenghts, they're not equal
+ } else {
return false;
}
- Sequence::const_iterator xseq_i = x.begin();
- Sequence::const_iterator yseq_i = y.begin();
- // since the length of the two sequences is equal, we only need to
- // test one.
- for(; xseq_i != x.end(); ++xseq_i, ++yseq_i) {
- if (toupper(*xseq_i) != toupper(*yseq_i)) {
- return false;
+}
+
+bool operator==(const Sequence& x, const Sequence& y)
+{
+ if (x.seq and y.seq) {
+ // both x and y are defined
+ if (SeqSpan::isFamily(x.seq, y.seq)) {
+ // both are part of the same SeqSpan tree
+ return *(x.seq) == *(y.seq);
+ } else {
+ // we'll have to do a real comparison
+ return sequence_insensitive_equality<SeqSpan::const_iterator, SeqSpan::const_iterator>(
+ x.begin(), x.end(),
+ y.begin(), y.end()
+ );
}
+ } else {
+ // true if they're both empty (with either a null SeqSpanRef or
+ // a zero length string
+ return (x.size() == y.size());
}
- return true;
}
bool operator!=(const Sequence& x, const Sequence& y)
{
return not operator==(x, y);
}
+
projects / mussa.git / blobdiff