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: Tobias Rausch <rausch@embl.de> 33 // ========================================================================== 34 35 #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_FRAGMENT_H_ 36 #define SEQAN_INCLUDE_SEQAN_ALIGN_FRAGMENT_H_ 37 38 namespace seqan 39 { 40 41 ////////////////////////////////////////////////////////////////////////////// 42 // Fragment Specs 43 ////////////////////////////////////////////////////////////////////////////// 44 45 /*! 46 * @class ExactFragment 47 * @extends Fragment 48 * @headerfile <seqan/align.h> 49 * @brief A type for ungapped, pairwise segment matches. 50 * 51 * @signature template <[typename TSize[, typename TSpec]]> 52 * class Fragment<TSize, ExactFragment<TSpec> >; 53 * 54 * @tparam TSize The size type of the underlying sequence. Default: <tt>Size<CharString>Type</tt>. 55 * @tparam TSpec Specializing type. Default: <tt>ExactFragment<></tt>. 56 */ 57 58 template<typename TSpec = Default> 59 struct ExactFragment; 60 61 62 /*! 63 * @class ExactReversableFragment 64 * @extends Fragment 65 * @headerfile <seqan/align.h> 66 * @brief A type for ungapped, pairwise segment matches that maybe in reverse orientation. 67 * 68 * Compared to the @link ExactFragment @endlink <tt>ExactReversableFragment</tt> is a specialized type of @link Fragment @endlink. A @link 69 * ExactReversableFragment @endlink stores an additional bool value to indicate whether a match is in reverse 70 * orientation or not. 71 * 72 * @signature template <[typename TSize[, typename TSpec]]> 73 * class Fragment<TSize, ExactReversableFragment<TSpec> >; 74 * 75 * @tparam TSize The size type of the underlying sequence. Default: <tt>Size<CharString>Type</tt>. 76 * @tparam TSpec Specializing type. Default: <tt>ExactFragment<></tt>. 77 */ 78 79 template<typename TSpec = Default> 80 struct ExactReversableFragment; 81 82 83 ////////////////////////////////////////////////////////////////////////////// 84 // Default Fragment is the exact one 85 ////////////////////////////////////////////////////////////////////////////// 86 87 /*! 88 * @class Fragment 89 * @headerfile <seqan/align.h> 90 * @brief A type for pairwise segment matches. 91 * 92 * @signature template <[typename TSize[, typename TSpec]]> 93 * class Fragment; 94 * 95 * @tparam TSize The size type of the underlying sequence. Default: <tt>Size<CharString>Type</tt>. 96 * @tparam TSpec Specializing type. Default: <tt>ExactFragment<></tt>. 97 * 98 * @section Examples 99 * 100 * @code{.cpp} 101 * // Construct fragment. 102 * unsigned seqId1 = 0, beg1 = 0, seqId2 = 32, beg2 = 42, len = 33; 103 * Fragment<> fragment(seqId1, beg1, seqId2, beg2, len); 104 * 105 * // Update fragment's properties. 106 * fragmentBegin(fragment, 0) = 10; 107 * fragmentBegin(fragment, 1) = 10; 108 * sequenceId(fragment, 0) = 33; 109 * sequenceId(fragment, 1) = 44; 110 * fragmentLength(fragment) += 42; 111 * @endcode 112 */ 113 114 115 template<typename TSize = typename Size<String<char> >::Type, typename TSpec = ExactFragment<> > 116 class Fragment; 117 118 119 ////////////////////////////////////////////////////////////////////////////// 120 // Size Metafunction 121 ////////////////////////////////////////////////////////////////////////////// 122 123 template<typename TSize, typename TSpec> 124 struct Size<Fragment<TSize, TSpec> > { 125 typedef TSize Type; 126 }; 127 128 129 template<typename TSize, typename TSpec> 130 struct Size<Fragment<TSize, TSpec> const> { 131 typedef TSize Type; 132 }; 133 134 135 ////////////////////////////////////////////////////////////////////////////// 136 // Exact Fragment 137 ////////////////////////////////////////////////////////////////////////////// 138 139 140 template<typename TSize, typename TSpec> 141 class Fragment<TSize, ExactFragment<TSpec> > { 142 public: 143 typedef typename Id<Fragment>::Type TId; 144 145 TId seqId1; 146 TSize begin1; 147 TId seqId2; 148 TSize begin2; 149 TSize len; 150 151 /*! 152 * @fn ExactFragment::Fragment 153 * @brief Constructor. 154 * 155 * @signature Fragment::Fragment(); 156 * @signature Fragment::Fragment(seqID1, beg1, seqID2, beg2, l); 157 * 158 * @param[in] seqID1 ID of the first sequence. Type: <tt>TId</tt>. 159 * @param[in] beg1 Begin position of segment match in first sequence. Type: <tt>TSize</tt>. 160 * @param[in] seqID2 ID of the second sequence. Type: <tt>TId</tt>. 161 * @param[in] beg2 Begin position of segment match in second sequence. Type: <tt>TSize</tt>. 162 * @param[in] l The length of the segment match. Type: <tt>TSize</tt>. 163 */ 164 165 Fragment() : seqId1(0), begin1(0), seqId2(0), begin2(0), len(0) {} 166 167 Fragment(TId sqId1, TSize beg1, TId sqId2, TSize beg2, TSize l) : 168 seqId1(sqId1), begin1(beg1), seqId2(sqId2), begin2(beg2), len(l) 169 {} 170 171 }; 172 173 template<typename TSize, typename TSpec> 174 inline bool 175 operator==(Fragment<TSize, ExactFragment<TSpec> > const & left, 176 Fragment<TSize, ExactFragment<TSpec> > const & right) 177 { 178 return (left.seqId1 == right.seqId1 && 179 left.begin1 == right.begin1 && 180 left.seqId2 == right.seqId2 && 181 left.begin2 == right.begin2 && 182 left.len == right.len); 183 } 184 185 template<typename TSize, typename TSpec> 186 inline bool 187 operator<(Fragment<TSize, ExactFragment<TSpec> > const & left, 188 Fragment<TSize, ExactFragment<TSpec> > const & right) 189 { 190 if (left.seqId1 < right.seqId1) 191 return true; 192 if (left.seqId1 > right.seqId1) 193 return false; 194 if (left.begin1 < right.begin1) 195 return true; 196 if (left.begin1 > right.begin1) 197 return false; 198 if (left.seqId2 < right.seqId2) 199 return true; 200 if (left.seqId2 > right.seqId2) 201 return false; 202 if (left.begin2 < right.begin2) 203 return true; 204 if (left.begin2 > right.begin2) 205 return false; 206 if (left.len < right.len) 207 return true; 208 // if (left.len > right.len) 209 // return false; 210 return false; 211 } 212 213 ////////////////////////////////////////////////////////////////////////////// 214 // Exact Fragment that is a forward or reverse match 215 ////////////////////////////////////////////////////////////////////////////// 216 217 218 template<typename TSize, typename TSpec> 219 class Fragment<TSize, ExactReversableFragment<TSpec> > { 220 public: 221 typedef typename Id<Fragment>::Type TId_; 222 223 TId_ seqId1; 224 TSize begin1; 225 TId_ seqId2; 226 TSize begin2; 227 TSize len; 228 bool reversed; 229 230 /*! 231 * @fn ExactReversableFragment::Fragment 232 * @brief Constructor. 233 * 234 * @signature Fragment::Fragment(); 235 * @signature Fragment::Fragment(seqID1, beg1, seqID2, beg2, l[, reversed]); 236 * 237 * @param[in] seqID1 ID of the first sequence. Type: <tt>TId</tt>. 238 * @param[in] beg1 Begin position of segment match in first sequence. Type: <tt>TSize</tt>. 239 * @param[in] seqID2 ID of the second sequence. Type: <tt>TId</tt>. 240 * @param[in] beg2 Begin position of segment match in second sequence. Type: <tt>TSize</tt>. 241 * @param[in] l The length of the segment match. Type: <tt>TSize</tt>. 242 * @param[in] reversed A bool; <tt>true</tt> if the segments match in reverse orientation, <tt>false</tt> otherwise. 243 */ 244 245 246 Fragment() : seqId1(0), begin1(0), seqId2(0), begin2(0), len(0), reversed(false) {} 247 248 Fragment(TId_ sqId1, TSize beg1, TId_ sqId2, TSize beg2, TSize l) : 249 seqId1(sqId1), begin1(beg1), seqId2(sqId2), begin2(beg2), len(l), reversed(false) 250 {} 251 252 Fragment(TId_ sqId1, TSize beg1, TId_ sqId2, TSize beg2, TSize l, bool rev) : 253 seqId1(sqId1), begin1(beg1), seqId2(sqId2), begin2(beg2), len(l), reversed(rev) 254 {} 255 }; 256 257 template<typename TSize, typename TSpec> 258 inline bool 259 operator==(Fragment<TSize, ExactReversableFragment<TSpec> > const & left, 260 Fragment<TSize, ExactReversableFragment<TSpec> > const & right) 261 { 262 return (left.seqId1 == right.seqId1 && 263 left.begin1 == right.begin1 && 264 left.seqId2 == right.seqId2 && 265 left.begin2 == right.begin2 && 266 left.len == right.len && 267 left.reversed == right.reversed); 268 } 269 270 template<typename TSize, typename TSpec> 271 inline bool 272 operator<(Fragment<TSize, ExactReversableFragment<TSpec> > const & left, 273 Fragment<TSize, ExactReversableFragment<TSpec> > const & right) 274 { 275 if (left.seqId1 < right.seqId1) 276 return true; 277 if (left.seqId1 > right.seqId1) 278 return false; 279 if (left.begin1 < right.begin1) 280 return true; 281 if (left.begin1 > right.begin1) 282 return false; 283 if (left.seqId2 < right.seqId2) 284 return true; 285 if (left.seqId2 > right.seqId2) 286 return false; 287 if (left.begin2 < right.begin2) 288 return true; 289 if (left.begin2 > right.begin2) 290 return false; 291 if (left.len < right.len) 292 return true; 293 if (left.len > right.len) 294 return false; 295 if (left.reversed < right.reversed) 296 return true; 297 // if (left.reversed > right.reversed) 298 // return false; 299 return false; 300 } 301 302 ////////////////////////////////////////////////////////////////////////////// 303 304 /*! 305 * @fn Fragment#label 306 * @brief Access to the Fragment's label. 307 * 308 * @signature TInfix label(frag, stringSet, seqID); 309 * 310 * @param[in] frag The Fragment to query. 311 * @param[in] stringSet The @link StringSet @endlink with the sequences. 312 * @param[in] seqID The id of the sequence for which the label should be retrieved. 313 */ 314 315 template<typename TSize, typename TSpec, typename TStringSet, typename TVal> 316 inline typename Infix<typename Value<TStringSet>::Type>::Type 317 label(Fragment<TSize, TSpec> const& f, 318 TStringSet& str, 319 TVal const seqId) 320 { 321 typedef typename Id<Fragment<TSize, TSpec> >::Type TId; 322 return ((TId) seqId == (f.seqId1)) ? infix(getValueById(str, (TId) seqId), f.begin1, f.begin1 + f.len) : infix(getValueById(str, (TId) seqId), f.begin2, f.begin2 + f.len); 323 } 324 325 ////////////////////////////////////////////////////////////////////////////// 326 327 /*! 328 * @fn Fragment#sequenceId 329 * @brief Access to the sequence ID of a fragment. 330 * 331 * @signature TId sequenceId(frag, seqNum); 332 * 333 * @param[in] frag A Fragment. 334 * @param[in] seqNum The sequence number for which the id should be retrieved. Note that @link Fragment @endlink 335 stores information about exactly two sequences which can be accessed with seqNum 0 or 1 but whose 336 ids may differ from their seqNum. 337 * 338 * @return TId Reference to the sequence fragment id member. 339 */ 340 341 template<typename TSize, typename TSpec, typename TVal> 342 inline typename Id<Fragment<TSize, TSpec> >::Type & 343 sequenceId(Fragment<TSize, TSpec> & f, 344 TVal const seqId) 345 { 346 return (seqId == (TVal)0) ? f.seqId1 : f.seqId2; 347 } 348 349 template<typename TSize, typename TSpec, typename TVal> 350 inline typename Id<Fragment<TSize, TSpec> >::Type const & 351 sequenceId(Fragment<TSize, TSpec> const & f, 352 TVal const seqId) 353 { 354 return (seqId == (TVal)0) ? f.seqId1 : f.seqId2; 355 } 356 357 ////////////////////////////////////////////////////////////////////////////// 358 359 /*! 360 * @fn Fragment#fragmentBegin 361 * @brief Return fragment begin. 362 * 363 * @signature TSize fragmentBegin(frag, seqId); 364 * 365 * @param[in] frag The Fragment to query. 366 * @param[in] seqId The id of the sequence to get the begin for. 367 * 368 * @return TSize Reference to the fragment begin position member. 369 */ 370 371 template<typename TSize, typename TSpec, typename TVal> 372 inline TSize& 373 fragmentBegin(Fragment<TSize, TSpec> const& f, 374 TVal const seqId) 375 { 376 typedef typename Id<Fragment<TSize, TSpec> >::Type TId; 377 return ((TId) seqId == f.seqId1) ? const_cast<TSize&>(f.begin1) : const_cast<TSize&>(f.begin2); 378 } 379 380 ////////////////////////////////////////////////////////////////////////////// 381 382 template<typename TSize, typename TSpec, typename TVal> 383 inline TSize& 384 fragmentLength(Fragment<TSize, TSpec> const& f, 385 TVal const) 386 { 387 return const_cast<TSize&>(f.len); 388 } 389 390 ////////////////////////////////////////////////////////////////////////////// 391 392 /*! 393 * @fn Fragment#fragmentLength 394 * @brief Return fragment begin. 395 * 396 * @signature TSize fragmentLength(frag); 397 * 398 * @param[in] frag The Fragment to query for its length. 399 * 400 * @return TSize Reference to the Fragment's length. 401 */ 402 403 template<typename TSize, typename TSpec> 404 inline TSize& 405 fragmentLength(Fragment<TSize, TSpec> const& f) 406 { 407 return const_cast<TSize&>(f.len); 408 } 409 410 ////////////////////////////////////////////////////////////////////////////// 411 412 template<typename TSize, typename TSpec, typename TId1, typename TPosition1, typename TId2, typename TPosition2> 413 inline void 414 getProjectedPosition(Fragment<TSize, ExactFragment<TSpec> > const& f, 415 TId1 const seqId, 416 TPosition1 const pos, 417 TId2& seqId2, 418 TPosition2& pos2) 419 { 420 typedef typename Id<Fragment<TSize, TSpec> >::Type TId; 421 422 if ((TId) seqId == f.seqId1) { 423 SEQAN_ASSERT((TPosition1)f.begin1<=pos); 424 SEQAN_ASSERT(pos - f.begin1 < f.len) ; 425 pos2 = f.begin2 + (pos - f.begin1); 426 seqId2 = f.seqId2; 427 return; 428 } else { 429 SEQAN_ASSERT((TPosition1)f.begin2<=pos); 430 SEQAN_ASSERT(pos - f.begin2 < f.len); 431 pos2 = f.begin1 + (pos - f.begin2); 432 seqId2 = f.seqId1; 433 return; 434 } 435 } 436 437 438 ////////////////////////////////////////////////////////////////////////////// 439 440 template<typename TSize, typename TSpec, typename TValue, typename TId1, typename TPosition1, typename TId2, typename TPosition2> 441 inline void 442 getProjectedPosition(Fragment<TSize, ExactFragment<TSpec> > const& f, 443 TValue seg_num, 444 TId1 const seqId, 445 TPosition1 const pos, 446 TId2& seqId2, 447 TPosition2& pos2) 448 { 449 (void) seqId; // When compiled without assertions. 450 SEQAN_ASSERT((seg_num == 0 && seqId == f.seqId1) || (seg_num == 1 && seqId == f.seqId2)); 451 452 if (seg_num == 0) { 453 SEQAN_ASSERT((TPosition1)f.begin1<=pos); 454 SEQAN_ASSERT(pos - f.begin1 < f.len) ; 455 pos2 = f.begin2 + (pos - f.begin1); 456 seqId2 = f.seqId2; 457 return; 458 } else { 459 SEQAN_ASSERT((TPosition1)f.begin2<=pos); 460 SEQAN_ASSERT(pos - f.begin2 < f.len); 461 pos2 = f.begin1 + (pos - f.begin2); 462 seqId2 = f.seqId1; 463 return; 464 } 465 } 466 467 468 469 ////////////////////////////////////////////////////////////////////////////// 470 471 template<typename TSize, typename TSpec, typename TId1, typename TPosition1, typename TId2, typename TPosition2> 472 inline void 473 getProjectedPosition(Fragment<TSize, ExactReversableFragment<TSpec> > const& f, 474 TId1 const seqId, 475 TPosition1 const pos, 476 TId2& seqId2, 477 TPosition2& pos2) 478 { 479 typedef typename Id<Fragment<TSize, TSpec> >::Type TId; 480 481 if ((TId) seqId == f.seqId1) { 482 SEQAN_ASSERT((TPosition1)f.begin1<=pos); 483 SEQAN_ASSERT(pos - f.begin1 < f.len) ; 484 if (f.reversed) pos2 = (f.begin2 + f.len - 1) - (pos - f.begin1); 485 else pos2 = f.begin2 + (pos - f.begin1); 486 seqId2 = f.seqId2; 487 return; 488 } else { 489 SEQAN_ASSERT((TPosition1)f.begin2<=pos); 490 SEQAN_ASSERT(pos - f.begin2 < f.len); 491 if (f.reversed) pos2 = (f.begin1 + f.len - 1) - (pos - f.begin2); 492 else pos2 = f.begin1 + (pos - f.begin2); 493 seqId2 = f.seqId1; 494 return; 495 } 496 } 497 498 499 ///////////////////////////////////////////////////////////// 500 501 template<typename TSize, typename TSpec, typename TValue, typename TId1, typename TPosition1, typename TId2, typename TPosition2> 502 inline void 503 getProjectedPosition(Fragment<TSize, ExactReversableFragment<TSpec> > const& f, 504 TValue seg_num, 505 TId1 const seqId, 506 TPosition1 const pos, 507 TId2& seqId2, 508 TPosition2& pos2) 509 { 510 (void) seqId; // When compiled without assertions. 511 SEQAN_ASSERT((seg_num == 0 && seqId==f.seqId1) || (seg_num == 1 && seqId==f.seqId2)); 512 513 if (seg_num == 0) { 514 SEQAN_ASSERT((TPosition1)f.begin1<=pos); 515 SEQAN_ASSERT(pos - f.begin1 < f.len) ; 516 if (f.reversed) pos2 = (f.begin2 + f.len - 1) - (pos - f.begin1); 517 else pos2 = f.begin2 + (pos - f.begin1); 518 seqId2 = f.seqId2; 519 return; 520 } else { 521 SEQAN_ASSERT((TPosition1)f.begin2<=pos); 522 SEQAN_ASSERT(pos - f.begin2 < f.len); 523 if (f.reversed) pos2 = (f.begin1 + f.len - 1) - (pos - f.begin2); 524 else pos2 = f.begin1 + (pos - f.begin2); 525 seqId2 = f.seqId1; 526 return; 527 } 528 } 529 530 ////////////////////////////////////////////////////////////////////////////// 531 532 /*! 533 * @fn ExactReversableFragment#isReversed 534 * @brief Return whether a segment match is in reverse orientation. 535 * 536 * @signature bool isReversed(frag); 537 * 538 * @param[in] frag The Fragment to query for reverseness. 539 * 540 * @return bool <tt>true</tt> if the fragment is reversed and <tt>false</tt> otherwise. 541 */ 542 543 template<typename TSize, typename TSpec> 544 inline bool 545 isReversed(Fragment<TSize, ExactReversableFragment<TSpec> > const& f) 546 { 547 return f.reversed; 548 } 549 550 // Compare lexicographically as tuple. 551 552 template<typename TSize, typename TSpec> 553 inline bool operator>(Fragment<TSize, ExactFragment<TSpec> > const & lhs, 554 Fragment<TSize, ExactFragment<TSpec> > const & rhs) 555 { 556 if (lhs.seqId1 > rhs.seqId1) 557 return true; 558 if (lhs.seqId1 == rhs.seqId1 && lhs.begin1 > rhs.begin1) 559 return true; 560 if (lhs.seqId1 == rhs.seqId1 && lhs.begin1 == rhs.begin1 && lhs.seqId2 > rhs.seqId2) 561 return true; 562 if (lhs.seqId1 == rhs.seqId1 && lhs.begin1 == rhs.begin1 && lhs.seqId2 == rhs.seqId2 && lhs.begin2 > rhs.begin2) 563 return true; 564 if (lhs.seqId1 == rhs.seqId1 && lhs.begin1 == rhs.begin1 && lhs.seqId2 == rhs.seqId2 && lhs.begin2 == rhs.begin2 && lhs.len > rhs.len) 565 return true; 566 return false; 567 } 568 569 } // namespace seqan 570 571 #endif // #ifndef SEQAN_INCLUDE_SEQAN_ALIGN_FRAGMENT_H_ 572