1 /** 2 Work with scaffold graphs. 3 4 Copyright: © 2018 Arne Ludwig <arne.ludwig@posteo.de> 5 License: Subject to the terms of the MIT license, as written in the 6 included LICENSE file. 7 Authors: Arne Ludwig <arne.ludwig@posteo.de> 8 */ 9 module dentist.common.scaffold; 10 11 import dentist.common.alignments : getType; 12 import dentist.util.log; 13 import dentist.util.math : 14 add, 15 bulkAdd, 16 filterEdges, 17 Graph, 18 MissingNodeException, 19 NaturalNumberSet; 20 import std.algorithm : 21 count, 22 equal, 23 filter, 24 fold, 25 joiner, 26 map, 27 minElement, 28 setDifference, 29 sort, 30 sum; 31 import std.array : appender, array; 32 import std.conv : to; 33 import std.functional : binaryFun; 34 import std.range : 35 enumerate, 36 iota, 37 isForwardRange, 38 only, 39 refRange, 40 retro, 41 save, 42 walkLength; 43 import std.typecons : Flag, No, Tuple, Yes; 44 import vibe.data.json : toJson = serializeToJson; 45 46 debug 47 { 48 import std.conv : to; 49 import std.stdio : writeln; 50 } 51 52 53 /// 54 unittest 55 { 56 // contig 1 contig 2 57 // 58 // o o o o 59 // / e1 e2 \ / e4 60 // o -- o ------- o -- o 61 // \ e3 \ 62 // \ \ e5 (strong evidence) 63 // e10 \ ____________ / 64 // \ / e6 \ / 65 // o -- o o -- o o -- o 66 // \ e7 e8 / \ e9 67 // o o o o o o 68 // 69 // contig 5 contig 4 contig 3 70 // 71 alias Payload = int[]; 72 alias J = Join!Payload; 73 alias S = Scaffold!Payload; 74 alias CN = ContigNode; 75 alias CP = ContigPart; 76 auto scaffold = buildScaffold!(concatenatePayloads!Payload, Payload)(5, [ 77 J(CN(1, CP.end), CN(1, CP.post ), [1]), // e1 78 J(CN(1, CP.end), CN(1, CP.post ), [1]), // e1 79 J(CN(2, CP.pre), CN(2, CP.begin), [1]), // e2 80 J(CN(1, CP.end), CN(2, CP.begin), [1]), // e3 81 J(CN(2, CP.end), CN(2, CP.post ), [1]), // e4 82 J(CN(2, CP.end), CN(3, CP.end ), [1, 1]), // e5 83 J(CN(4, CP.end), CN(3, CP.end ), [1]), // e6 84 J(CN(4, CP.end), CN(4, CP.post ), [1]), // e7 85 J(CN(3, CP.pre), CN(3, CP.begin), [1]), // e8 86 J(CN(3, CP.end), CN(3, CP.post ), [1]), // e9 87 J(CN(4, CP.end), CN(1, CP.begin), [1]), // e10 88 ]).discardAmbiguousJoins!Payload(1.5); 89 // 90 // contig 1 contig 2 91 // 92 // o o o o 93 // / e1 e2 \ / e4 94 // o -- o ------- o -- o 95 // e3 96 // 97 // o -- o o -- o o -- o 98 // \ e7 e8 / \ e9 99 // o o o o o o 100 // 101 // contig 5 contig 4 contig 3 102 103 assert(J(CN(1, CP.end), CN(1, CP.post)) in scaffold); // e1 104 assert(J(CN(2, CP.pre), CN(2, CP.begin)) in scaffold); // e2 105 assert(J(CN(1, CP.end), CN(2, CP.begin)) in scaffold); // e3 106 assert(J(CN(2, CP.end), CN(2, CP.post)) in scaffold); // e4 107 assert(J(CN(2, CP.end), CN(3, CP.end)) in scaffold); // e5 108 assert(J(CN(4, CP.end), CN(3, CP.end)) !in scaffold); // e6 109 assert(J(CN(4, CP.end), CN(4, CP.post)) in scaffold); // e7 110 assert(J(CN(3, CP.pre), CN(3, CP.begin)) in scaffold); // e8 111 assert(J(CN(3, CP.end), CN(3, CP.post)) in scaffold); // e9 112 assert(J(CN(4, CP.end), CN(1, CP.begin)) !in scaffold); // e10 113 114 assert(scaffold.get(J(CN(1, CP.end), CN(1, CP.post))).payload == [1, 1]); // e1 115 } 116 117 /// Each contig has four designated parts where joins can start or end. 118 static enum ContigPart : ubyte 119 { 120 /// Designates a transcendent point *before* the contig where 121 /// front extensions end. 122 pre, 123 /// Designates the begin of the contig. 124 begin, 125 /// Designates the end of the contig. 126 end, 127 /// Designates a transcendent point *after* the contig where 128 /// back extensions end. 129 post, 130 } 131 132 bool isReal(in ContigPart contigPart) pure nothrow 133 { 134 return contigPart == ContigPart.begin || contigPart == ContigPart.end; 135 } 136 137 bool isTranscendent(in ContigPart contigPart) pure nothrow 138 { 139 return contigPart == ContigPart.pre || contigPart == ContigPart.post; 140 } 141 142 /** 143 A contig is represented by four `ContigNodes` in the scaffold graph. 144 145 See_Also: ContigPart 146 */ 147 alias ContigNode = Tuple!(size_t, "contigId", ContigPart, "contigPart"); 148 /// Represents a set of joins which conclude possibly several scaffolding 149 /// variants. 150 alias Scaffold(T) = Graph!(ContigNode, void, No.isDirected, T); 151 alias Join(T) = Scaffold!T.Edge; 152 153 Join!T sumPayloads(T)(Join!T[] joins...) nothrow 154 { 155 assert(joins.length > 0); 156 157 auto mergedJoin = joins[0]; 158 159 mergedJoin.payload = joins 160 .map!"a.payload" 161 .sum; 162 163 return mergedJoin; 164 } 165 166 Join!T concatenatePayloads(T)(Join!T[] joins...) nothrow 167 { 168 assert(joins.length > 0); 169 170 auto mergedJoin = joins[0]; 171 172 mergedJoin.payload = joins 173 .map!"a.payload" 174 .joiner 175 .array; 176 177 return mergedJoin; 178 } 179 180 J dontJoin(J)(J j1, J j2) pure nothrow 181 { 182 j1.payload = T.init; 183 184 return j1; 185 } 186 187 /// Returns true iff join is a default edge of the scaffold graph. 188 bool isDefault(J)(in J join) pure nothrow 189 { 190 return join.start.contigPart == ContigPart.begin && 191 join.end.contigPart == ContigPart.end && 192 join.start.contigId == join.end.contigId; 193 } 194 195 /// Returns true iff join is a unknown edge, ie. an edge for unknown sequence 196 /// (`n`s) of the scaffold graph. 197 bool isUnkown(J)(in J join) pure nothrow 198 { 199 return join.start.contigId != join.end.contigId && 200 (join.start.contigPart != join.end.contigPart) && 201 join.start.contigPart.isTranscendent && 202 join.end.contigPart.isTranscendent; 203 } 204 205 /// Returns true iff join is a gap edge of the scaffold graph. 206 bool isGap(J)(in J join) pure nothrow 207 { 208 return join.start.contigId != join.end.contigId && 209 (join.start.contigPart.isReal) && 210 (join.end.contigPart.isReal); 211 } 212 213 /// Returns true iff join is a gap edge and anti-parallel. 214 bool isAntiParallel(J)(in J join) pure nothrow 215 { 216 return join.isGap && join.start.contigPart == join.end.contigPart; 217 } 218 219 /// Returns true iff join is a gap edge and parallel. 220 bool isParallel(J)(in J join) pure nothrow 221 { 222 return join.isGap && join.start.contigPart != join.end.contigPart; 223 } 224 225 /// Returns true iff join is an extension edge of the scaffold graph. 226 bool isExtension(J)(in J join) pure nothrow 227 { 228 return isFrontExtension(join) ^ isBackExtension(join); 229 } 230 231 /// Returns true iff join is a front extension edge of the scaffold graph. 232 bool isFrontExtension(J)(in J join) pure nothrow 233 { 234 return join.start.contigId == join.end.contigId && 235 join.start.contigPart == ContigPart.pre && 236 join.end.contigPart == ContigPart.begin; 237 } 238 239 /// Returns true iff join is a back extension edge of the scaffold graph. 240 bool isBackExtension(J)(in J join) pure nothrow 241 { 242 return join.start.contigId == join.end.contigId && 243 join.start.contigPart == ContigPart.end && 244 join.end.contigPart == ContigPart.post; 245 } 246 247 /// Returns true iff join is a valid edge of the scaffold graph. 248 bool isValid(J)(in J join) pure nothrow 249 { 250 return isDefault(join) ^ isGap(join) ^ isExtension(join) ^ isUnkown(join); 251 } 252 253 /// Build a scaffold graph using `rawJoins`. This creates default edges and 254 /// inserts the rawJoins. 255 Scaffold!T buildScaffold(alias mergeMultiEdges, T, R)(in size_t numReferenceContigs, R rawJoins) 256 { 257 auto scaffold = initScaffold!T(numReferenceContigs) 258 .addJoins!(mergeMultiEdges, T)(rawJoins); 259 260 return scaffold; 261 } 262 263 /// Creates a scaffold with all the default edges. Optionally specify a 264 /// function that produces the payloads. 265 /// 266 /// See_Also: `getDefaultJoin` 267 Scaffold!T initScaffold(alias getPayload, T)(in size_t numReferenceContigs) 268 { 269 auto contigIds = iota(1, numReferenceContigs + 1); 270 auto contigNodes = contigIds 271 .map!(contigId => only( 272 ContigNode(contigId, ContigPart.pre), 273 ContigNode(contigId, ContigPart.begin), 274 ContigNode(contigId, ContigPart.end), 275 ContigNode(contigId, ContigPart.post), 276 )) 277 .joiner 278 .array; 279 auto initialScaffold = Scaffold!T(contigNodes); 280 281 static if (__traits(compiles, getPayload is null) && (getPayload is null)) 282 { 283 alias createDefaultJoin = getDefaultJoin!T; 284 } 285 else 286 { 287 alias createDefaultJoin = getDefaultJoin!(getPayload, T); 288 } 289 290 initialScaffold.bulkAddForce(contigIds.map!createDefaultJoin.array); 291 292 return initialScaffold; 293 } 294 295 /// ditto 296 Scaffold!T initScaffold(T)(in size_t numReferenceContigs) 297 { 298 return initScaffold!(null, T)(numReferenceContigs); 299 } 300 301 /** 302 Get the default join for contigId. Initialize payload with 303 `getPayload(contigId)` if given. 304 */ 305 Join!T getDefaultJoin(T)(size_t contigId) pure nothrow 306 { 307 return Join!T( 308 ContigNode(contigId, ContigPart.begin), 309 ContigNode(contigId, ContigPart.end), 310 ); 311 } 312 313 /// ditto 314 Join!T getDefaultJoin(alias getPayload, T)(size_t contigId) pure nothrow 315 { 316 return Join!T( 317 ContigNode(contigId, ContigPart.begin), 318 ContigNode(contigId, ContigPart.end), 319 getPayload(contigId), 320 ); 321 } 322 323 private Scaffold!T addJoins(alias mergeMultiEdges, T, R)(Scaffold!T scaffold, R rawJoins) if (isForwardRange!R) 324 { 325 version (assert) 326 { 327 foreach (ref join; rawJoins.save) 328 { 329 assert(join.isValid && !join.isDefault); 330 } 331 } 332 333 scaffold.bulkAdd!mergeMultiEdges(rawJoins); 334 335 return scaffold; 336 } 337 338 /// This removes ambiguous gap insertions. 339 Scaffold!T discardAmbiguousJoins(T)(Scaffold!T scaffold, in double bestPileUpMargin) 340 { 341 static if (__traits(compiles, getType(join.payload))) 342 alias joinToJson = (join) => shouldLog(LogLevel.debug_) 343 ? join.toJson 344 : [ 345 "start": join.start.toJson, 346 "end": join.end.toJson, 347 "payload": [ 348 "type": join.payload.getType.to!string.toJson, 349 "readAlignments": shouldLog(LogLevel.debug_) 350 ? join.payload.map!"a[]".array.toJson 351 : toJson(null), 352 ].toJson, 353 ].toJson; 354 else 355 alias joinToJson = (join) => join.toJson; 356 357 auto incidentEdgesCache = scaffold.allIncidentEdges(); 358 359 foreach (contigNode; scaffold.nodes) 360 { 361 assert(!contigNode.contigPart.isTranscendent || incidentEdgesCache[contigNode].length <= 1); 362 363 if (contigNode.contigPart.isReal && incidentEdgesCache[contigNode].length > 2) 364 { 365 auto incidentGapJoins = incidentEdgesCache[contigNode].filter!isGap.array; 366 auto correctGapJoinIdx = incidentGapJoins.findCorrectGapJoin!T(bestPileUpMargin); 367 368 if (correctGapJoinIdx < incidentGapJoins.length) 369 { 370 // Keep correct gap join for diagnostic output 371 auto correctGapJoin = incidentGapJoins[correctGapJoinIdx]; 372 373 // Remove correct gap join from the list 374 incidentGapJoins = incidentGapJoins[0 .. correctGapJoinIdx] ~ 375 incidentGapJoins[correctGapJoinIdx + 1 .. $]; 376 377 logJsonDiagnostic( 378 "info", "removing bad gap pile ups", 379 "sourceContigNode", contigNode.toJson, 380 "correctGapJoin", joinToJson(correctGapJoin), 381 "removedGapJoins", incidentGapJoins.map!joinToJson.array.toJson, 382 ); 383 } 384 else 385 { 386 logJsonDiagnostic( 387 "info", "skipping ambiguous gap pile ups", 388 "sourceContigNode", contigNode.toJson, 389 "removedGapJoins", incidentGapJoins.map!joinToJson.array.toJson, 390 ); 391 } 392 393 // Mark bad/ambiguous pile ups for removal 394 foreach (gapJoin; incidentGapJoins) 395 { 396 gapJoin.payload = T.init; 397 scaffold.add!(scaffold.ConflictStrategy.replace)(gapJoin); 398 } 399 } 400 } 401 402 return removeNoneJoins!T(scaffold); 403 } 404 405 size_t findCorrectGapJoin(T)(Join!T[] incidentGapJoins, in double bestPileUpMargin) 406 { 407 if (incidentGapJoins.length < 2) 408 // There is no ambiguity in the first place. 409 return 0; 410 411 auto pileUpsLengths = incidentGapJoins 412 .map!"a.payload.length" 413 .map!(to!double) 414 .enumerate 415 .array; 416 pileUpsLengths.sort!"a.value > b.value"; 417 auto largestPileUp = pileUpsLengths[0]; 418 auto sndLargestPileUp = pileUpsLengths[1]; 419 420 if (sndLargestPileUp.value * bestPileUpMargin < largestPileUp.value) 421 return largestPileUp.index; 422 else 423 return size_t.max; 424 } 425 426 /// Get join for a stretch of unknown sequence (`n`s). 427 Join!T getUnkownJoin(T)(size_t preContigId, size_t postContigId, T payload) pure nothrow 428 { 429 assert(preContigId != postContigId); 430 return Join!T( 431 ContigNode(preContigId, ContigPart.post), 432 ContigNode(postContigId, ContigPart.pre), 433 payload, 434 ); 435 } 436 437 /// Normalizes unknown joins such that they join contigs or are removed as 438 /// applicable. 439 Scaffold!T normalizeUnkownJoins(T)(Scaffold!T scaffold) 440 { 441 mixin(traceExecution); 442 443 auto numUnkownJoins = scaffold.edges.count!isUnkown; 444 auto newJoins = appender!(Join!T[]); 445 newJoins.reserve(numUnkownJoins); 446 auto removalAcc = appender!(Join!T[]); 447 removalAcc.reserve(numUnkownJoins); 448 auto degreesCache = scaffold.allDegrees(); 449 450 foreach (unkownJoin; scaffold.edges.filter!isUnkown) 451 { 452 auto preContigId = unkownJoin.start.contigId; 453 auto preContigEnd = ContigNode(preContigId, ContigPart.end); 454 auto postContigId = unkownJoin.end.contigId; 455 auto postContigBegin = ContigNode(postContigId, ContigPart.begin); 456 457 bool isPreContigUnconnected = degreesCache[preContigEnd] == 1; 458 bool hasPreContigExtension = scaffold.has(Join!T( 459 preContigEnd, 460 unkownJoin.start, 461 )); 462 bool hasPreContigGap = !isPreContigUnconnected && !hasPreContigExtension; 463 bool isPostContigUnconnected = degreesCache[postContigBegin] == 1; 464 bool hasPostContigExtension = scaffold.has(Join!T( 465 unkownJoin.end, 466 postContigBegin, 467 )); 468 bool hasPostContigGap = !isPostContigUnconnected && !hasPostContigExtension; 469 470 if (isPreContigUnconnected && isPostContigUnconnected) 471 { 472 newJoins ~= Join!T( 473 preContigEnd, 474 postContigBegin, 475 unkownJoin.payload, 476 ); 477 478 unkownJoin.payload = T.init; 479 removalAcc ~= unkownJoin; 480 } 481 else if (isPreContigUnconnected && hasPostContigExtension) 482 { 483 newJoins ~= Join!T( 484 preContigEnd, 485 unkownJoin.end, 486 unkownJoin.payload, 487 ); 488 489 unkownJoin.payload = T.init; 490 removalAcc ~= unkownJoin; 491 } 492 else if (hasPreContigExtension && isPostContigUnconnected) 493 { 494 newJoins ~= Join!T( 495 unkownJoin.start, 496 postContigBegin, 497 unkownJoin.payload, 498 ); 499 500 unkownJoin.payload = T.init; 501 removalAcc ~= unkownJoin; 502 } 503 else if (hasPreContigGap || hasPostContigGap) 504 { 505 unkownJoin.payload = T.init; 506 removalAcc ~= unkownJoin; 507 } 508 } 509 510 scaffold.bulkAddForce(newJoins.data); 511 foreach (noneJoin; removalAcc.data) 512 { 513 scaffold.add!(scaffold.ConflictStrategy.replace)(noneJoin); 514 } 515 516 return removeNoneJoins!T(scaffold); 517 } 518 519 /// 520 unittest 521 { 522 alias J = Join!int; 523 alias S = Scaffold!int; 524 alias CN = ContigNode; 525 alias CP = ContigPart; 526 527 // Case 1: 528 // 529 // o oxxxxo o => o o o o 530 // => 531 // o -- o o -- o => o -- oxxxxxxxxo -- o 532 auto scaffold1 = buildScaffold!(sumPayloads!int, int)(2, [ 533 getUnkownJoin!int(1, 2, 1), 534 ]).normalizeUnkownJoins!int; 535 assert(getDefaultJoin!int(1) in scaffold1); 536 assert(getDefaultJoin!int(2) in scaffold1); 537 assert(J(CN(1, CP.end), CN(2, CP.begin)) in scaffold1); 538 assert(scaffold1.edges.walkLength == 3); 539 540 // Case 2: 541 // 542 // o oxxxxo o => o o xxxo o 543 // \ => / \ 544 // o -- o o -- o => o -- oxxxx o -- o 545 auto scaffold2 = buildScaffold!(sumPayloads!int, int)(2, [ 546 getUnkownJoin!int(1, 2, 1), 547 J(CN(2, CP.pre), CN(2, CP.begin), 1), 548 ]).normalizeUnkownJoins!int; 549 assert(getDefaultJoin!int(1) in scaffold2); 550 assert(getDefaultJoin!int(2) in scaffold2); 551 assert(J(CN(1, CP.end), CN(2, CP.pre)) in scaffold2); 552 assert(J(CN(2, CP.pre), CN(2, CP.begin)) in scaffold2); 553 assert(scaffold2.edges.walkLength == 4); 554 555 // Case 3: 556 // 557 // o oxxxxo o => o o o o 558 // => 559 // o -- o o -- o => o -- o o -- o 560 // | => | 561 // o -- o => o -- o 562 // => 563 // o o => o o 564 auto scaffold3 = buildScaffold!(sumPayloads!int, int)(3, [ 565 getUnkownJoin!int(1, 2, 1), 566 J(CN(2, CP.begin), CN(3, CP.begin), 1), 567 ]).normalizeUnkownJoins!int; 568 assert(getDefaultJoin!int(1) in scaffold3); 569 assert(getDefaultJoin!int(2) in scaffold3); 570 assert(getDefaultJoin!int(3) in scaffold3); 571 assert(J(CN(2, CP.begin), CN(3, CP.begin)) in scaffold3); 572 assert(scaffold3.edges.walkLength == 4); 573 574 // Case 4: 575 // 576 // o oxxxxo o => o oxxx o o 577 // / => / \ 578 // o -- o o -- o => o -- o xxxxo -- o 579 auto scaffold4 = buildScaffold!(sumPayloads!int, int)(2, [ 580 getUnkownJoin!int(1, 2, 1), 581 J(CN(1, CP.end), CN(1, CP.post), 1), 582 ]).normalizeUnkownJoins!int; 583 assert(getDefaultJoin!int(1) in scaffold4); 584 assert(getDefaultJoin!int(2) in scaffold4); 585 assert(J(CN(1, CP.post), CN(2, CP.begin)) in scaffold4); 586 assert(J(CN(1, CP.end), CN(1, CP.post)) in scaffold4); 587 assert(scaffold4.edges.walkLength == 4); 588 589 // Case 5: 590 // 591 // o oxxxxo o 592 // / \ 593 // o -- o o -- o 594 auto scaffold5 = buildScaffold!(sumPayloads!int, int)(2, [ 595 getUnkownJoin!int(1, 2, 1), 596 J(CN(1, CP.end), CN(1, CP.post), 1), 597 J(CN(2, CP.pre), CN(2, CP.begin), 1), 598 ]).normalizeUnkownJoins!int; 599 assert(getDefaultJoin!int(1) in scaffold5); 600 assert(getDefaultJoin!int(2) in scaffold5); 601 assert(getUnkownJoin!int(1, 2, 1) in scaffold5); 602 assert(J(CN(1, CP.end), CN(1, CP.post)) in scaffold5); 603 assert(J(CN(2, CP.pre), CN(2, CP.begin)) in scaffold5); 604 assert(scaffold5.edges.walkLength == 5); 605 606 // Case 6: 607 // 608 // o oxxxxo o => o o o o 609 // / => / 610 // o -- o o -- o => o -- o o -- o 611 // | => | 612 // o -- o => o -- o 613 // => 614 // o o => o o 615 auto scaffold6 = buildScaffold!(sumPayloads!int, int)(3, [ 616 getUnkownJoin!int(1, 2, 1), 617 J(CN(1, CP.end), CN(1, CP.post), 1), 618 J(CN(2, CP.begin), CN(3, CP.begin), 1), 619 ]).normalizeUnkownJoins!int; 620 assert(getDefaultJoin!int(1) in scaffold6); 621 assert(getDefaultJoin!int(2) in scaffold6); 622 assert(getDefaultJoin!int(3) in scaffold6); 623 assert(J(CN(1, CP.end), CN(1, CP.post)) in scaffold6); 624 assert(J(CN(2, CP.begin), CN(3, CP.begin)) in scaffold6); 625 assert(scaffold6.edges.walkLength == 5); 626 627 // Case 7: 628 // 629 // o oxxxxo o => o o o o 630 // => 631 // o -- o o -- o => o -- o o -- o 632 // | => | 633 // o -- o => o -- o 634 // => 635 // o o => o o 636 auto scaffold7 = buildScaffold!(sumPayloads!int, int)(3, [ 637 getUnkownJoin!int(1, 2, 1), 638 J(CN(1, CP.end), CN(3, CP.end), 1), 639 ]).normalizeUnkownJoins!int; 640 assert(getDefaultJoin!int(1) in scaffold7); 641 assert(getDefaultJoin!int(2) in scaffold7); 642 assert(getDefaultJoin!int(3) in scaffold7); 643 assert(J(CN(1, CP.end), CN(3, CP.end)) in scaffold7); 644 assert(scaffold7.edges.walkLength == 4); 645 646 // Case 8: 647 // 648 // o oxxxxo o => o o o o 649 // \ => \ 650 // o -- o o -- o => o -- o o -- o 651 // | => | 652 // o -- o => o -- o 653 // => 654 // o o => o o 655 auto scaffold8 = buildScaffold!(sumPayloads!int, int)(3, [ 656 getUnkownJoin!int(1, 2, 1), 657 J(CN(1, CP.end), CN(3, CP.end), 1), 658 J(CN(2, CP.pre), CN(2, CP.begin), 1), 659 ]).normalizeUnkownJoins!int; 660 assert(getDefaultJoin!int(1) in scaffold8); 661 assert(getDefaultJoin!int(2) in scaffold8); 662 assert(getDefaultJoin!int(3) in scaffold8); 663 assert(J(CN(1, CP.end), CN(3, CP.end)) in scaffold8); 664 assert(J(CN(2, CP.pre), CN(2, CP.begin)) in scaffold8); 665 assert(scaffold8.edges.walkLength == 5); 666 667 // Case 9: 668 // 669 // o oxxxxo o => o o o o 670 // => 671 // o -- o o -- o => o -- o o -- o 672 // | | => | | 673 // o ------ o => o ------ o 674 // => 675 // o o => o o 676 auto scaffold9 = buildScaffold!(sumPayloads!int, int)(3, [ 677 getUnkownJoin!int(1, 2, 1), 678 J(CN(1, CP.end), CN(3, CP.begin), 1), 679 J(CN(2, CP.begin), CN(3, CP.end), 1), 680 ]).normalizeUnkownJoins!int; 681 assert(getDefaultJoin!int(1) in scaffold9); 682 assert(getDefaultJoin!int(2) in scaffold9); 683 assert(getDefaultJoin!int(3) in scaffold9); 684 assert(J(CN(1, CP.end), CN(3, CP.begin)) in scaffold9); 685 assert(J(CN(2, CP.begin), CN(3, CP.end)) in scaffold9); 686 assert(scaffold9.edges.walkLength == 5); 687 } 688 689 /// Determine which kinds of joins are allowed. 690 enum JoinPolicy 691 { 692 /// Only join gaps inside of scaffolds (marked by `n`s in FASTA). 693 scaffoldGaps, 694 /// Join gaps inside of scaffolds (marked by `n`s in FASTA) and try to 695 /// join scaffolds. 696 scaffolds, 697 /// Break input into contigs and re-scaffold everything; maintains scaffold gaps where new 698 /// scaffolds are consistent. 699 contigs, 700 } 701 702 /// Enforce joinPolicy in scaffold. 703 Scaffold!T enforceJoinPolicy(T)(Scaffold!T scaffold, in JoinPolicy joinPolicy) 704 { 705 mixin(traceExecution); 706 707 if (joinPolicy == JoinPolicy.contigs) 708 { 709 // Just do nothing; this is the default mode of operation. 710 return scaffold; 711 } 712 713 alias orderByNodes = Scaffold!T.orderByNodes; 714 715 assert(joinPolicy == JoinPolicy.scaffoldGaps || joinPolicy == JoinPolicy.scaffolds); 716 717 auto allowedJoins = scaffold 718 .edges 719 .filter!isUnkown 720 .map!(join => only( 721 Join!T( 722 ContigNode(join.start.contigId, ContigPart.end), 723 ContigNode(join.start.contigId, ContigPart.post), 724 ), 725 Join!T( 726 ContigNode(join.start.contigId, ContigPart.end), 727 ContigNode(join.end.contigId, ContigPart.begin), 728 ), 729 Join!T( 730 ContigNode(join.end.contigId, ContigPart.pre), 731 ContigNode(join.end.contigId, ContigPart.begin), 732 ), 733 )) 734 .joiner; 735 auto gapJoins = scaffold.edges.filter!isGap; 736 auto forbiddenJoins = setDifference!orderByNodes(gapJoins, allowedJoins).array; 737 738 // NOTE: joins between scaffolds will be re-included into the graph later 739 // if joinPolicy == scaffolds. 740 foreach (forbiddenJoin; forbiddenJoins) 741 { 742 forbiddenJoin.payload = T.init; 743 scaffold.add!(scaffold.ConflictStrategy.replace)(forbiddenJoin); 744 } 745 746 scaffold = removeNoneJoins!T(scaffold); 747 748 if (joinPolicy == JoinPolicy.scaffolds) 749 { 750 scaffold = normalizeUnkownJoins!T(scaffold); 751 752 alias validJoin = (candidate) => scaffold.degree(candidate.start) == 1 753 && scaffold.degree(candidate.end) == 1; 754 foreach (scaffoldJoin; forbiddenJoins.filter!validJoin) 755 { 756 scaffold.add(scaffoldJoin); 757 } 758 759 if (shouldLog(LogLevel.info)) 760 { 761 forbiddenJoins = forbiddenJoins.filter!(j => !validJoin(j)).array; 762 } 763 } 764 765 logJsonInfo("forbiddenJoins", forbiddenJoins 766 .filter!(gapJoin => scaffold.degree(gapJoin.start) == 1 && scaffold.degree(gapJoin.end) == 1) 767 .map!(join => [ 768 "start": join.start.toJson, 769 "end": join.end.toJson, 770 ]) 771 .array 772 .toJson); 773 774 return scaffold; 775 } 776 777 /// Enforce joinPolicy in scaffold. 778 Scaffold!T removeExtensions(T)(Scaffold!T scaffold) 779 { 780 auto extensionJoins = scaffold.edges.filter!isExtension; 781 782 foreach (extensionJoin; extensionJoins) 783 { 784 extensionJoin.payload = T.init; 785 scaffold.add!(scaffold.ConflictStrategy.replace)(extensionJoin); 786 } 787 788 return removeNoneJoins!T(scaffold); 789 } 790 791 /// Remove marked edges from the graph. This always keeps the default edges. 792 Scaffold!T removeNoneJoins(T)(Scaffold!T scaffold) 793 { 794 scaffold.filterEdges!(noneJoinFilter!T); 795 796 return scaffold; 797 } 798 799 bool noneJoinFilter(T)(Join!T join) 800 { 801 return isDefault(join) || join.payload != T.init; 802 } 803 804 /// Remove extension edges were they coincide with a gap edge combining their 805 /// payload. This is intended to build pile ups with all reads that contribute 806 /// to each gap. 807 Scaffold!T mergeExtensionsWithGaps(alias mergePayloads, T)(Scaffold!T scaffold) 808 { 809 foreach (contigNode; scaffold.nodes) 810 { 811 assert(!contigNode.contigPart.isTranscendent || scaffold.degree(contigNode) <= 1); 812 assert(scaffold.degree(contigNode) <= 3); 813 814 if (contigNode.contigPart.isReal && scaffold.degree(contigNode) == 3) 815 { 816 auto incidentJoins = scaffold.incidentEdges(contigNode) 817 .filter!(j => !isDefault(j)).array; 818 assert(incidentJoins.length == 2); 819 // The gap join has real `contigPart`s on both ends. 820 int gapJoinIdx = incidentJoins[0].target(contigNode).contigPart.isReal ? 0 : 1; 821 auto gapJoin = incidentJoins[gapJoinIdx]; 822 auto extensionJoin = incidentJoins[$ - gapJoinIdx - 1]; 823 824 gapJoin.payload = binaryFun!mergePayloads(gapJoin.payload, extensionJoin.payload); 825 extensionJoin.payload = T.init; 826 827 scaffold.add!(scaffold.ConflictStrategy.replace)(gapJoin); 828 scaffold.add!(scaffold.ConflictStrategy.replace)(extensionJoin); 829 } 830 } 831 832 return removeNoneJoins!T(scaffold); 833 } 834 835 /// 836 unittest 837 { 838 alias J = Join!int; 839 alias S = Scaffold!int; 840 alias CN = ContigNode; 841 alias CP = ContigPart; 842 // contig 1 contig 2 843 // 844 // o o o o 845 // / e1 e2 \ / e4 846 // o -- o ------- o -- o 847 // e3 848 // 849 // o -- o o -- o o -- o 850 // \ e7 e8 / \ e9 851 // o o o o o o 852 // 853 // contig 5 contig 4 contig 3 854 // 855 auto scaffold = buildScaffold!(sumPayloads!int, int)(5, [ 856 J(CN(1, CP.end), CN(1, CP.post ), 1), // e1 857 J(CN(1, CP.end), CN(1, CP.post ), 1), // e1 858 J(CN(2, CP.pre), CN(2, CP.begin), 1), // e2 859 J(CN(1, CP.end), CN(2, CP.begin), 1), // e3 860 J(CN(2, CP.end), CN(2, CP.post ), 1), // e4 861 J(CN(4, CP.end), CN(4, CP.post ), 1), // e7 862 J(CN(3, CP.pre), CN(3, CP.begin), 1), // e8 863 J(CN(3, CP.end), CN(3, CP.post ), 1), // e9 864 ]).mergeExtensionsWithGaps!("a + b", int); 865 866 assert(J(CN(1, CP.end), CN(1, CP.post)) !in scaffold); // e1 867 assert(J(CN(2, CP.pre), CN(2, CP.begin)) !in scaffold); // e2 868 assert(J(CN(1, CP.end), CN(2, CP.begin)) in scaffold); // e3 869 assert(J(CN(2, CP.end), CN(2, CP.post)) in scaffold); // e4 870 assert(J(CN(4, CP.end), CN(4, CP.post)) in scaffold); // e7 871 assert(J(CN(3, CP.pre), CN(3, CP.begin)) in scaffold); // e8 872 assert(J(CN(3, CP.end), CN(3, CP.post)) in scaffold); // e9 873 874 assert(scaffold.get(J(CN(1, CP.end), CN(2, CP.begin))).payload == 4); // merged 2 * e1 + e2 + e3 875 } 876 877 /** 878 Performs a linear walk through a scaffold graph starting in startNode. 879 A linear walk is a sequence of adjacent joins where no node is visited 880 twice unless the graph is cyclic in which case the first node will appear 881 twice. The implementation requires the graph to have linear components, 882 ie. for every node the degree must be at most two. If the component of 883 `startNode` is cyclic then the walk will in `startNode` and the `isCyclic` 884 flag will be set. 885 886 The direction of the walk can be influenced by giving `firstJoin`. 887 888 **Note:** if one wants to read the `isCyclic` flag it is required to use 889 `std.range.refRange` in most cases. 890 891 Returns: range of joins in the scaffold graph. 892 Throws: MissingNodeException if any node is encountered that is not part 893 of the graph. 894 */ 895 LinearWalk!T linearWalk(T)( 896 Scaffold!T scaffold, 897 ContigNode startNode, 898 Scaffold!T.IncidentEdgesCache incidentEdgesCache = Scaffold!T.IncidentEdgesCache.init, 899 ) 900 { 901 return LinearWalk!T(scaffold, startNode, incidentEdgesCache); 902 } 903 904 /// ditto 905 LinearWalk!T linearWalk(T)( 906 Scaffold!T scaffold, 907 ContigNode startNode, 908 Join!T firstJoin, 909 Scaffold!T.IncidentEdgesCache incidentEdgesCache = Scaffold!T.IncidentEdgesCache.init, 910 ) 911 { 912 return LinearWalk!T(scaffold, startNode, firstJoin, incidentEdgesCache); 913 } 914 915 /// 916 unittest 917 { 918 alias Payload = int; 919 alias J = Join!Payload; 920 alias S = Scaffold!Payload; 921 alias CN = ContigNode; 922 alias CP = ContigPart; 923 // contig 1 contig 2 924 // 925 // o o o o 926 // / e4 927 // o -- o ------- o -- o 928 // e3 929 // 930 // o -- o o -- o o -- o 931 // \ e7 e8 / \ e9 932 // o o o o o o 933 // 934 // contig 5 contig 4 contig 3 935 // 936 auto joins1 = [ 937 J(CN(1, CP.end), CN(2, CP.begin), 1), // e3 938 J(CN(2, CP.end), CN(2, CP.post ), 1), // e4 939 J(CN(4, CP.end), CN(4, CP.post ), 1), // e7 940 J(CN(3, CP.pre), CN(3, CP.begin), 1), // e8 941 J(CN(3, CP.end), CN(3, CP.post ), 1), // e9 942 ]; 943 auto scaffold1 = buildScaffold!(sumPayloads!Payload, Payload)(5, joins1); 944 auto scaffold1Cache = scaffold1.allIncidentEdges(); 945 auto walks1 = [ 946 [ 947 getDefaultJoin!Payload(1), 948 joins1[0], 949 getDefaultJoin!Payload(2), 950 joins1[1], 951 ], 952 [ 953 getDefaultJoin!Payload(4), 954 joins1[2], 955 ], 956 [ 957 joins1[3], 958 getDefaultJoin!Payload(3), 959 joins1[4], 960 ], 961 ]; 962 963 alias getWalkStart = (walk) => walk[0].source(walk[0].getConnectingNode(walk[1])); 964 965 foreach (walk; walks1) 966 { 967 auto reverseWalk = walk.retro.array; 968 auto computedWalk = linearWalk!Payload(scaffold1, getWalkStart(walk)); 969 auto computedReverseWalk = linearWalk!Payload(scaffold1, getWalkStart(reverseWalk)); 970 971 assert(equal(walk[], refRange(&computedWalk))); 972 assert(!computedWalk.isCyclic); 973 assert(equal(reverseWalk[], refRange(&computedReverseWalk))); 974 assert(!computedReverseWalk.isCyclic); 975 } 976 977 foreach (walk; walks1) 978 { 979 auto reverseWalk = walk.retro.array; 980 auto computedWalk = linearWalk!Payload(scaffold1, getWalkStart(walk), scaffold1Cache); 981 auto computedReverseWalk = linearWalk!Payload(scaffold1, getWalkStart(reverseWalk), scaffold1Cache); 982 983 assert(equal(walk[], refRange(&computedWalk))); 984 assert(!computedWalk.isCyclic); 985 assert(equal(reverseWalk[], refRange(&computedReverseWalk))); 986 assert(!computedReverseWalk.isCyclic); 987 } 988 989 // contig 1 contig 2 990 // 991 // o o o o 992 // 993 // e1 994 // o -- o ------- o -- o 995 // \_________________/ 996 // e2 997 // 998 auto joins2 = [ 999 J(CN(1, CP.end), CN(2, CP.begin), 1), // e1 1000 J(CN(2, CP.end), CN(1, CP.begin ), 1), // e2 1001 ]; 1002 auto scaffold2 = buildScaffold!(sumPayloads!Payload, Payload)(2, joins2); 1003 auto scaffold2Cache = scaffold2.allIncidentEdges(); 1004 auto walk2 = [ 1005 getDefaultJoin!Payload(1), 1006 joins2[0], 1007 getDefaultJoin!Payload(2), 1008 joins2[1], 1009 ]; 1010 1011 { 1012 auto computedWalk = linearWalk!Payload(scaffold2, getWalkStart(walk2), walk2[0]); 1013 auto reverseWalk2 = walk2.retro.array; 1014 auto computedReverseWalk2 = linearWalk!Payload(scaffold2, 1015 getWalkStart(reverseWalk2), reverseWalk2[0]); 1016 1017 assert(equal(walk2[], refRange(&computedWalk))); 1018 assert(computedWalk.isCyclic); 1019 assert(equal(reverseWalk2[], refRange(&computedReverseWalk2))); 1020 assert(computedWalk.isCyclic); 1021 } 1022 1023 { 1024 auto computedWalk = linearWalk!Payload(scaffold2, getWalkStart(walk2), walk2[0], scaffold2Cache); 1025 auto reverseWalk2 = walk2.retro.array; 1026 auto computedReverseWalk2 = linearWalk!Payload(scaffold2, 1027 getWalkStart(reverseWalk2), reverseWalk2[0], scaffold2Cache); 1028 1029 assert(equal(walk2[], refRange(&computedWalk))); 1030 assert(computedWalk.isCyclic); 1031 assert(equal(reverseWalk2[], refRange(&computedReverseWalk2))); 1032 assert(computedWalk.isCyclic); 1033 } 1034 } 1035 1036 struct LinearWalk(T) 1037 { 1038 alias IncidentEdgesCache = Scaffold!T.IncidentEdgesCache; 1039 static enum emptyIncidentEdgesCache = IncidentEdgesCache.init; 1040 1041 private Scaffold!T scaffold; 1042 private IncidentEdgesCache incidentEdgesCache; 1043 private size_t currentNodeIdx; 1044 private Join!T currentJoin; 1045 private bool isEmpty = false; 1046 private Flag!"isCyclic" isCyclic = No.isCyclic; 1047 private NaturalNumberSet visitedNodes; 1048 1049 /// Start linear walk through a scaffold graph in startNode. 1050 this( 1051 Scaffold!T scaffold, 1052 ContigNode startNode, 1053 IncidentEdgesCache incidentEdgesCache = emptyIncidentEdgesCache, 1054 ) 1055 { 1056 this.scaffold = scaffold; 1057 this.incidentEdgesCache = incidentEdgesCache == emptyIncidentEdgesCache 1058 ? scaffold.allIncidentEdges() 1059 : incidentEdgesCache; 1060 this.currentNode = startNode; 1061 this.visitedNodes.reserveFor(this.scaffold.nodes.length - 1); 1062 this.markVisited(this.currentNodeIdx); 1063 this.popFront(); 1064 } 1065 1066 /// Start linear walk through a scaffold graph in startNode. 1067 this( 1068 Scaffold!T scaffold, 1069 ContigNode startNode, 1070 Join!T firstJoin, 1071 IncidentEdgesCache incidentEdgesCache = emptyIncidentEdgesCache, 1072 ) 1073 { 1074 this.scaffold = scaffold; 1075 this.incidentEdgesCache = incidentEdgesCache == emptyIncidentEdgesCache 1076 ? scaffold.allIncidentEdges() 1077 : incidentEdgesCache; 1078 this.currentNode = startNode; 1079 this.visitedNodes.reserveFor(this.scaffold.nodes.length - 1); 1080 this.markVisited(this.currentNodeIdx); 1081 this.currentJoin = firstJoin; 1082 currentNode = currentJoin.target(currentNode); 1083 this.markVisited(this.currentNodeIdx); 1084 } 1085 1086 void popFront() 1087 { 1088 assert(!empty, "Attempting to popFront an empty LinearWalk"); 1089 assert(scaffold.degree(currentNode) <= 2, "fork in linear walk"); 1090 1091 // If isCyclic is set then the last edge of the cycle is being popped. 1092 // Thus we stop. 1093 if (isCyclic) 1094 { 1095 return endOfWalk(); 1096 } 1097 1098 auto candidateEdges = incidentEdgesCache[currentNode] 1099 .filter!(join => !visitedNodes.has(scaffold.indexOf(join.target(currentNode)))); 1100 bool noSuccessorNodes = candidateEdges.empty; 1101 1102 if (noSuccessorNodes) 1103 { 1104 // Iff the current node has more than one neighbor the graph has 1105 // a cycle. 1106 if (scaffold.degree(currentNode) > 1) 1107 { 1108 assert(scaffold.degree(currentNode) == 2); 1109 1110 return lastEdgeOfCycle(); 1111 } 1112 else 1113 { 1114 return endOfWalk(); 1115 } 1116 } 1117 1118 currentJoin = candidateEdges.front; 1119 currentNode = currentJoin.target(currentNode); 1120 markVisited(currentNodeIdx); 1121 } 1122 1123 @property Join!T front() 1124 { 1125 assert(!empty, "Attempting to fetch the front of an empty LinearWalk"); 1126 return currentJoin; 1127 } 1128 1129 @property bool empty() 1130 { 1131 return isEmpty; 1132 } 1133 1134 private void lastEdgeOfCycle() 1135 { 1136 isCyclic = Yes.isCyclic; 1137 // Find the missing edge. 1138 currentJoin = scaffold 1139 .incidentEdges(currentNode) 1140 .filter!(join => join != currentJoin) 1141 .front; 1142 } 1143 1144 private void endOfWalk() 1145 { 1146 currentNodeIdx = scaffold.nodes.length; 1147 isEmpty = true; 1148 } 1149 1150 private @property ContigNode currentNode() 1151 { 1152 return scaffold.nodes[currentNodeIdx]; 1153 } 1154 1155 private @property void currentNode(ContigNode node) 1156 { 1157 this.currentNodeIdx = this.scaffold.indexOf(node); 1158 } 1159 1160 private void markVisited(size_t nodeIdx) 1161 { 1162 this.visitedNodes.add(nodeIdx); 1163 } 1164 } 1165 1166 /// Get a range of `ContigNode`s where full contig walks should start. 1167 auto scaffoldStarts(T)( 1168 Scaffold!T scaffold, 1169 Scaffold!T.IncidentEdgesCache incidentEdgesCache = Scaffold!T.IncidentEdgesCache.init, 1170 ) 1171 { 1172 static struct ContigStarts 1173 { 1174 alias IncidentEdgesCache = Scaffold!T.IncidentEdgesCache; 1175 static enum emptyIncidentEdgesCache = IncidentEdgesCache.init; 1176 1177 Scaffold!T scaffold; 1178 IncidentEdgesCache incidentEdgesCache; 1179 bool _empty = false; 1180 NaturalNumberSet unvisitedNodes; 1181 ContigNode currentContigStart; 1182 1183 this(Scaffold!T scaffold, IncidentEdgesCache incidentEdgesCache = emptyIncidentEdgesCache) 1184 { 1185 this.scaffold = scaffold; 1186 this.incidentEdgesCache = incidentEdgesCache == emptyIncidentEdgesCache 1187 ? scaffold.allIncidentEdges() 1188 : incidentEdgesCache; 1189 unvisitedNodes.reserveFor(scaffold.nodes.length); 1190 1191 foreach (nodeIdx; iota(scaffold.nodes.length)) 1192 { 1193 unvisitedNodes.add(nodeIdx); 1194 } 1195 1196 popFront(); 1197 } 1198 1199 void popFront() 1200 { 1201 ContigNode walkToEndNode(ContigNode lastNode, Join!T walkEdge) 1202 { 1203 auto nextNode = walkEdge.target(lastNode); 1204 unvisitedNodes.remove(scaffold.indexOf(nextNode)); 1205 1206 return nextNode; 1207 } 1208 1209 assert(!empty, "Attempting to popFront an empty ContigStarts"); 1210 1211 if (unvisitedNodes.empty) 1212 { 1213 _empty = true; 1214 1215 return; 1216 } 1217 1218 auto unvisitedNodeIdx = unvisitedNodes.minElement(); 1219 auto unvisitedNode = scaffold.nodes[unvisitedNodeIdx]; 1220 auto unvisitedNodeDegree = incidentEdgesCache[unvisitedNode].length; 1221 unvisitedNodes.remove(unvisitedNodeIdx); 1222 1223 // Ignore unconnected nodes. 1224 if (unvisitedNodeDegree > 0) 1225 { 1226 auto endNodes = incidentEdgesCache[unvisitedNode] 1227 .map!(firstEdge => linearWalk!T(scaffold, unvisitedNode, firstEdge, incidentEdgesCache) 1228 .fold!walkToEndNode(unvisitedNode)); 1229 currentContigStart = unvisitedNodeDegree == 1 1230 // If the start node has only one edge it is itself an end node. 1231 ? minElement(endNodes, unvisitedNode) 1232 : minElement(endNodes); 1233 } 1234 else 1235 { 1236 popFront(); 1237 } 1238 } 1239 1240 @property ContigNode front() pure nothrow 1241 { 1242 assert(!empty, "Attempting to fetch the front of an empty ContigStarts"); 1243 return currentContigStart; 1244 } 1245 1246 @property bool empty() const pure nothrow 1247 { 1248 return _empty; 1249 } 1250 } 1251 1252 return ContigStarts(scaffold, incidentEdgesCache); 1253 } 1254 1255 /// 1256 unittest 1257 { 1258 alias Payload = int; 1259 alias J = Join!Payload; 1260 alias S = Scaffold!Payload; 1261 alias CN = ContigNode; 1262 alias CP = ContigPart; 1263 // contig 1 contig 2 1264 // 1265 // o o o o 1266 // / e4 1267 // o -- o ------- o -- o 1268 // e3 1269 // 1270 // o -- o o -- o o -- o 1271 // \ e7 e8 / \ e9 1272 // o o o o o o 1273 // 1274 // contig 5 contig 4 contig 3 1275 // 1276 auto scaffold1 = buildScaffold!(sumPayloads!Payload, Payload)(5, [ 1277 J(CN(1, CP.end), CN(2, CP.begin), 1), // e3 1278 J(CN(2, CP.end), CN(2, CP.post ), 1), // e4 1279 J(CN(4, CP.end), CN(4, CP.post ), 1), // e7 1280 J(CN(3, CP.pre), CN(3, CP.begin), 1), // e8 1281 J(CN(3, CP.end), CN(3, CP.post ), 1), // e9 1282 ]); 1283 1284 assert(equal(scaffoldStarts!Payload(scaffold1), [ 1285 CN(1, CP.begin), 1286 CN(3, CP.pre), 1287 CN(4, CP.begin), 1288 CN(5, CP.begin), 1289 ])); 1290 1291 // contig 1 contig 2 1292 // 1293 // o o o o 1294 // 1295 // e1 1296 // o -- o ------- o -- o 1297 // \_________________/ 1298 // e2 1299 // 1300 auto scaffold2 = buildScaffold!(sumPayloads!Payload, Payload)(2, [ 1301 J(CN(1, CP.end), CN(2, CP.begin), 1), // e1 1302 J(CN(2, CP.end), CN(1, CP.begin ), 1), // e2 1303 ]); 1304 1305 assert(equal(scaffoldStarts!Payload(scaffold2), [ 1306 CN(1, CP.begin), 1307 ])); 1308 }