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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
|
//![header]
/*==========================================================================
SeqAn - The Library for Sequence Analysis
http://www.seqan.de
============================================================================
Copyright (C) 2010
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 3 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
============================================================================
Author: Manuel Holtgrewe <manuel.holtgrewe@fu-berlin.de>
============================================================================
Demonstration on how to initialize a scoring matrix programatically with:
- one of the built-in matrices, here BLOSUM30
- arbitrary values
- a new, built-in matrix.
==========================================================================*/
//![header]
//![includes]
#include <iostream>
#include <seqan/basic.h>
#include <seqan/stream.h> // For printing strings.
#include <seqan/score.h> // The module score.
using namespace seqan;
//![includes]
//![user-defined-matrix]
// Extend SeqAn by a user-define scoring matrix.
namespace seqan {
// We have to create a new specialization of the ScoringMatrix_ class
// for the DNA alphabet. For this, we first create a new tag.
struct UserDefinedMatrix {};
// Then, we specialize the class ScoringMatrix_ for the Dna5 alphabet.
template <>
struct ScoringMatrixData_<int, Dna5, UserDefinedMatrix>
{
enum
{
VALUE_SIZE = ValueSize<Dna5>::VALUE,
TAB_SIZE = VALUE_SIZE * VALUE_SIZE
};
static inline int const * getData()
{
// The user defined data table. In this case, we use the data from BLOSUM-30.
static int const _data[TAB_SIZE] =
{
1, 0, 0, 0, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 1, 0,
0, 0, 0, 0, 0
};
return _data;
}
};
} // namespace seqan
//![user-defined-matrix]
//![show-scoring-matrix]
// Print a scoring scheme matrix to stdout.
template <typename TScoreValue, typename TSequenceValue, typename TSpec>
void showScoringMatrix(Score<TScoreValue, ScoreMatrix<TSequenceValue, TSpec> > const & scoringScheme)
{
// Print top row.
for (unsigned i = 0; i < ValueSize<TSequenceValue>::VALUE; ++i)
std::cout << "\t" << TSequenceValue(i);
std::cout << std::endl;
// Print each row.
for (unsigned i = 0; i < ValueSize<TSequenceValue>::VALUE; ++i)
{
std::cout << TSequenceValue(i);
for (unsigned j = 0; j < ValueSize<TSequenceValue>::VALUE; ++j)
{
std::cout << "\t" << score(scoringScheme, TSequenceValue(i), TSequenceValue(j));
}
std::cout << std::endl;
}
}
//![show-scoring-matrix]
//![main]
int main()
{
// 1. Define type and constants.
//
// Define types for the score value and the scoring scheme.
typedef int TValue;
typedef Score<TValue, ScoreMatrix<Dna5, Default> > TScoringScheme;
// Define our gap scores in some constants.
int const gapOpenScore = -1;
int const gapExtendScore = -1;
// 2. Construct scoring scheme with default/empty matrix.
//
// Construct new scoring scheme, alternatively only give one score
// that is used for both opening and extension.
TScoringScheme scoringScheme(gapExtendScore, gapOpenScore);
// 3. Fill the now-existing ScoreMatrix
//
// The scoring scheme now already has a matrix of the size
// ValueSize<Dna5>::VALUE x ValueSize<Dna5>::VALUE which
// we can now fill.
// 3.1 We fill the scoring scheme with the product of the coordinates.
std::cout << std::endl << "Coordinate Products" << std::endl;
for (unsigned i = 0; i < ValueSize<Dna5>::VALUE; ++i)
{
for (unsigned j = 0; j < ValueSize<Dna5>::VALUE; ++j)
{
setScore(scoringScheme, Dna5(i), Dna5(j), i * j);
}
}
showScoringMatrix(scoringScheme);
// 3.2 Now, we fill it with the user defined matrix above.
std::cout << "User defined matrix (also Dna5 scoring matrix)..." << std::endl;
setDefaultScoreMatrix(scoringScheme, UserDefinedMatrix());
showScoringMatrix(scoringScheme);
// 4. Show our user-defined Dna5 scoring matrix.
std::cout << "User DNA scoring scheme..." << std::endl;
Score<TValue, ScoreMatrix<Dna5, UserDefinedMatrix> > userScoringSchemeDna;
showScoringMatrix(userScoringSchemeDna);
return 0;
}
//![main]
|
