1+ #include < algorithm>
2+
3+ #include " SequenceAlignment.h"
4+
5+ SequenceAlignment::SequenceAlignment ():
6+ m_numScores(0 ),
7+ m_gapScore(0 )
8+ {
9+ }
10+
11+ SequenceAlignment::~SequenceAlignment ()
12+ {
13+ }
14+
15+ void SequenceAlignment::clear ()
16+ {
17+ m_charToIndexMap.clear ();
18+ m_numScores = 0 ;
19+ m_gapScore = 0 ;
20+ m_scores.clear ();
21+ }
22+
23+ bool SequenceAlignment::read (std::istream& in)
24+ {
25+ clear ();
26+
27+ // read alphabet
28+ std::string alphabet;
29+ in >> m_alphabet;
30+
31+ uint32_t idx = 0 ;
32+ for (auto c : m_alphabet)
33+ {
34+ m_charToIndexMap.insert (std::make_pair (c, idx++));
35+ }
36+
37+ m_numScores = m_alphabet.size ();
38+ m_scores.resize (m_numScores * m_numScores, 0 );
39+
40+ for (uint32_t i = 0 ; i < m_numScores; ++i)
41+ {
42+ for (uint32_t j = 0 ; j < m_numScores; ++j)
43+ {
44+ in >> m_scores[i * m_numScores + j];
45+ }
46+ }
47+ in >> m_gapScore;
48+
49+ in >> m_sequnces[0 ] >> m_sequnces[1 ];
50+
51+ return true ;
52+ }
53+
54+ bool SequenceAlignment::run (
55+ uint32_t &scoreNeedlemanWunsch,
56+ std::string (&seq)[2])
57+ {
58+ uint32_t size[2 ] = { m_sequnces[0 ].size () + 1 , m_sequnces[1 ].size () + 1 };
59+ std::vector<char > nwScore (size[0 ] * size[1 ], 0 );
60+ for (uint32_t i = 0 ; i < size[1 ]; ++i)
61+ {
62+ nwScore[i] = (i + 1 ) * m_gapScore;
63+ }
64+ for (uint32_t i = 0 ; i < size[0 ]; ++i)
65+ {
66+ nwScore[i * size[1 ]] = (i + 1 ) * m_gapScore;
67+ }
68+
69+ uint32_t tmpScore[3 ];
70+ for (uint32_t i = 1 ; i < size[0 ]; ++i)
71+ {
72+ for (uint32_t j = 1 ; j < size[1 ]; ++j)
73+ {
74+ uint32_t idx1 = m_charToIndexMap[m_sequnces[0 ][i - 1 ]];
75+ uint32_t idx2 = m_charToIndexMap[m_sequnces[1 ][j - 1 ]];
76+ tmpScore[0 ] = nwScore[(i - 1 ) * size[1 ] + (j - 1 )] + m_scores[idx1 * m_numScores + idx2];
77+ tmpScore[1 ] = nwScore[i * size[1 ] + (j - 1 )] + m_gapScore;
78+ tmpScore[2 ] = nwScore[(i - 1 ) * size[1 ] + j] + m_gapScore;
79+ nwScore[i * size[1 ] + j] = std::min (tmpScore[0 ], std::min (tmpScore[1 ], tmpScore[2 ]));
80+ }
81+ }
82+
83+ scoreNeedlemanWunsch = nwScore[size[0 ] * size[1 ] - 1 ];
84+
85+ for (uint32_t i = size[0 ] - 1 , j = size[1 ] - 1 ; (i > 0 ) && (j > 0 );)
86+ {
87+ uint32_t idx1 = m_charToIndexMap[m_sequnces[0 ][i - 1 ]];
88+ uint32_t idx2 = m_charToIndexMap[m_sequnces[1 ][j - 1 ]];
89+ tmpScore[0 ] = nwScore[(i - 1 ) * size[1 ] + (j - 1 )] + m_scores[idx1 * m_numScores + idx2];
90+ tmpScore[1 ] = nwScore[i * size[1 ] + (j - 1 )] + m_gapScore;
91+ tmpScore[2 ] = nwScore[(i - 1 ) * size[1 ] + j] + m_gapScore;
92+ if (nwScore[i * size[1 ] + j] == tmpScore[0 ])
93+ {
94+ seq[0 ].push_back (m_sequnces[0 ][i - 1 ]);
95+ seq[1 ].push_back (m_sequnces[1 ][j - 1 ]);
96+ --i, --j;
97+ }
98+ else if (nwScore[i * size[1 ] + j] == tmpScore[1 ])
99+ {
100+ seq[0 ].push_back (' ' );
101+ seq[1 ].push_back (m_sequnces[1 ][j - 1 ]);
102+ --j;
103+ }
104+ else if (nwScore[i * size[1 ] + j] == tmpScore[2 ])
105+ {
106+ seq[0 ].push_back (m_sequnces[0 ][i - 1 ]);
107+ seq[1 ].push_back (' ' );
108+ --i;
109+ }
110+ }
111+ std::reverse (seq[0 ].begin (), seq[0 ].end ());
112+ std::reverse (seq[1 ].begin (), seq[1 ].end ());
113+
114+ return true ;
115+ }
0 commit comments