#include "mussa_exceptions.hpp"
#include <string>
+#include <stdexcept>
#include <iostream>
#include <sstream>
+#include <set>
annot::annot()
: begin(0),
{
}
-const std::string Sequence::dna_alphabet("AaCcGgTtNn\012\015");
-const std::string Sequence::rna_alphabet("AaCcGgNnUu\012\015");
- //! this is the general iupac alphabet for nucleotides
-const std::string Sequence::nucleic_iupac_alphabet("AaCcGgTtUuRrYyMmKkSsWwBbDdHhVvNn\012\015");
- //! the protein alphabet
-const std::string Sequence::protein_alphabet("AaCcDdEeFfGgHhIiKkLlMmNnPpQqRrSsTtVvWwYy\012\015");
-Sequence::Sequence()
+Sequence::Sequence(alphabet_ref alphabet_)
: parent(0),
+ alphabet(alphabet_),
seq_start(0),
seq_count(0),
strand(UnknownStrand)
{
}
-Sequence::Sequence(const char *seq)
+Sequence::Sequence(const char *seq, alphabet_ref alphabet_)
: parent(0),
+ alphabet(alphabet_),
seq_start(0),
seq_count(0),
strand(UnknownStrand),
header(""),
species("")
{
- set_filtered_sequence(seq);
+ set_filtered_sequence(seq, alphabet);
}
-Sequence::Sequence(const std::string& seq)
+Sequence::Sequence(const std::string& seq, alphabet_ref alphabet_)
: parent(0),
+ alphabet(alphabet_),
seq_start(0),
seq_count(0),
strand(UnknownStrand),
header(""),
species("")
{
- set_filtered_sequence(seq);
+ 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),
strand(o.strand),
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;
}
}
+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 a fasta file into a sequence
-/*!
- * \param file_path the location of the fasta file in the filesystem
- * \param seq_num which sequence in the file to load
- * \param start_index starting position in the fasta sequence, 0 for beginning
- * \param end_index ending position in the fasta sequence, 0 for end
- * \return error message, empty string if no error. (gag!)
- */
-void Sequence::load_fasta(fs::path file_path, int seq_num,
- int start_index, int end_index)
+void Sequence::load_fasta(fs::path file_path, alphabet_ref a,
+ int seq_num, int start_index, int end_index)
{
fs::fstream data_file;
data_file.open(file_path, std::ios::in);
throw mussa_load_error("Sequence File: "+file_path.string()+" not found");
} else {
try {
- load_fasta(data_file, seq_num, start_index, end_index);
+ load_fasta(data_file, a, seq_num, start_index, end_index);
} catch(sequence_empty_error e) {
// there doesn't appear to be any sequence
// catch and rethrow to include the filename
}
}
+void Sequence::load_fasta(std::iostream& file,
+ int seq_num, int start_index, int end_index)
+{
+ load_fasta(file, alphabet, seq_num, start_index, end_index);
+}
+
void
-Sequence::load_fasta(std::iostream& data_file, int seq_num,
+Sequence::load_fasta(std::iostream& data_file, alphabet_ref a,
+ int seq_num,
int start_index, int end_index)
{
std::string file_data_line;
std::string rev_comp;
std::string sequence_raw;
std::string seq_tmp; // holds sequence during basic filtering
+ const Alphabet &alpha = get_alphabet(a);
if (seq_num == 0) {
throw mussa_load_error("fasta sequence number is 1 based (can't be 0)");
multiplatform_getline(data_file,file_data_line);
if (file_data_line.substr(0,1) == ">")
read_seq = false;
- else sequence_raw += file_data_line;
+ else {
+ for (std::string::const_iterator line_i = file_data_line.begin();
+ line_i != file_data_line.end();
+ ++line_i)
+ {
+ if(alpha.exists(*line_i)) {
+ sequence_raw += *line_i;
+ } else {
+ throw sequence_invalid_load_error("Unrecognized characters in fasta sequence");
+ }
+ }
+ }
}
// Lastly, if subselection of the sequence was specified we keep cut out
std::string msg("The selected sequence appears to be empty");
throw sequence_empty_error(msg);
}
- set_filtered_sequence(sequence_raw, start_index, end_index-start_index);
+ set_filtered_sequence(sequence_raw, a, start_index, end_index-start_index);
} else {
std::string errormsg("There were no fasta sequences");
throw sequence_empty_file_error(errormsg);
}
}
-void Sequence::set_filtered_sequence(const std::string &old_seq,
+void Sequence::set_filtered_sequence(const std::string &in_seq,
+ alphabet_ref alphabet_,
size_type start,
size_type count,
strand_type strand_)
{
- char conversionTable[257];
-
+ alphabet = alphabet_;
if ( count == npos)
- count = old_seq.size() - start;
+ count = in_seq.size() - start;
boost::shared_ptr<seq_string> new_seq(new seq_string);
new_seq->reserve(count);
- // Make a conversion table
-
- // everything we don't specify below will become 'N'
- for(int table_i=0; table_i < 256; table_i++)
- {
- conversionTable[table_i] = 'N';
- }
- // add end of string character for printing out table for testing purposes
- conversionTable[256] = '\0';
-
- // we want these to map to themselves - ie not to change
- conversionTable[(int)'A'] = 'A';
- conversionTable[(int)'T'] = 'T';
- conversionTable[(int)'G'] = 'G';
- conversionTable[(int)'C'] = 'C';
- // this is to upcase
- conversionTable[(int)'a'] = 'A';
- conversionTable[(int)'t'] = 'T';
- conversionTable[(int)'g'] = 'G';
- conversionTable[(int)'c'] = 'C';
-
// finally, the actual conversion loop
- for(std::string::size_type seq_index = 0; seq_index < count; seq_index++)
+ const Alphabet& alpha_impl = get_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)
{
- new_seq->append(1, conversionTable[ (int)old_seq[seq_index+start]]);
+ if (alpha_impl.exists(*seq_i)) {
+ new_seq->append(1, *seq_i);
+ } else {
+ new_seq->append(1, 'N');
+ }
}
parent = 0;
seq = new_seq;
((spirit::ch_p('>')|spirit::str_p(">")) >>
(*(spirit::print_p))[spirit::assign_a(name)] >>
spirit::eol_p >>
- (+(spirit::chset<>(nucleic_iupac_alphabet.c_str())))[spirit::assign_a(seq)]
+ (+(spirit::chset<>(Alphabet::nucleic_alphabet.c_str())))[spirit::assign_a(seq)]
)[push_back_seq(query_seqs, name, seq)]
) >>
*spirit::space_p
return new_seq;
}
-std::string
-Sequence::rev_comp() const
+std::string Sequence::create_reverse_map() const
{
- std::string rev_comp;
- char conversionTable[257];
- int seq_i, table_i, len;
-
- len = length();
- rev_comp.reserve(len);
- // make a conversion table
- // init all parts of conversion table to '~' character
- // '~' I doubt will ever appear in a sequence file (jeez, I hope)
- // and may the fleas of 1000 camels infest the genitals of any biologist (and
- // seven generations of their progeny) who decides to make it mean
- // something special!!!
- // PS - double the curse for any smartass non-biologist who tries it as well
- for(table_i=0; table_i < 256; table_i++)
- {
- conversionTable[table_i] = '~';
- }
- // add end of string character for printing out table for testing purposes
- conversionTable[256] = '\0';
-
- // add in the characters for the bases we want to convert
- conversionTable[(int)'A'] = 'T';
- conversionTable[(int)'T'] = 'A';
- conversionTable[(int)'G'] = 'C';
- conversionTable[(int)'C'] = 'G';
- conversionTable[(int)'N'] = 'N';
+ 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;
+}
- // finally, the actual conversion loop
- for(seq_i = len - 1; seq_i >= 0; seq_i--)
+Sequence Sequence::rev_comp() const
+{
+ std::string rev_comp;
+ rev_comp.reserve(length());
+
+ std::string rc_map = create_reverse_map();
+
+ // reverse and convert
+ seq_string::const_reverse_iterator seq_i;
+ seq_string::const_reverse_iterator seq_end = seq->rend();
+ for(seq_i = seq->rbegin();
+ seq_i != seq_end;
+ ++seq_i)
{
- table_i = (int) at(seq_i);
- rev_comp += conversionTable[table_i];
+ rev_comp.append(1, rc_map[*seq_i]);
}
-
- return rev_comp;
+ return Sequence(rev_comp, alphabet);
}
void Sequence::set_fasta_header(std::string header_)
return "";
}
-void Sequence::set_sequence(const std::string& s)
+const Alphabet& Sequence::get_alphabet() const
+{
+ return get_alphabet(alphabet);
+}
+
+const Alphabet& Sequence::get_alphabet(alphabet_ref alpha) const
{
- set_filtered_sequence(s);
+ 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)
+{
+ set_filtered_sequence(s, a);
}
std::string Sequence::get_sequence() const
}
getline(load_file, file_data_line);
// looks like the sequence is written as a single line
- set_filtered_sequence(file_data_line);
+ set_filtered_sequence(file_data_line, reduced_dna_alphabet);
getline(load_file, file_data_line);
getline(load_file, file_data_line);
if (file_data_line == "<Annotations>")
load_file.close();
}
-std::string
-Sequence::rc_motif(std::string a_motif) const
-{
- std::string rev_comp;
- char conversionTable[257];
- int seq_i, table_i, len;
-
- len = a_motif.length();
- rev_comp.reserve(len);
-
- for(table_i=0; table_i < 256; table_i++)
- {
- conversionTable[table_i] = '~';
- }
- // add end of std::string character for printing out table for testing purposes
- conversionTable[256] = '\0';
-
- // add in the characters for the bases we want to convert (IUPAC)
- conversionTable[(int)'A'] = 'T';
- conversionTable[(int)'T'] = 'A';
- conversionTable[(int)'G'] = 'C';
- conversionTable[(int)'C'] = 'G';
- conversionTable[(int)'N'] = 'N';
- conversionTable[(int)'M'] = 'K';
- conversionTable[(int)'R'] = 'Y';
- conversionTable[(int)'W'] = 'W';
- conversionTable[(int)'S'] = 'S';
- conversionTable[(int)'Y'] = 'R';
- conversionTable[(int)'K'] = 'M';
- conversionTable[(int)'V'] = 'B';
- conversionTable[(int)'H'] = 'D';
- conversionTable[(int)'D'] = 'H';
- conversionTable[(int)'B'] = 'V';
-
- // finally, the actual conversion loop
- for(seq_i = len - 1; seq_i >= 0; seq_i--)
- {
- //std::cout << "** i = " << seq_i << " bp = " <<
- table_i = (int) a_motif[seq_i];
- rev_comp += conversionTable[table_i];
- }
-
- //std::cout << "seq: " << a_motif << std::endl;
- //std::cout << "rc: " << rev_comp << std::endl;
-
- return rev_comp;
-}
-
-std::string
-Sequence::motif_normalize(const std::string& a_motif)
-{
- std::string valid_motif;
- int seq_i, len;
-
- len = a_motif.length();
- valid_motif.reserve(len);
-
- // this just upcases IUPAC symbols. Eventually should return an error if non IUPAC is present.
- // current nonIUPAC symbols are omitted, which is not reported atm
- for(seq_i = 0; seq_i < len; seq_i++)
- {
- if ((a_motif[seq_i] == 'a') || (a_motif[seq_i] == 'A'))
- valid_motif += 'A';
- else if ((a_motif[seq_i] == 't') || (a_motif[seq_i] == 'T'))
- valid_motif += 'T';
- else if ((a_motif[seq_i] == 'g') || (a_motif[seq_i] == 'G'))
- valid_motif += 'G';
- else if ((a_motif[seq_i] == 'c') || (a_motif[seq_i] == 'C'))
- valid_motif += 'C';
- else if ((a_motif[seq_i] == 'n') || (a_motif[seq_i] == 'N'))
- valid_motif += 'N';
- else if ((a_motif[seq_i] == 'm') || (a_motif[seq_i] == 'M'))
- valid_motif += 'M';
- else if ((a_motif[seq_i] == 'r') || (a_motif[seq_i] == 'R'))
- valid_motif += 'R';
- else if ((a_motif[seq_i] == 'w') || (a_motif[seq_i] == 'W'))
- valid_motif += 'W';
- else if ((a_motif[seq_i] == 's') || (a_motif[seq_i] == 'S'))
- valid_motif += 'S';
- else if ((a_motif[seq_i] == 'y') || (a_motif[seq_i] == 'Y'))
- valid_motif += 'Y';
- else if ((a_motif[seq_i] == 'k') || (a_motif[seq_i] == 'K'))
- valid_motif += 'G';
- else if ((a_motif[seq_i] == 'v') || (a_motif[seq_i] == 'V'))
- valid_motif += 'V';
- else if ((a_motif[seq_i] == 'h') || (a_motif[seq_i] == 'H'))
- valid_motif += 'H';
- else if ((a_motif[seq_i] == 'd') || (a_motif[seq_i] == 'D'))
- valid_motif += 'D';
- else if ((a_motif[seq_i] == 'b') || (a_motif[seq_i] == 'B'))
- valid_motif += 'B';
- else {
- std::string msg = "Letter ";
- msg += a_motif[seq_i];
- msg += " is not a valid IUPAC symbol";
- throw motif_normalize_error(msg);
- }
- }
- //std::cout << "valid_motif is: " << valid_motif << std::endl;
- return valid_motif;
-}
void Sequence::add_motif(const Sequence& a_motif)
{
}
std::vector<int>
-Sequence::find_motif(const std::string& a_motif) const
+Sequence::find_motif(const Sequence& a_motif) const
{
std::vector<int> motif_match_starts;
- std::string norm_motif_rc;
+ Sequence norm_motif_rc;
motif_match_starts.clear();
+ // std::cout << "motif is: " << norm_motif << std::endl;
- //std::cout << "motif is: " << a_motif << std::endl;
- std::string norm_motif = motif_normalize(a_motif);
- //std::cout << "motif is: " << a_motif << std::endl;
-
- if (norm_motif.size() > 0)
+ if (a_motif.size() > 0)
{
//std::cout << "Sequence: none blank motif\n";
- motif_scan(norm_motif, &motif_match_starts);
+ motif_scan(a_motif, &motif_match_starts);
- norm_motif_rc = rc_motif(a_motif);
+ norm_motif_rc = a_motif.rev_comp();;
// make sure not to do search again if it is a palindrome
- if (norm_motif_rc != norm_motif) {
+ if (norm_motif_rc != a_motif) {
motif_scan(norm_motif_rc, &motif_match_starts);
}
}
return motif_match_starts;
}
-std::vector<int>
-Sequence::find_motif(const Sequence& a_motif) const
-{
- return find_motif(a_motif.get_sequence());
-}
-
void
-Sequence::motif_scan(std::string a_motif, std::vector<int> * motif_match_starts) const
+Sequence::motif_scan(const Sequence& a_motif, std::vector<int> * motif_match_starts) const
{
// if there's no sequence we can't scan for it?
// should this throw an exception?
if (!seq) return;
- std::string::const_iterator seq_c = seq->begin();
- std::string::size_type seq_i;
- int motif_i, motif_len;
+ std::string::size_type seq_i = 0;
+ Sequence::size_type motif_i = 0;
+ Sequence::size_type motif_len = a_motif.length();
+ Sequence::value_type motif_char;
+ Sequence::value_type seq_char;
- //std::cout << "Sequence: motif, seq len = " << sequence.length() << std::endl;
- motif_len = a_motif.length();
-
- //std::cout << "motif_length: " << motif_len << std::endl;
- //std::cout << "RAAARRRRR\n";
-
- motif_i = 0;
-
- //std::cout << "motif: " << a_motif << std::endl;
-
- //std::cout << "Sequence: motif, length= " << length << std::endl;
- seq_i = 0;
- while (seq_i < length())
+ while (seq_i < seq->length())
{
- //std::cout << seq_c[seq_i];
- //std::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
- if (a_motif[motif_i] =='N')
+ motif_char = toupper(a_motif[motif_i]);
+ seq_char = toupper(seq->at(seq_i));
+ if (motif_char =='N')
motif_i++;
- else if (a_motif[motif_i] == seq_c[seq_i])
+ else if (motif_char == seq_char)
motif_i++;
- else if ((a_motif[motif_i] =='M') &&
- ((seq_c[seq_i]=='A') || (seq_c[seq_i]=='C')))
+ else if ((motif_char =='M') && ((seq_char=='A') || (seq_char=='C')))
motif_i++;
- else if ((a_motif[motif_i] =='R') &&
- ((seq_c[seq_i]=='A') || (seq_c[seq_i]=='G')))
+ else if ((motif_char =='R') && ((seq_char=='A') || (seq_char=='G')))
motif_i++;
- else if ((a_motif[motif_i] =='W') &&
- ((seq_c[seq_i]=='A') || (seq_c[seq_i]=='T')))
+ else if ((motif_char =='W') && ((seq_char=='A') || (seq_char=='T')))
motif_i++;
- else if ((a_motif[motif_i] =='S') &&
- ((seq_c[seq_i]=='C') || (seq_c[seq_i]=='G')))
+ else if ((motif_char =='S') && ((seq_char=='C') || (seq_char=='G')))
motif_i++;
- else if ((a_motif[motif_i] =='Y') &&
- ((seq_c[seq_i]=='C') || (seq_c[seq_i]=='T')))
+ else if ((motif_char =='Y') && ((seq_char=='C') || (seq_char=='T')))
motif_i++;
- else if ((a_motif[motif_i] =='K') &&
- ((seq_c[seq_i]=='G') || (seq_c[seq_i]=='T')))
+ else if ((motif_char =='K') && ((seq_char=='G') || (seq_char=='T')))
motif_i++;
- else if ((a_motif[motif_i] =='V') &&
- ((seq_c[seq_i]=='A') || (seq_c[seq_i]=='C') ||
- (seq_c[seq_i]=='G')))
+ else if ((motif_char =='V') &&
+ ((seq_char=='A') || (seq_char=='C') || (seq_char=='G')))
motif_i++;
- else if ((a_motif[seq_i] =='H') &&
- ((seq_c[seq_i]=='A') || (seq_c[seq_i]=='C') ||
- (seq_c[seq_i]=='T')))
+ else if ((motif_char =='H') &&
+ ((seq_char=='A') || (seq_char=='C') || (seq_char=='T')))
motif_i++;
- else if ((a_motif[motif_i] =='D') &&
- ((seq_c[seq_i]=='A') || (seq_c[seq_i]=='G') ||
- (seq_c[seq_i]=='T')))
+ else if ((motif_char =='D') &&
+ ((seq_char=='A') || (seq_char=='G') || (seq_char=='T')))
motif_i++;
- else if ((a_motif[motif_i] =='B') &&
- ((seq_c[seq_i]=='C') || (seq_c[seq_i]=='G') ||
- (seq_c[seq_i]=='T')))
+ else if ((motif_char =='B') &&
+ ((seq_char=='C') || (seq_char=='G') || (seq_char=='T')))
motif_i++;
-
else
{
+ // if a motif doesn't match, erase our current trial and try again
seq_i -= motif_i;
motif_i = 0;
}
// end Nora stuff, now we see if a match is found this pass
if (motif_i == motif_len)
{
- //std::cout << "!!";
annot new_motif;
motif_match_starts->push_back(seq_i - motif_len + 1);
motif_i = 0;
// 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 (*xseq_i != *yseq_i) {
+ if (toupper(*xseq_i) != toupper(*yseq_i)) {
return false;
}
}
return true;
}
+
+bool operator!=(const Sequence& x, const Sequence& y)
+{
+ return not operator==(x, y);
+}
projects / mussa.git / blobdiff