1 // ==========================================================================
2 // SeqAn - The Library for Sequence Analysis
3 // ==========================================================================
4 // Copyright (c) 2006-2018, Knut Reinert, FU Berlin
5 // All rights reserved.
6 //
7 // Redistribution and use in source and binary forms, with or without
8 // modification, are permitted provided that the following conditions are met:
9 //
10 // * Redistributions of source code must retain the above copyright
11 // notice, this list of conditions and the following disclaimer.
12 // * Redistributions in binary form must reproduce the above copyright
13 // notice, this list of conditions and the following disclaimer in the
14 // documentation and/or other materials provided with the distribution.
15 // * Neither the name of Knut Reinert or the FU Berlin nor the names of
16 // its contributors may be used to endorse or promote products derived
17 // from this software without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 // ARE DISCLAIMED. IN NO EVENT SHALL KNUT REINERT OR THE FU BERLIN BE LIABLE
23 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
25 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
26 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27 // LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28 // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
29 // DAMAGE.
30 //
31 // ==========================================================================
32 // Author: Manuel Holtgrewe <manuel.holtgrewe@fu-berlin.de>
33 // ==========================================================================
34 // Interface functions for banded local alignment.
35 // ==========================================================================
36
37 #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_LOCAL_ALIGNMENT_BANDED_H_
38 #define SEQAN_INCLUDE_SEQAN_ALIGN_LOCAL_ALIGNMENT_BANDED_H_
39
40 namespace seqan {
41
42 // ============================================================================
43 // Forwards
44 // ============================================================================
45
46 // ============================================================================
47 // Tags, Classes, Enums
48 // ============================================================================
49
50 // ============================================================================
51 // Metafunctions
52 // ============================================================================
53
54 // ============================================================================
55 // Functions
56 // ============================================================================
57
58 // ----------------------------------------------------------------------------
59 // Function localAlignment() [banded, Align]
60 // ----------------------------------------------------------------------------
61
62 template <typename TSequence, typename TAlignSpec, typename TScoreValue, typename TScoreSpec, typename TTag>
localAlignment(Align<TSequence,TAlignSpec> & align,Score<TScoreValue,TScoreSpec> const & scoringScheme,int lowerDiag,int upperDiag,TTag const & tag)63 TScoreValue localAlignment(Align<TSequence, TAlignSpec> & align,
64 Score<TScoreValue, TScoreSpec> const & scoringScheme,
65 int lowerDiag,
66 int upperDiag,
67 TTag const & tag)
68 {
69 typedef Align<TSequence, TAlignSpec> TAlign;
70 typedef typename Size<TAlign>::Type TSize;
71 typedef typename Position<TAlign>::Type TPosition;
72 typedef TraceSegment_<TPosition, TSize> TTraceSegment;
73 typedef AlignConfig2<DPLocal, DPBandConfig<BandOn>, FreeEndGaps_<> > TAlignConfig2;
74
75 SEQAN_ASSERT_EQ(length(rows(align)), 2u);
76
77 String<TTraceSegment> trace;
78 DPScoutState_<Default> dpScoutState;
79 TScoreValue res = _setUpAndRunAlignment(trace, dpScoutState, source(row(align, 0)), source(row(align, 1)),
80 scoringScheme, TAlignConfig2(lowerDiag, upperDiag), tag);
81
82 _adaptTraceSegmentsTo(row(align, 0), row(align, 1), trace);
83 return res;
84 }
85
86 template <typename TSequence, typename TAlignSpec, typename TScoreValue, typename TScoreSpec>
localAlignment(Align<TSequence,TAlignSpec> & align,Score<TScoreValue,TScoreSpec> const & scoringScheme,int lowerDiag,int upperDiag)87 TScoreValue localAlignment(Align<TSequence, TAlignSpec> & align,
88 Score<TScoreValue, TScoreSpec> const & scoringScheme,
89 int lowerDiag,
90 int upperDiag)
91 {
92 SEQAN_ASSERT(length(rows(align)) == 2u);
93 if (_usesAffineGaps(scoringScheme, source(row(align, 0)), source(row(align, 1))))
94 return localAlignment(align, scoringScheme, lowerDiag, upperDiag, AffineGaps());
95 else
96 return localAlignment(align, scoringScheme, lowerDiag, upperDiag, LinearGaps());
97 }
98
99 // ----------------------------------------------------------------------------
100 // Function localAlignment() [banded, Gaps]
101 // ----------------------------------------------------------------------------
102
103 template <typename TSequenceH, typename TGapsSpecH, typename TSequenceV, typename TGapsSpecV,
104 typename TScoreValue, typename TScoreSpec, typename TTag>
localAlignment(Gaps<TSequenceH,TGapsSpecH> & gapsH,Gaps<TSequenceV,TGapsSpecV> & gapsV,Score<TScoreValue,TScoreSpec> const & scoringScheme,int lowerDiag,int upperDiag,TTag const & tag)105 TScoreValue localAlignment(Gaps<TSequenceH, TGapsSpecH> & gapsH,
106 Gaps<TSequenceV, TGapsSpecV> & gapsV,
107 Score<TScoreValue, TScoreSpec> const & scoringScheme,
108 int lowerDiag,
109 int upperDiag,
110 TTag const & tag)
111 {
112 typedef typename Size<TSequenceH>::Type TSize;
113 typedef typename Position<TSequenceH>::Type TPosition;
114 typedef TraceSegment_<TPosition, TSize> TTraceSegment;
115 typedef AlignConfig2<DPLocal, DPBandConfig<BandOn>, FreeEndGaps_<> > TAlignConfig2;
116
117 String<TTraceSegment> trace;
118 DPScoutState_<Default> dpScoutState;
119 TScoreValue res = _setUpAndRunAlignment(trace, dpScoutState, source(gapsH), source(gapsV), scoringScheme,
120 TAlignConfig2(lowerDiag, upperDiag), tag);
121 _adaptTraceSegmentsTo(gapsH, gapsV, trace);
122 return res;
123 }
124
125 template <typename TSequenceH, typename TGapsSpecH, typename TSequenceV, typename TGapsSpecV,
126 typename TScoreValue, typename TScoreSpec>
localAlignment(Gaps<TSequenceH,TGapsSpecH> & gapsH,Gaps<TSequenceV,TGapsSpecV> & gapsV,Score<TScoreValue,TScoreSpec> const & scoringScheme,int lowerDiag,int upperDiag)127 TScoreValue localAlignment(Gaps<TSequenceH, TGapsSpecH> & gapsH,
128 Gaps<TSequenceV, TGapsSpecV> & gapsV,
129 Score<TScoreValue, TScoreSpec> const & scoringScheme,
130 int lowerDiag,
131 int upperDiag)
132 {
133 if (_usesAffineGaps(scoringScheme, source(gapsH), source(gapsV)))
134 return localAlignment(gapsH, gapsV, scoringScheme, lowerDiag, upperDiag, AffineGaps());
135 else
136 return localAlignment(gapsH, gapsV, scoringScheme, lowerDiag, upperDiag, LinearGaps());
137 }
138
139 // ----------------------------------------------------------------------------
140 // Function localAlignment() [banded, Graph<Alignment<> >]
141 // ----------------------------------------------------------------------------
142
143 template <typename TStringSet, typename TCargo, typename TGraphSpec, typename TScoreValue, typename TScoreSpec,
144 typename TTag>
localAlignment(Graph<Alignment<TStringSet,TCargo,TGraphSpec>> & alignmentGraph,Score<TScoreValue,TScoreSpec> const & scoringScheme,int lowerDiag,int upperDiag,TTag const & tag)145 TScoreValue localAlignment(Graph<Alignment<TStringSet, TCargo, TGraphSpec> > & alignmentGraph,
146 Score<TScoreValue, TScoreSpec> const & scoringScheme,
147 int lowerDiag,
148 int upperDiag,
149 TTag const & tag)
150 {
151 typedef Graph<Alignment<TStringSet, TCargo, TGraphSpec> > TGraph;
152 typedef typename Size<TGraph>::Type TSize;
153 typedef typename Position<TGraph>::Type TPosition;
154 typedef TraceSegment_<TPosition, TSize> TTraceSegment;
155 typedef AlignConfig2<DPLocal, DPBandConfig<BandOn>, FreeEndGaps_<> > TAlignConfig2;
156
157 String<TTraceSegment> trace;
158 DPScoutState_<Default> dpScoutState;
159 TScoreValue res = _setUpAndRunAlignment(trace, dpScoutState, value(stringSet(alignmentGraph), 0),
160 value(stringSet(alignmentGraph), 1), scoringScheme,
161 TAlignConfig2(lowerDiag, upperDiag), tag);
162
163 _adaptTraceSegmentsTo(alignmentGraph, positionToId(stringSet(alignmentGraph), 0),
164 positionToId(stringSet(alignmentGraph), 1), trace);
165 return res;
166 }
167
168 template <typename TStringSet, typename TCargo, typename TGraphSpec,
169 typename TScoreValue, typename TScoreSpec>
localAlignment(Graph<Alignment<TStringSet,TCargo,TGraphSpec>> & alignmentGraph,Score<TScoreValue,TScoreSpec> const & scoringScheme,int lowerDiag,int upperDiag)170 TScoreValue localAlignment(Graph<Alignment<TStringSet, TCargo, TGraphSpec> > & alignmentGraph,
171 Score<TScoreValue, TScoreSpec> const & scoringScheme,
172 int lowerDiag,
173 int upperDiag)
174 {
175 SEQAN_ASSERT(length(stringSet(alignmentGraph)) == 2u);
176
177 if (_usesAffineGaps(scoringScheme, stringSet(alignmentGraph)[0], stringSet(alignmentGraph)[1]))
178 return localAlignment(alignmentGraph, scoringScheme, lowerDiag, upperDiag, AffineGaps());
179 else
180 return localAlignment(alignmentGraph, scoringScheme, lowerDiag, upperDiag, LinearGaps());
181 }
182
183 // ----------------------------------------------------------------------------
184 // Function localAlignment() [banded, String<Fragment<> >]
185 // ----------------------------------------------------------------------------
186
187 // Full interface.
188 template <typename TSize, typename TFragmentSpec, typename TStringSpec, typename TSequence, typename TStringSetSpec,
189 typename TScoreValue, typename TScoreSpec, typename TTag>
localAlignment(String<Fragment<TSize,TFragmentSpec>,TStringSpec> & fragmentString,StringSet<TSequence,TStringSetSpec> const & strings,Score<TScoreValue,TScoreSpec> const & scoringScheme,int lowerDiag,int upperDiag,TTag const & tag)190 TScoreValue localAlignment(String<Fragment<TSize, TFragmentSpec>, TStringSpec> & fragmentString,
191 StringSet<TSequence, TStringSetSpec> const & strings,
192 Score<TScoreValue, TScoreSpec> const & scoringScheme,
193 int lowerDiag,
194 int upperDiag,
195 TTag const & tag)
196 {
197 typedef String<Fragment<TSize, TFragmentSpec>, TStringSpec> TFragments;
198 typedef typename Position<TFragments>::Type TPosition;
199 typedef TraceSegment_<TPosition, TSize> TTraceSegment;
200 typedef AlignConfig2<DPLocal, DPBandConfig<BandOn>, FreeEndGaps_<> > TAlignConfig2;
201
202 String<TTraceSegment> trace;
203 DPScoutState_<Default> dpScoutState;
204 TScoreValue res = _setUpAndRunAlignment(trace, dpScoutState, value(strings, 0), value(strings, 1), scoringScheme,
205 TAlignConfig2(lowerDiag, upperDiag), tag);
206 _adaptTraceSegmentsTo(fragmentString, positionToId(strings, 0), positionToId(strings, 1), trace);
207 return res;
208 }
209
210 template <typename TSize, typename TFragmentSpec, typename TStringSpec,
211 typename TSequence, typename TStringSetSpec,
212 typename TScoreValue, typename TScoreSpec>
localAlignment(String<Fragment<TSize,TFragmentSpec>,TStringSpec> & fragmentString,StringSet<TSequence,TStringSetSpec> const & strings,Score<TScoreValue,TScoreSpec> const & scoringScheme,int lowerDiag,int upperDiag)213 TScoreValue localAlignment(String<Fragment<TSize, TFragmentSpec>, TStringSpec> & fragmentString,
214 StringSet<TSequence, TStringSetSpec> const & strings,
215 Score<TScoreValue, TScoreSpec> const & scoringScheme,
216 int lowerDiag,
217 int upperDiag)
218 {
219 SEQAN_ASSERT(length(strings) == 2u);
220
221 if (_usesAffineGaps(scoringScheme, strings[0], strings[1]))
222 return localAlignment(fragmentString, strings, scoringScheme, lowerDiag, upperDiag, AffineGaps());
223 else
224 return localAlignment(fragmentString, strings, scoringScheme, lowerDiag, upperDiag, LinearGaps());
225 }
226
227 // ============================================================================
228 // Many-vs-Many align interfaces.
229 // ============================================================================
230
231 // ----------------------------------------------------------------------------
232 // Function localAlignment() [banded, SIMD version, GapsH, GapsV]
233 // ----------------------------------------------------------------------------
234
235 template <typename TGapSequenceH, typename TSetSpecH,
236 typename TGapSequenceV, typename TSetSpecV,
237 typename TScoreValue, typename TScoreSpec,
238 typename TAlgoTag>
239 inline auto
localAlignment(StringSet<TGapSequenceH,TSetSpecH> & gapSeqSetH,StringSet<TGapSequenceV,TSetSpecV> & gapSeqSetV,Score<TScoreValue,TScoreSpec> const & scoringScheme,int const lowerDiag,int const upperDiag,TAlgoTag const &)240 localAlignment(StringSet<TGapSequenceH, TSetSpecH> & gapSeqSetH,
241 StringSet<TGapSequenceV, TSetSpecV> & gapSeqSetV,
242 Score<TScoreValue, TScoreSpec> const & scoringScheme,
243 int const lowerDiag,
244 int const upperDiag,
245 TAlgoTag const & /*algoTag*/)
246 {
247 typedef AlignConfig2<DPLocal, DPBandConfig<BandOn>, FreeEndGaps_<> > TAlignConfig2;
248 typedef typename SubstituteAlgoTag_<TAlgoTag>::Type TGapModel;
249
250 return _alignWrapper(gapSeqSetH, gapSeqSetV, scoringScheme, TAlignConfig2(lowerDiag, upperDiag), TGapModel());
251 }
252
253 // ----------------------------------------------------------------------------
254 // Function localAlignment() [banded, SIMD version, StringSet<Align>]
255 // ----------------------------------------------------------------------------
256
257 template <typename TSequence, typename TAlignSpec, typename TSetSpec,
258 typename TScoreValue, typename TScoreSpec,
259 typename TAlgoTag>
260 inline String<TScoreValue>
localAlignment(StringSet<Align<TSequence,TAlignSpec>,TSetSpec> & alignSet,Score<TScoreValue,TScoreSpec> const & scoringScheme,int const lowerDiag,int const upperDiag,TAlgoTag const & algoTag)261 localAlignment(StringSet<Align<TSequence, TAlignSpec>, TSetSpec> & alignSet,
262 Score<TScoreValue, TScoreSpec> const & scoringScheme,
263 int const lowerDiag,
264 int const upperDiag,
265 TAlgoTag const & algoTag)
266 {
267 typedef Align<TSequence, TAlignSpec> TAlign;
268 typedef typename Row<TAlign>::Type TGapSequence;
269
270 StringSet<TGapSequence, Dependent<> > gapSetH;
271 StringSet<TGapSequence, Dependent<> > gapSetV;
272 reserve(gapSetH, length(alignSet));
273 reserve(gapSetV, length(alignSet));
274
275 for (auto & align : alignSet)
276 {
277 appendValue(gapSetH, row(align, 0));
278 appendValue(gapSetV, row(align, 1));
279 }
280
281 return localAlignment(gapSetH, gapSetV, scoringScheme, lowerDiag, upperDiag, algoTag);
282 }
283
284 template <typename TSequence, typename TAlignSpec,
285 typename TScoreValue, typename TScoreSpec>
localAlignment(StringSet<Align<TSequence,TAlignSpec>> & align,Score<TScoreValue,TScoreSpec> const & scoringScheme,int lowerDiag,int upperDiag)286 String<TScoreValue> localAlignment(StringSet<Align<TSequence, TAlignSpec> > & align,
287 Score<TScoreValue, TScoreSpec> const & scoringScheme,
288 int lowerDiag, int upperDiag)
289 {
290 if (_usesAffineGaps(scoringScheme, source(row(align[0], 0)), source(row(align[0], 1))))
291 return localAlignment(align, scoringScheme, lowerDiag, upperDiag, AffineGaps());
292 else
293 return localAlignment(align, scoringScheme, lowerDiag, upperDiag, LinearGaps());
294 }
295
296 // ----------------------------------------------------------------------------
297 // Function localAlignmentScore() [banded, no-simd, TSequence]
298 // ----------------------------------------------------------------------------
299
300 template <typename TSequenceH,
301 typename TSequenceV,
302 typename TScoreValue, typename TScoreSpec,
303 typename TAlgoTag>
SEQAN_FUNC_DISABLE_IF(And<And<Is<ContainerConcept<TSequenceH>>,Is<ContainerConcept<typename Value<TSequenceH>::Type>>>,And<Is<ContainerConcept<TSequenceV>>,Is<ContainerConcept<typename Value<TSequenceV>::Type>>>>,TScoreValue)304 SEQAN_FUNC_DISABLE_IF(And<And<Is<ContainerConcept<TSequenceH>>, Is<ContainerConcept<typename Value<TSequenceH>::Type>>>,
305 And<Is<ContainerConcept<TSequenceV>>, Is<ContainerConcept<typename Value<TSequenceV>::Type>>>
306 >, TScoreValue)
307 localAlignmentScore(TSequenceH const & seqH,
308 TSequenceV const & seqV,
309 Score<TScoreValue, TScoreSpec> const & scoringScheme,
310 int const lowerDiag,
311 int const upperDiag,
312 TAlgoTag const & /*algoTag*/)
313 {
314 typedef AlignConfig2<DPLocal, DPBandConfig<BandOn>, FreeEndGaps_<>, TracebackOff> TAlignConfig2;
315 typedef typename SubstituteAlgoTag_<TAlgoTag>::Type TGapModel;
316
317 DPScoutState_<Default> dpScoutState;
318 String<TraceSegment_<unsigned, unsigned> > traceSegments; // Dummy segments.
319 return _setUpAndRunAlignment(traceSegments, dpScoutState, seqH, seqV, scoringScheme,
320 TAlignConfig2{lowerDiag, upperDiag}, TGapModel{});
321 }
322
323 // ----------------------------------------------------------------------------
324 // Function localAlignmentScore() [banded, Simd, TSequence]
325 // ----------------------------------------------------------------------------
326
327 template <typename TSeqH,
328 typename TSeqV,
329 typename TScoreValue, typename TScoreSpec,
330 typename TAlgoTag>
331 inline
SEQAN_FUNC_ENABLE_IF(And<And<Is<ContainerConcept<TSeqH>>,Is<ContainerConcept<typename Value<TSeqH>::Type>>>,And<Is<ContainerConcept<TSeqV>>,Is<ContainerConcept<typename Value<TSeqV>::Type>>>>,String<TScoreValue>)332 SEQAN_FUNC_ENABLE_IF(And<And<Is<ContainerConcept<TSeqH>>, Is<ContainerConcept<typename Value<TSeqH>::Type>>>,
333 And<Is<ContainerConcept<TSeqV>>, Is<ContainerConcept<typename Value<TSeqV>::Type>>>
334 >, String<TScoreValue>)
335 localAlignmentScore(TSeqH const & stringsH,
336 TSeqV const & stringsV,
337 Score<TScoreValue, TScoreSpec> const & scoringScheme,
338 int const lowerDiag,
339 int const upperDiag,
340 TAlgoTag const & /*algoTag*/)
341 {
342 typedef AlignConfig2<DPLocal, DPBandConfig<BandOn>, FreeEndGaps_<>, TracebackOff> TAlignConfig2;
343 typedef typename SubstituteAlgoTag_<TAlgoTag>::Type TGapModel;
344
345 SEQAN_ASSERT_EQ(length(stringsH), length(stringsV));
346 return _alignWrapper(stringsH, stringsV, scoringScheme, TAlignConfig2{lowerDiag, upperDiag}, TGapModel());
347 }
348
349 // Interface without algorithm tag.
350 template <typename TSequenceH,
351 typename TSequenceV,
352 typename TScoreValue, typename TScoreSpec>
localAlignmentScore(TSequenceH const & seqH,TSequenceV const & seqV,Score<TScoreValue,TScoreSpec> const & scoringScheme,int const lowerDiag,int const upperDiag)353 auto localAlignmentScore(TSequenceH const & seqH,
354 TSequenceV const & seqV,
355 Score<TScoreValue, TScoreSpec> const & scoringScheme,
356 int const lowerDiag,
357 int const upperDiag)
358 {
359 if (scoreGapOpen(scoringScheme) == scoreGapExtend(scoringScheme))
360 return localAlignmentScore(seqH, seqV, scoringScheme, lowerDiag, upperDiag, NeedlemanWunsch());
361 else
362 return localAlignmentScore(seqH, seqV, scoringScheme, lowerDiag, upperDiag, Gotoh());
363 }
364
365 } // namespace seqan
366
367 #endif // #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_LOCAL_ALIGNMENT_BANDED_H_
368