1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
|
#include <range/v3/view/zip.hpp>
#include <seqan3/alignment/scoring/nucleotide_scoring_scheme.hpp>
#include <seqan3/alphabet/nucleotide/dna5.hpp>
#include <seqan3/alphabet/nucleotide/rna15.hpp>
#include <seqan3/core/debug_stream.hpp>
#include <seqan3/range/views/zip.hpp>
int main()
{
using seqan3::operator""_dna5;
using seqan3::operator""_dna15;
using seqan3::operator""_rna15;
// You can score two letters:
seqan3::nucleotide_scoring_scheme scheme; // hamming is default
seqan3::debug_stream << "Score between DNA5 A and G: " << (int) scheme.score('A'_dna5, 'G'_dna5) << "\n"; // == -1
seqan3::debug_stream << "Score between DNA5 A and A: " << (int) scheme.score('A'_dna5, 'A'_dna5) << "\n"; // == 0
// You can also score differenct nucleotides:
scheme.set_simple_scheme(seqan3::match_score{3}, seqan3::mismatch_score{-2});
seqan3::debug_stream << "Score between DNA5 A and RNA15 G: " << (int) scheme.score('A'_dna5, 'G'_rna15) << "\n"; // == -2
seqan3::debug_stream << "Score between DNA5 A and RNA15 A: " << (int) scheme.score('A'_dna5, 'A'_rna15) << "\n"; // == 3
// You can "edit" a given matrix directly:
seqan3::nucleotide_scoring_scheme scheme2; // hamming distance is default
seqan3::debug_stream << "Score between DNA A and G before edit: "
<< (int) scheme2.score('A'_dna15, 'G'_dna15) << "\n"; // == -1
scheme2.score('A'_dna15, 'G'_dna15) = 3;
seqan3::debug_stream << "Score after editing: " << (int) scheme2.score('A'_dna15, 'G'_dna15) << "\n"; // == 3
// You can score two sequences:
std::vector<seqan3::dna15> one = "AGAATA"_dna15;
std::vector<seqan3::dna15> two = "ATACTA"_dna15;
seqan3::nucleotide_scoring_scheme scheme3; // hamming distance is default
int score = 0;
for (auto pair : seqan3::views::zip(one, two))
score += scheme3.score(std::get<0>(pair), std::get<1>(pair));
seqan3::debug_stream << "Score: " << score << "\n"; // == 0 - 1 + 0 - 1 + 0 + 0 = -2
}
|
