// ---------- 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),
- alphabet(alphabet_),
- seq_start(0),
- seq_count(0),
+Sequence::Sequence(AlphabetRef alphabet)
+ : seq(new SeqSpan("", alphabet)),
strand(UnknownStrand)
{
}
{
}
-Sequence::Sequence(const char *seq, alphabet_ref alphabet_)
- : parent(0),
- alphabet(alphabet_),
- seq_start(0),
- seq_count(0),
- strand(UnknownStrand),
+Sequence::Sequence(const char *seq, AlphabetRef alphabet_)
+ : strand(UnknownStrand),
header(""),
species("")
{
- set_filtered_sequence(seq, alphabet);
+ set_filtered_sequence(seq, alphabet_);
}
-Sequence::Sequence(const std::string& seq, alphabet_ref alphabet_)
- : parent(0),
- alphabet(alphabet_),
- seq_start(0),
- seq_count(0),
- strand(UnknownStrand),
+Sequence::Sequence(const std::string& seq, AlphabetRef alphabet_)
+ : strand(UnknownStrand),
header(""),
species("")
{
- set_filtered_sequence(seq, alphabet);
+ set_filtered_sequence(seq, alphabet_);
}
Sequence::Sequence(const Sequence& o)
- : parent(o.parent),
- seq(o.seq),
- alphabet(o.alphabet),
- seq_start(o.seq_start),
- seq_count(o.seq_count),
+ : seq(o.seq),
strand(o.strand),
header(o.header),
species(o.species),
{
}
+Sequence::Sequence(const Sequence* o)
+ : seq(o->seq),
+ strand(o->strand),
+ header(o->header),
+ species(o->species),
+ annots(o->annots),
+ motif_list(o->motif_list)
+{
+}
+
+Sequence::Sequence(const SeqSpanRef& seq_ref)
+ : seq(seq_ref),
+ 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;
void Sequence::load_fasta(fs::path file_path, int seq_num, int start_index, int end_index)
{
- load_fasta(file_path, alphabet, seq_num, start_index, end_index);
+ load_fasta(file_path, reduced_nucleic_alphabet, seq_num, start_index, end_index);
}
//! load a fasta file into a sequence
-void Sequence::load_fasta(fs::path file_path, alphabet_ref a,
+void Sequence::load_fasta(fs::path file_path, AlphabetRef a,
int seq_num, int start_index, int end_index)
{
fs::fstream data_file;
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());
}
}
}
void Sequence::load_fasta(std::istream& file,
int seq_num, int start_index, int end_index)
{
- load_fasta(file, alphabet, seq_num, start_index, end_index);
+ load_fasta(file, reduced_nucleic_alphabet, seq_num, start_index, end_index);
}
void
-Sequence::load_fasta(std::istream& data_file, alphabet_ref a,
+Sequence::load_fasta(std::istream& data_file, AlphabetRef a,
int seq_num,
int start_index, int end_index)
{
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;
std::string seq_tmp; // holds sequence during basic filtering
- const Alphabet &alpha = get_alphabet(a);
+ const Alphabet &alpha = get_alphabet();
if (seq_num == 0) {
throw mussa_load_error("fasta sequence number is 1 based (can't be 0)");
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());
}
}
}
}
void Sequence::set_filtered_sequence(const std::string &in_seq,
- alphabet_ref alphabet_,
+ AlphabetRef alphabet_,
size_type start,
size_type count,
strand_type strand_)
{
- 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
+ const Alphabet& alpha_impl = Alphabet::get_alphabet(alphabet_); // go get one of our actual alphabets
std::string::const_iterator seq_i = in_seq.begin()+start;
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, alphabet_));
+ 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::create_reverse_map() const
-{
- std::string rc_map(256, '~');
- // if we're rna, use U instead of T
- // we might want to add an "is_rna" to sequence at somepoint
- char TU = (alphabet == reduced_rna_alphabet) ? 'U' : 'T';
- char tu = (alphabet == reduced_rna_alphabet) ? 'u' : 't';
- rc_map['A'] = TU ; rc_map['a'] = tu ;
- rc_map['T'] = 'A'; rc_map['t'] = 'a';
- rc_map['U'] = 'A'; rc_map['u'] = 'a';
- rc_map['G'] = 'C'; rc_map['g'] = 'c';
- rc_map['C'] = 'G'; rc_map['c'] = 'g';
- rc_map['M'] = 'K'; rc_map['m'] = 'k';
- rc_map['R'] = 'Y'; rc_map['r'] = 'y';
- rc_map['W'] = 'W'; rc_map['w'] = 'w';
- rc_map['S'] = 'S'; rc_map['s'] = 's';
- rc_map['Y'] = 'R'; rc_map['y'] = 'r';
- rc_map['K'] = 'M'; rc_map['k'] = 'm';
- rc_map['V'] = 'B'; rc_map['v'] = 'b';
- rc_map['H'] = 'D'; rc_map['h'] = 'd';
- rc_map['D'] = 'H'; rc_map['d'] = 'h';
- rc_map['B'] = 'V'; rc_map['b'] = 'v';
- rc_map['N'] = 'N'; rc_map['n'] = 'n';
- rc_map['X'] = 'X'; rc_map['x'] = 'x';
- rc_map['?'] = '?';
- rc_map['.'] = '.';
- rc_map['-'] = '-';
- rc_map['~'] = '~'; // not really needed, but perhaps it's clearer.
- return rc_map;
-}
+// FIXME: This needs to be moved into SeqSpan
Sequence Sequence::rev_comp() const
{
std::string rev_comp;
rev_comp.reserve(length());
- std::string rc_map = create_reverse_map();
+ std::string rc_map = seq->seq->create_reverse_complement_map();
// reverse and convert
Sequence::const_reverse_iterator seq_i;
{
rev_comp.append(1, rc_map[*seq_i]);
}
- return Sequence(rev_comp, alphabet);
+ return Sequence(rev_comp, seq->seq->get_alphabet_ref());
+}
+
+const Alphabet& Sequence::get_alphabet() const
+{
+ return (seq) ? seq->get_alphabet() : Alphabet::empty_alphabet();
}
void Sequence::set_fasta_header(std::string header_)
return "";
}
-const Alphabet& Sequence::get_alphabet() const
-{
- return get_alphabet(alphabet);
-}
-
-const Alphabet& Sequence::get_alphabet(alphabet_ref alpha) const
-{
- switch (alpha) {
- case reduced_dna_alphabet:
- return Alphabet::reduced_dna_alphabet();
- case reduced_rna_alphabet:
- return Alphabet::reduced_rna_alphabet();
- case reduced_nucleic_alphabet:
- return Alphabet::reduced_nucleic_alphabet();
- case nucleic_alphabet:
- return Alphabet::nucleic_alphabet();
- case protein_alphabet:
- return Alphabet::protein_alphabet();
- default:
- throw std::runtime_error("unrecognized alphabet type");
- break;
- }
-}
-
-void Sequence::set_sequence(const std::string& s, alphabet_ref a)
+void Sequence::set_sequence(const std::string& s, AlphabetRef a)
{
set_filtered_sequence(s, a);
}
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)
{
// this is pretty much a straight translation of Nora's python code
// to match iupac letter codes
motif_char = toupper(a_motif[motif_i]);
- seq_char = toupper(seq->at(seq_start+seq_i));
+ seq_char = toupper(seq->at(seq_i));
if (motif_char =='N')
motif_i++;
else if (motif_char == seq_char)
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