//-------+---------+---------+---------+---------+---------+---------+--------= // // File: gapped_extend.c // //---------- // // gapped_extend-- // Support for extending anchors to alignments (gaps allowed). // // The Y-drop variant of dynamic programming is applied to create a gapped // alignment by extending in both directions from each 'anchor point". We use // the term Y-drop to distinguish this from the similar X-drop technique for // ungapped alignments. // // The underlying DP algorithm here is the one shown in figure 4 of reference // [1]. However, where the algorithm in [1] computes the entire DP matrix, we // only compute horizontal slices of the DP matrix, with the bounds of each row // determined by (1) any neighboring alignment segments, and (2) cells scoring // less than Y from the max score. // // Throughout this module, we consider sequence 1 ("target") on the vertical // edge of the dynamic programming matrix, and sequence 2 ("query") horizontal. // Deletions are vertical, insertions horizontal. // // More detailed information can be found in the headers of the functions. // Specifically, the function ydrop_one_sided_align() is the heart of the DP // algorithm, and contains more details. // // References: // [1] Approximate Matching of Regular Expressions. Myers and Miller, Bull. // Math. Biol. 51 (1989), pp. 5-37. // //---------- //---------- // // other files // //---------- #include // standard C stuff #define true 1 #define false 0 #include // standard C string stuff #include // standard C value limit stuff #include "build_options.h" // build options #include "utilities.h" // utility stuff #include "dna_utilities.h" // dna/scoring stuff #include "sequences.h" // sequence stuff #include "segment.h" // segment table management stuff #include "edit_script.h" // alignment edit script stuff #include "diag_hash.h" // diagonals hashing stuff #include "identity_dist.h" // identity distribution stuff #include "coverage_dist.h" // query coverage distribution stuff #include "continuity_dist.h" // query continuity distribution stuff #include "output.h" // alignment outout format stuff #define gapped_extend_owner // (make this the owner of its globals) #include "gapped_extend.h" // interface to this module // debugging defines //#define snoopAnchors // if this is defined, extra code is added to // .. track anchors through the alignment // .. process //#define snoopAnchorToGapped // if this is defined, extra code is added to // .. track anchors through the alignment // .. process (different than snoopAnchors) //#define snoopBlocks // if this is defined, extra code is added to // .. track alignment blocks through the process //#define snoopSubprobs // if this is defined, extra code is added to // .. report alignment sub-problems //#define snoopAlgorithm // if this is defined, extra code is added to // .. track the dynamic programming algorithm //#define snoopAlgorithmTrap // if this is defined, extra code is added to // .. allow trapping the dynamic programming // .. algorithm in a debugger, at a specific // .. cell //#define debugPosA1 15384592 // if defined (and snoopAlgorithm is also), //#define debugPosA2 11569819 // .. only information about extension of this //#define debugPosB1 130 // .. particular seed is output; note that //#define debugPosB2 136 // .. this position is origin-zero //#define snoopTraceback // if this is defined, extra code is added to // .. track the dynamic programming algorithm's // .. traceback process //#define debugHspImmediate // if this is defined, extra code is added to // .. aid debugging of gappily_extend_hsps() //#define snoopBatches // if this is defined, extra code is added to // .. track the processing of HSPs in batches //#define snoopEditScripts // if this is defined, extra code is added to // .. examine the edit scripts created by // .. format_alignment() //#define snoopAlignioInput // if this is defined, extra code is added to // .. examine the alignio data structure, which // .. is used to pass data into ydrop_align() //#define snoopAlignioOutput // if this is defined, extra code is added to // .. examine the results in the alignio data // .. structure, which is used to pass data out // .. of ydrop_align() //#define snoopSpecialHsp // if this is defined, extra code is added to // // .. so that only one particular HSP is // // .. processed //#define specialPosA 439 // if defined (and snoopSpecialHsp is also), //#define specialPosB 16 // .. only this HSP is processed; note that // // .. this position is origin-zero //#define snoopBounds // if this is defined, extra code is added to // .. examine bounding alignments //---------- // // stats to augment crude profiling // //---------- //--- this module's contribution to profiling for the whole program -- #ifndef dbgTiming #define dbg_timing_set_stat(field,val) ; #define dbg_timing_count_stat(field) ; #define dbg_timing_report_stat(field,name) ; #endif // not dbgTiming #ifdef dbgTiming struct { int numExtensions; } gappedExtendTimingStats; #define dbg_timing_set_stat(field,val) (gappedExtendTimingStats.field = val) #define dbg_timing_count_stat(field) ++gappedExtendTimingStats.field #define dbg_timing_report_stat(field,name) fprintf(stderr,"%-26s %d\n",name":",gappedExtendTimingStats.field) #endif // dbgTiming //--- profiling for this module only -- #ifndef dbgTimingGappedExtend #define dbg_timing_gapped_extend_sub(v) ; #define dbg_timing_gapped_extend_add(v) ; #define dbg_timing_gapped_extend_copy(dst,src) ; #define dbg_timing_gapped_extend_report(f,v,s) ; #endif // not dbgTimingGappedExtend #ifdef dbgTimingGappedExtend #define read_clock() microsec_clock() #define clocksPerSec 1000000 u64 microsec_clock (void); // (from lastz.c) s64 debugClockAboveBelow = 0, debugClockLeftRight = 0, debugClockYdropAlign = 0, debugClockYdropOneSidedAlign = 0, debugClockUpdateLrBounds = 0, debugClockNextSweepSeg = 0, debugClockPrevSweepSeg = 0, debugClockUpdateActiveSegs = 0, debugClockFilterActiveSegs = 0; #define dbg_timing_gapped_extend_sub(v) { v -= (s64) read_clock(); } #define dbg_timing_gapped_extend_add(v) { v += (s64) read_clock(); } #define dbg_timing_gapped_extend_copy(dst,src) { dst = src; } #define dbg_timing_gapped_extend_report(f,v,s) { fprintf(f,"%-26s %.3f\n",s":",((float)(v))/clocksPerSec); } #endif // dbgTimingGappedExtend //---------- // // private data // //---------- // miscellany #define negInf negInfinity #define anchorPeakLen 31 #undef min #define min(x,y) ((x)<(y)?(x):(y)) #undef max #define max(x,y) ((x)>(y)?(x):(y)) #define signed_difference(a,b) (((sgnpos)(a))-((sgnpos)(b))) // straight alignment segments typedef struct aliseg { char type; // one of diagSeg, horzSeg, vertSeg unspos b1, b2, e1, e2; struct aliseg* nextSeg; struct aliseg* prevSeg; } aliseg; #define diagSeg 0 #define horzSeg 1 // parallel to sequence 2 #define vertSeg 2 // parallel to sequence 1 // gapped alignment; this starts as a single point anchor, and is then extended // into an alignment typedef struct galign { unspos pos1, pos2; // anchor location or start of alignment (these // .. are origin-zero) unspos end1, end2; // end of alignment (these are inclusive) u64 hspId; // (for debugging) uniquely identifies the // .. hsp that led to this alignment attempt aliseg* firstSeg; // the alignment, in diagonal, vertical and aliseg* lastSeg; // .. horizontal segments alignel* align; // second form of the alignment, for external // .. consumption only struct galign *leftAlign1; // at the alignment's beginning and ending struct galign *rightAlign1; // .. points, we keep pointers to the alignments struct galign *leftAlign2; // .. immediately left and right of the terminus struct galign *rightAlign2; aliseg* leftSeg1; // these correspond to leftAlign1, etc., and aliseg* rightSeg1; // .. are the closest segments to the alignment aliseg* leftSeg2; // .. (e.g. leftSeg1 is a segment in leftAlign1) aliseg* rightSeg2; struct galign* next; // alignments are linked both by increasing struct galign* prev; // .. start-point and by decreasing end-point } galign; // input/output for ydrop_align() typedef struct alignio { // the following must be supplied by caller u8* seq1, *seq2; // target and query sequences u8* rev1, *rev2; // reverse of target and query sequences (NOT // .. reverse complement) unspos len1, len2; // total length of sequences unspos low1, low2; // limits of sub-interval in each sequence; low unspos high1, high2; // .. is the leftmost position allowed; high is // .. one past the rightmost position allowed; // .. both are origin-0 indexes into seq->v[] unspos anchor1, anchor2; // position to start the alignment u64 hspId; // (for debugging) uniquely identifies the // .. hsp that led to this alignment attempt scoreset* scoring; // substitution and gap scores to use score yDrop; // value of Y-dropoff parameter int trimToPeak; // whether y-drop should be trimmed (see // description in ydrop_one_sided_align) tback* tb; // block of memory for alignment trace-back galign* leftAlign; // closest alignment to left of anchor pt galign* rightAlign; // closest alignment to right of anchor pt aliseg* leftSeg; // specific segment of leftAlign aliseg* rightSeg; // specific segment of rightAlign galign* aboveList; // alignments starting above the anchor point galign* belowList; // alignments ending below the anchor point // the following are returned by ydrop_align() score s; // alignment's score unspos start1, start2; // alignment's start on target and query unspos stop1, stop2; // alignment's end on target and query editscript* script; // alignment's edit script } alignio; // dynamic programming structure // // We (try to) declare the basic cell with a power-of-two size; doing so MAY // improve speed on some platforms, presumably because random indexing into the // array of dpCells would involve a shift rather than a multiply; however, if // we can't figure out what size the padding should be (or if we don't need any, // as would be the case for dpCell_padding_sz==16), we don't add any #define unpadded_dpCell_sz (score_sz+score_sz+unspos_sz) #if (unpadded_dpCell_sz == 12) #define dpCell_padding_sz 4 #elif (unpadded_dpCell_sz == 20) #define dpCell_padding_sz 12 #elif (unpadded_dpCell_sz == 24) #define dpCell_padding_sz 8 #endif typedef struct dpCell { score DD, CC; unspos mask; // mask out grid points that are part of previous // .. alignments #ifdef dpCell_padding_sz char padding[dpCell_padding_sz]; #endif } dpCell; typedef struct dpMatrix { dpCell* p; u32 len; } dpMatrix; // linked list of active segments typedef struct activeseg { aliseg* seg; unspos x; // column position where segment intersects sweep row; // .. this is relative to the current slice of the DP // .. matrix, e.g. ranging from LY to RY in the context // .. of the alignment sweep unspos lastRow; // (also relative to the alignment sweep) char type; // one of diagSeg, horzSeg, vertSeg char filter; struct activeseg* next; } activeseg; #ifdef snoopSubprobs seq* snoopSubprobsSeq1; seq* snoopSubprobsSeq2; #endif // snoopSubprobs // segment batches-- partitioning of a segment table typedef struct segbatch { u32 start; // index (into a segment table) of the first // .. entry in a batch u32 end; // index (into a segment table) of the first // .. entry NOT in a batch (i.e. the one after // .. the last entry). partition* part; // sequence partition that "contains" this // .. batch; this can be NULL if we aren't // .. dealing with partitions } segbatch; typedef struct sbtable { u32 size; // the number of entries allocated for batch[] u32 len; // the number of batches (the number of entries // .. actually used) segbatch batch[1]; // the batch table (variable-length array) } sbtable; #define sbtable_bytes(size) (sizeof(sbtable) + (((size)-1)*sizeof(segbatch))) // flag related to reprting of "truncated alignments" static int haveReportedTruncation = false; //---------- // // prototypes for private functions // //---------- static unspos segment_peak (u8* s1, u8* s2, unspos segLength, scoreset* scoring); static sbtable* batched_segments (segtable* anchors, seqpartition* sp); static galign** init_from_anchors (segtable* anchors, u32 numExtraSlots); static int identical_sequences (seq* seq1, seq* seq2, scoreset* scoring, score* s); static int identical_partitioned_sequences (seq* seq1, seq* seq2); static int identical_partition_of_sequence (seq* seq1, seq* seq2); static score score_identical_partition (seq* seq1, seq* seq2, partition* p1, partition* p2, scoreset* scoring); static score score_identical_partition_of (seq* seq1, seq* seq2, partition* p1, scoreset* scoring); static void ydrop_align (alignio* io); static score ydrop_one_sided_align (alignio* io, int reversed, u8* A, u8* B, unspos M, unspos N, int trimToPeak, editscript** script, unspos* end1, unspos* end2); static void dp_ready (dpMatrix* dynProg, unspos needed); static int msp_left_right (galign* obi, galign* m); static void get_above_below (alignio* io, galign* obi, galign* oed); static void align_left_right (galign* obi, galign* m); static void insert_align (galign* m, galign** obi, galign** oed); static void update_LR_bounds (int reversed, aliseg** rightSeg, aliseg** leftSeg, galign** rightAlign, galign** leftAlign, unspos row, unspos anchor1, unspos anchor2, sgnpos* L, sgnpos* R, unspos* LY, unspos* RY); static sgnpos next_sweep_seg (int lookRight, aliseg** bp, galign** mp, unspos row, unspos anchor1, unspos anchor2); static sgnpos prev_sweep_seg (int lookRight, aliseg** bp, galign** mp, unspos row, unspos anchor1, unspos anchor2); static void update_active_segs (int reversed, activeseg** active, galign** alignList, dpCell* dp, unspos row, unspos anchor1, unspos anchor2, unspos LY, unspos RY); static void build_active_seg (int reversed, activeseg* act, dpCell* dp, unspos row, unspos anchor1, unspos anchor2, unspos LY, unspos RY); static activeseg* add_new_active (int reversed, activeseg* active, galign* alignList, dpCell* dp, unspos row, unspos anchor1, unspos anchor2, unspos LY, unspos RY); static void filter_active_segs (activeseg** active, int filter); static alignel* format_alignment (alignio* io, galign* m); static void save_seg (galign* m, unspos b1, unspos b2, unspos e1, unspos e2); #ifdef snoopAlignioInput static void dump_alignio_input (FILE* f, alignio* io); #endif // snoopAlignioInput #if ((defined(snoopAlignioOutput)) || (defined(snoopEditScripts))) static void dump_alignio_output(FILE* f, alignio* io); #endif // snoopAlignioOutput OR snoopEditScripts static u64 count_paired_bases (galign* mp); static void warn_for_paired_bases_limit (seq* seq2, u64 maxPairedBases, int overlyPairedKeep); //---------- // // reduce_to_points-- // Convert each segment in a table to it's "peak". The definition of peak is // the midpoint of the highest scoring subsegment of a given length. // //---------- // // Arguments: // seq* seq1: The first sequence. // seq* seq2: The second sequence. // scoreset* scoring: The scoring scheme to use. // segtable* anchors: The segment table to modify. // // Returns: // (nothing) // //---------- void reduce_to_points (seq* seq1, seq* seq2, scoreset* scoring, segtable* anchors) { u32 ix; segment* seg; unspos peak; int reportAnchor = false; if (gapped_extend_dbgShowAnchors) { if (anchors->len == 0) fprintf (stderr, "reduce_to_points: no anchors\n"); else fprintf (stderr, "reduce_to_points: %s anchors\n", ucommatize(anchors->len)); } for (ix=0,seg=anchors->seg ; ixlen ; ix++,seg++) { if (gapped_extend_dbgShowAnchors) { reportAnchor = ((gapped_extend_dbgShowAnchorsHowOften == 0) || (ix % gapped_extend_dbgShowAnchorsHowOften == 0)); if (reportAnchor) fprintf (stderr, "reduce_to_points: reducing [%s] (" unsposSlashFmt " " unsposFmt " diag=" sgnposFmt " score=" scoreFmtSimple ")", commatize(ix), seg->pos1, seg->pos2, seg->length, diagNumber(seg->pos1,seg->pos2), seg->s); } peak = segment_peak (seq1->v+seg->pos1, seq2->v+seg->pos2, seg->length, scoring); seg->pos1 += peak; seg->pos2 += peak; seg->length = 0; if (reportAnchor) fprintf (stderr, " to (" unsposSlashFmt " " unsposFmt ")\n", seg->pos1, seg->pos2, seg->length); } gapped_extend_add_stat (numAnchors, anchors->len); } static unspos segment_peak (u8* s1, u8* s2, unspos segLength, scoreset* scoring) { u8* t1 = s1; u8* t2 = s2; score similarity, best; unspos ix, peak; if (segLength <= anchorPeakLen) peak = segLength / 2; else { similarity = 0; for (ix=0 ; ixsub[*t1++][*t2++]; //fprintf (stderr, "%c %c " scoreFmtSimple " " scoreFmtSimple "\n", // t1[-1], t2[-1], scoring->sub[t1[-1]][t2[-1]], similarity); } best = similarity; peak = anchorPeakLen / 2; for ( ; ixsub[*s1++][*s2++]; similarity += scoring->sub[*t1++][*t2++]; //fprintf (stderr, "%c %c " scoreFmtSimple " " scoreFmtSimple "\n", // t1[-1], t2[-1], scoring->sub[t1[-1]][t2[-1]], similarity); if (similarity > best) { best = similarity; peak = ix - (anchorPeakLen / 2); } } gapped_extend_count_stat (numPeaks); gapped_extend_add_stat (totalPeakScore, best); } return peak; } //---------- // // gapped_extend-- // performed gapped extension given a set of anchor segments. // //---------- // // Arguments: // seq* seq1: The sequence being searched. // u8* rev1: The reverse (NOT reverse complement) of seq1, // .. as a zero-terminated string. // seq* seq2: The sequence being searched for. // u8* rev2: The reverse of the seq2 (analagous to rev1). // int inhibitTrivial: true => don't output the trivial self-alignment. // scoreset* scoring: The scoring scheme to use. // segtable* anchors: The anchor segments. // void* tb: Memory in which to track gapped alignment // .. traceback. // int allBounds: true => bound gapped alignments by *all* gapped // .. extensions of higher-scoring HSPs (a // .. la blastz) // false => bound gapped alignments only by gapped // .. extensions that meet the score // .. threshold // score yDrop: Threshold to stop gapped extensions; if the // .. score drops off by more than yDrop, extension // .. stops // int trimToPeak: Whether y-drop should be trimmed (see // .. description in ydrop_one_sided_align). // sthresh scoreThresh: Minimum score required; gapped alignments are // .. discarded if they score less than this. // u64 maxPairedBases: Maximum number of "paired bases" we'll allow. // .. A paired base is a match or substitution in // .. the DP matrix. If we exceed this limit, we // .. abort processing. However, any gapped // .. alignments we find up to that point are part // .. of what we return. A value of zero inidcates // .. there is no limit. // int overlyPairedWarn: true => write a warning message when a query // .. .. exceeds the maxPairedBases // .. .. threshold. // int overlyPairedKeep: // How we should treat alignments for queries that // .. exceed the maxPairedBases threshold. // .. false => we discard all the alignments // .. true => we output whatever alignments we // .. .. happened to find prior to // .. .. exceeding the limit // // Returns: // A linked list of alignments. The caller is responsible for disposing of // these. The caller must also deallocate other memory; see note (1) below. // //---------- // // Notes: // // (1) This routine calls other routines that allocate long-term memory. The // Calling program is repsonsible for making sure that this memory is // disposed of at the appropriate time (usually at program exit). This // should be done by calling free_segment_batches(). // //---------- //=== stuff for gapped_extend_verbosity === #define debugGappedExtendVerbosity_1 \ if (gapped_extend_verbosity >= 2) \ { \ pos1 = mp->pos1; \ pos2 = mp->pos2; \ \ if (sp1->p != NULL) \ { \ part = lookup_partition (seq1, pos1); \ pos1 += part->sepBefore + 1; \ } \ if (sp2->p != NULL) \ { \ part = lookup_partition (seq2, pos2); \ pos2 += part->sepBefore + 1; \ } \ \ pos1 += seq1->startLoc; \ pos2 += seq2->startLoc; \ \ fprintf (stderr, "processing anchor #%u (of %u)" \ " hspId=" u64Fmt \ " (" unsposSlashFmt ")" \ " " unsposSlashFmt "\n", \ i+1, anchors->len, \ mp->hspId, \ mp->pos1, mp->pos2, \ pos1, pos2); \ } #define debugGappedExtendVerbosity_2 \ if (gapped_extend_verbosity >= 2) \ { \ pos1 = mp->pos1; len1 = mp->end1 - pos1; \ pos2 = mp->pos2; len2 = mp->end2 - pos2; \ \ if (sp1->p != NULL) \ { \ part = lookup_partition (seq1, pos1); \ pos1 += part->sepBefore + 1; \ } \ if (sp2->p != NULL) \ { \ part = lookup_partition (seq2, pos2); \ pos2 += part->sepBefore + 1; \ } \ \ pos1 += seq1->startLoc; \ pos2 += seq2->startLoc; \ \ fprintf (stderr, "alignment block" \ " score=" scoreFmtSimple \ " at (" unsposSlashFmt ") " unsposSlashFmt \ " length " unsposSlashFmt "\n", \ mp->align->s, \ mp->pos1, mp->pos2, pos1, pos2, \ len1, len2); \ } //=== stuff for gapped_extend_dbgShowIdentity === // nota bene: positions reported below are 1-based (not 0-based) #define debugGappedExtendDbgShowIdentity_1 \ if (gapped_extend_dbgShowIdentity) \ printf ("discarding " unsposSlashSFmt " score=" scoreFmtSimple "\n", \ mp->pos1+1, ((seq1->revCompFlags & rcf_rev) != 0)? "-" : "+", \ mp->pos2+1, ((seq2->revCompFlags & rcf_rev) != 0)? "-" : "+", \ mp->align->s); #define debugGappedExtendDbgShowIdentity_2 \ if (gapped_extend_dbgShowIdentity) \ printf ("discarding " unsposSlashSFmt " score=" scoreFmtSimple "\n", \ mp->pos1+1, ((seq1->revCompFlags & rcf_rev) != 0)? "-" : "+", \ mp->pos2+1, ((seq2->revCompFlags & rcf_rev) != 0)? "-" : "+", \ mp->align->s); #define debugGappedExtendDbgShowIdentity_3 \ if (gapped_extend_dbgShowIdentity) \ printf ("keeping " unsposSlashSFmt " score=" scoreFmtSimple "\n", \ mp->pos1+1, ((seq1->revCompFlags & rcf_rev) != 0)? "-" : "+", \ mp->pos2+1, ((seq2->revCompFlags & rcf_rev) != 0)? "-" : "+", \ mp->align->s); //=== stuff for snoopAnchors === #ifndef snoopAnchors #define debugSnoopAnchors_1 ; #endif // not snoopAnchors #ifdef snoopAnchors #define debugSnoopAnchors_1 \ { \ u32 ix; \ fprintf (stderr,"===== msp =====\n"); \ for (ix=0 ; ixlen ; ix++) \ { \ fprintf (stderr,"anchors[%3d,hspid=" u64Fmt "] " unsposSlashFmt " " unsposFmt " " scoreFmt, \ ix, anchors->seg[ix].hspId, \ anchors->seg[ix].pos1+1, \ anchors->seg[ix].pos2+1, \ anchors->seg[ix].length, \ anchors->seg[ix].s); \ fprintf (stderr," msp[%3d] " unsposSlashFmt " " unsposSlashFmt "\n", \ ix, msp[ix]->pos1+1, msp[ix]->pos2+1, \ msp[ix]->end1, msp[ix]->end2); \ } \ } #endif // snoopAnchors //=== stuff for snoopAnchorToGapped === #ifndef snoopAnchorToGapped #define debugSnoopAnchorToGapped_1 ; #define debugSnoopAnchorToGapped_2 ; #define debugSnoopAnchorToGapped_3 ; #define debugSnoopAnchorToGapped_4 ; #define debugSnoopAnchorToGapped_5 ; #endif // not snoopAnchorToGapped #ifdef snoopAnchorToGapped #define debugSnoopAnchorToGapped_1 \ fprintf (stderr, "processing anchor #%u (of %u)\n", \ i+1, anchors->len); #define debugSnoopAnchorToGapped_2 \ fprintf (stderr, "anchor: " unsposSlashFmt " (diag " sgnposFmt ")\n", \ mp->pos1, mp->pos2, diagNumber(mp->pos1,mp->pos2)); #define debugSnoopAnchorToGapped_3 \ fprintf (stderr, " gapped: " unsposSlashFmt " " unsposSlashFmt " " scoreFmt "\n", \ io.start1, io.start2, io.stop1, io.stop2, io.s); #define debugSnoopAnchorToGapped_4 \ fprintf (stderr, "finished processing %u anchors\n", anchors->len); \ fprintf (stderr, " head=%p\n", head); \ for (a=head,i=0 ; a!=NULL ; a=a->next,i++) \ fprintf (stderr, " [%u] %p" \ " %s:" unsposDotsFmt " %s:" unsposDotsFmt "\n", \ i, a, \ seq1->header, a->beg1, a->end1, \ seq2->header, a->beg2, a->end2); #define debugSnoopAnchorToGapped_5 \ fprintf (stderr, "finished processing %u anchors\n", anchors->len); \ fprintf (stderr, " paired threshold was exceeded\n"); #endif // snoopAnchorToGapped //=== stuff for snoopBlocks === #ifndef snoopBlocks #define debugSnoopBlocks_1 ; #define debugSnoopBlocks_2 ; #define debugSnoopBlocks_3 ; #define debugSnoopBlocks_3b ; #define debugSnoopBlocks_4 ; #define debugSnoopBlocks_5 ; #endif // not snoopBlocks #ifdef snoopBlocks #define debugSnoopBlocks_1 \ fprintf (stderr, "===== searching for alignment blocks =====\n"); #define debugSnoopBlocks_2 \ fprintf (stderr, "discarding alignment block [%8p]" \ " b " unsposFmt " " unsposFmt \ " e " unsposFmt " " unsposFmt \ " s " scoreFmtSimple "\n", \ mp, mp->pos1, mp->pos2, mp->end1, mp->end2, mp->align->s); #define debugSnoopBlocks_3 \ fprintf (stderr, "===== collecting alignment blocks =====\n"); #define debugSnoopBlocks_3b \ fprintf (stderr, "===== discarding alignment blocks =====\n"); #define debugSnoopBlocks_4 \ fprintf (stderr, "keeping alignment block [%8p -> %8p]" \ " b " unsposFmt " " unsposFmt \ " e " unsposFmt " " unsposFmt \ " s " scoreFmtSimple "\n", \ mp, mp->align, \ mp->pos1, mp->pos2, mp->end1, mp->end2, mp->align->s); #define debugSnoopBlocks_5 \ fprintf (stderr, "discarding alignment block [%8p]" \ " b " unsposFmt " " unsposFmt \ " e " unsposFmt " " unsposFmt \ " s " scoreFmtSimple "\n", \ mp, mp->pos1, mp->pos2, mp->end1, mp->end2, mp->align->s); #endif // snoopBlocks //=== stuff for snoopSubprobs === #ifndef snoopSubprobs #define debugSnoopSubprobs_1 ; #endif // not snoopSubprobs #ifdef snoopSubprobs #define debugSnoopSubprobs_1 \ snoopSubprobsSeq1 = seq1; \ snoopSubprobsSeq2 = seq2; #endif // snoopSubprobs //=== stuff for snoopBatches === #ifndef snoopBatches #define debugSnoopBatches_1 ; #endif // not snoopBatches #ifdef snoopBatches #define debugSnoopBatches_1 \ if (doHspsInBatches) \ { \ partition* batPart = segBatches->batch[batIx].part; \ fprintf (stderr, "batch[%u] %u..%u", \ batIx, startSegIx, endSegIx-1); \ if (batPart != NULL) \ fprintf (stderr, " " unsposFmt ".." unsposFmt, \ batPart->sepBefore+1, batPart->sepAfter); \ fprintf (stderr, "\n"); \ } #endif // snoopBatches //=== stuff for tryout === #ifndef tryout #define debugTriviality_1 ; #define debugTriviality_2 ; #define debugTriviality_3 ; #define debugTriviality_4 ; #define debugTriviality_5 ; #define debugTriviality_6 ; #endif // not tryout #ifdef tryout #define debugTriviality_1 \ if (gapped_extend_dbgTriviality) \ { \ fprintf (stderr, "trivial?: \"%s\" " unsposFmt " \"%s\" " unsposFmt "\n", \ name1, len1, name2, len2); \ } #define debugTriviality_2 \ if (gapped_extend_dbgTriviality) \ { \ fprintf (stderr, " sequence lengths differ\n"); \ } #define debugTriviality_3 \ if (gapped_extend_dbgTriviality) \ { \ fprintf (stderr, " alignment lengths differ\n"); \ } #define debugTriviality_4 \ if (gapped_extend_dbgTriviality) \ { \ fprintf (stderr, " sequence names differ\n"); \ } #define debugTriviality_5 \ if (gapped_extend_dbgTriviality) \ { \ fprintf (stderr, " alignment content differs\n"); \ } #define debugTriviality_6 \ if (gapped_extend_dbgTriviality) \ { \ fprintf (stderr, " it's trivial!\n"); \ } #endif // tryout //=== stuff for dbgTimingGappedExtend === #ifdef dbgTimingGappedExtend void gapped_extend_timing_report (arg_dont_complain(FILE* f)) { dbg_timing_gapped_extend_report (f, debugClockAboveBelow, "total time in above_below()"); dbg_timing_gapped_extend_report (f, debugClockLeftRight, "total time in left_right()"); dbg_timing_gapped_extend_report (f, debugClockYdropAlign, "total time in ydrop_align()"); dbg_timing_gapped_extend_report (f, debugClockYdropOneSidedAlign, " ydrop_one_sided_align()"); dbg_timing_gapped_extend_report (f, debugClockUpdateLrBounds, " update_lr_bounds()"); dbg_timing_gapped_extend_report (f, debugClockNextSweepSeg, " next_sweep_seg()"); dbg_timing_gapped_extend_report (f, debugClockPrevSweepSeg, " prev_sweep_seg()"); dbg_timing_gapped_extend_report (f, debugClockUpdateActiveSegs, " update_active_segs()"); dbg_timing_gapped_extend_report (f, debugClockFilterActiveSegs, " filter_active_segs()"); } #endif // dbgTimingGappedExtend //=== stuff for snoopEditScripts === #ifndef snoopEditScripts #define debugSnoopEditScripts_1 ; #define debugSnoopEditScripts_2 ; #define debugSnoopEditScripts_3 ; #define debugSnoopEditScripts_4 ; #define debugSnoopEditScripts_5 ; #endif // not snoopEditScripts #ifdef snoopEditScripts #define debugSnoopEditScripts_1 \ fprintf (stderr, "(adding full length trivial alignment)\n"); #define debugSnoopEditScripts_2 \ fprintf (stderr, "(adding full length trivial alignment vs. partition %u)\n", \ trivialPartIx); #define debugSnoopEditScripts_3 \ fprintf (stderr, "(adding trivial alignment for partition %u)\n", \ partIx); #define debugSnoopEditScripts_4 \ { \ char* opName[4] = { "?", "I", "D", "S" }; \ alignel* aa = mp->align; \ u32 opIx, op, rpt; \ fprintf (stderr, " " unsposDotsFmt " vs " unsposDotsFmt \ " (diag " sgnposFmt ")" \ " score=" scoreFmt, \ aa->beg1, aa->end1, aa->beg2, aa->end2, \ diagNumber(aa->beg1,aa->beg2), \ aa->s); \ if (aa->isTrivial) fprintf (stderr, " (trivial)"); \ fprintf (stderr, "\n "); \ for (opIx=0 ; opIxscript->len ; opIx++) \ { \ op = edit_op_operation (aa->script->op[opIx]); \ rpt = edit_op_repeat (aa->script->op[opIx]); \ fprintf (stderr, " %dx%s", rpt, opName[op]); \ } \ fprintf (stderr, "\n"); \ } #define debugSnoopEditScripts_5 \ fprintf (stderr, " (alignment is empty)\n"); #endif // snoopEditScripts //=== stuff for snoopSpecialHsp === #ifndef snoopSpecialHsp #define debugsnoopSpecialHsp_1 ; #endif // not snoopSpecialHsp #ifdef snoopSpecialHsp #define debugsnoopSpecialHsp_1 \ { \ if ((mp->pos1 != specialPosA) || (mp->pos2 != specialPosB)) \ { \ fprintf (stderr," Ignoring msp[%3d] " unsposSlashFmt " " unsposSlashFmt "\n", \ i, mp->pos1+1, mp->pos2+1, mp->end1, mp->end2); \ continue; \ } \ } #endif // snoopSpecialHsp //=== finally, the actual function gapped_extend() === alignel* gapped_extend (seq* seq1, u8* rev1, seq* seq2, u8* rev2, int inhibitTrivial, scoreset* scoring, segtable* anchors, tback* tb, int allBounds, score yDrop, int trimToPeak, sthresh scoreThresh, u64 maxPairedBases, int overlyPairedWarn, int overlyPairedKeep) { seqpartition* sp1 = &seq1->partition; seqpartition* sp2 = &seq2->partition; partition* p1, *p2; galign** msp; alignel* a; score s; alignio io; galign* mp, *mpNext; galign* orderBegInc, *orderEndDec; alignel* head, *last; aliseg* bp, *bq; u32 i; unspos pos1, pos2, len1, len2; partition* part; int partitionedTriviality, delayedCheckForTrivial; int doHspsInBatches; sbtable* segBatches; u32 batIx, startSegIx, endSegIx, partIx; int trivialPartIx = -1; u32 freeSlot; u64 pairedBases = 0, newPairedBases; orderBegInc = orderEndDec = NULL; // (compiler appeasement) if (scoreThresh.t != 'S') suicidef ("gapped_extend can't handle score threshold %s", score_thresh_to_string (&scoreThresh)); // create a gapped alignment table containing one entry for each HSP (plus // an additional slot for each possible trivial self-alignment-- see // note (1) in init_from_anchors); note that batched_segments sorts the // HSPs in descreasing score order (per batch) doHspsInBatches = false; if (gapped_extend_dbgAllowBatches) doHspsInBatches = (sp1->p != NULL); // (seq1 is partitioned) if (!doHspsInBatches) segBatches = batched_segments (anchors, NULL); else segBatches = batched_segments (anchors, sp1); if ((sp1->p == NULL) || (sp2->p == NULL)) msp = init_from_anchors (anchors, 1); else { u32 numExtraSlots = sp1->len; if (sp2->len < numExtraSlots) numExtraSlots = sp2->len; msp = init_from_anchors (anchors, numExtraSlots); } debugSnoopAnchors_1; debugSnoopBlocks_1; debugSnoopSubprobs_1; // set up the "io" block, which is used for communication with lower-level // routines. io.seq1 = seq1->v; io.seq2 = seq2->v; io.rev1 = rev1; io.rev2 = rev2; io.low1 = 0; io.len1 = io.high1 = seq1->len; io.low2 = 0; io.len2 = io.high2 = seq2->len; io.scoring = scoring; io.yDrop = yDrop; io.trimToPeak = trimToPeak; if (tb == NULL) suicide ("gapped_extend was given a NULL traceback pointer."); io.tb = tb; // special case check for trivial alignments of identical sequences; if the // two sequences are identical, we add the trivial self alignment to the // table immediately (so it can prevent some other alignment near the main // diagonal from merging with major portions of the trivial alignment); // note that if trivial alignments aren't desired we will discard them // later, after we've extended all the other anchors // // see note (1) in init_from_anchors-- the gapped alignment table has // extra slots to allow for the addition of these trivial self-alignments; // there is one extra slot for min(nTarget,nQuery), where e.g. nTarget is // the number of partitions in the target (or 1, if it is not partitioned). // // the check for trivialty becomes more complicated if we are dealing with // partitioned sequences // // also note that the more complicated test provides *no guarantee* that // it will prevent other other alignments from merging into the diagonal; // it is possible that two identical sequences will end up with an aligment // that only has part of the main diagonal, and no alignment covering the // complete main diagonal; the author does not consider this a significant // failure because it can only happen when the user does *not* tell the // program she is doing a self-alignment partitionedTriviality = false; delayedCheckForTrivial = ((inhibitTrivial) && ((sp1->p != NULL) || (sp2->p != NULL))); if ((sp1->p != NULL) // (seq1 is partitioned) && (sp2->p == NULL)) // (seq2 is not partitioned) { // $$$ this could be modified to check seq2 vs *every* partition in // $$$ .. seq1, but would require some restructuring of the code; we // $$$ .. don't expect there are several identical sequences in seq1 trivialPartIx = identical_partition_of_sequence (seq1, seq2); partitionedTriviality = (trivialPartIx != -1); delayedCheckForTrivial = (inhibitTrivial) && (!partitionedTriviality); } else if ((sp1->p != NULL) // (seq1 is partitioned) && (sp2->p != NULL)) // (seq2 is partitioned) { // $$$ this could be modified to check every partition in seq1 vs every // $$$ .. partition in seq2, but would require some restructuring of the // $$$ .. code; we don't expect seq2 to be partitioned in normal use; // $$$ .. moreover, this would require *many* more slots in the gapped // $$$ .. alignment table partitionedTriviality = identical_partitioned_sequences (seq1, seq2); delayedCheckForTrivial = (inhibitTrivial) && (!partitionedTriviality); } if (!doHspsInBatches) // empty the list of bounding orderBegInc = orderEndDec = NULL; // .. alignments if ((sp1->p == NULL) // (seq1 is not partitioned) && (sp2->p == NULL) // (seq2 is not partitioned) && (identical_sequences (seq1, seq2, scoring, &s))) { if (doHspsInBatches) suicidef ("internal error, attempt to add trivial self-alignment with batches in play"); freeSlot = anchors->len; mp = msp[freeSlot]; mp->pos1 = mp->pos2 = 0; mp->end1 = mp->end2 = seq1->len-1; mp->leftAlign1 = mp->leftAlign2 = mp->rightAlign1 = mp->rightAlign2 = NULL; mp->leftSeg1 = mp->leftSeg2 = mp->rightSeg1 = mp->rightSeg2 = NULL; mp->firstSeg = NULL; save_seg (mp, mp->pos1, mp->pos2, mp->end1, mp->end2); insert_align(mp, &orderBegInc, &orderEndDec); mp->lastSeg = mp->firstSeg; mp->firstSeg->prevSeg = mp->lastSeg->nextSeg = NULL; a = mp->align = malloc_or_die ("gapped_extend", sizeof(alignel)); a->script = edit_script_new(); edit_script_sub (&a->script, seq1->len); a->beg1 = a->beg2 = 1; a->end1 = a->end2 = seq1->len; a->seq1 = seq1->v; a->seq2 = seq2->v; if (s < scoreThresh.s) a->s = scoreThresh.s; // so it won't be discarded else a->s = s; a->next = NULL; a->isTrivial = true; debugSnoopEditScripts_1; debugSnoopEditScripts_4; } else if ((partitionedTriviality) && (sp2->p == NULL)) // (seq2 is not partitioned) { if (doHspsInBatches) suicidef ("internal error, attempt to add trivial self-alignment with batches in play"); freeSlot = anchors->len; p1 = &sp1->p[trivialPartIx]; s = score_identical_partition_of (seq1, seq2, p1, scoring); mp = msp[freeSlot]; mp->pos1 = p1->sepBefore + 1; mp->pos2 = 0; mp->end1 = p1->sepAfter - 1; mp->end2 = seq2->len - 1; mp->leftAlign1 = mp->leftAlign2 = mp->rightAlign1 = mp->rightAlign2 = NULL; mp->leftSeg1 = mp->leftSeg2 = mp->rightSeg1 = mp->rightSeg2 = NULL; mp->firstSeg = NULL; save_seg (mp, mp->pos1, mp->pos2, mp->end1, mp->end2); insert_align(mp, &orderBegInc, &orderEndDec); mp->lastSeg = mp->firstSeg; mp->firstSeg->prevSeg = mp->lastSeg->nextSeg = NULL; a = mp->align = malloc_or_die ("gapped_extend", sizeof(alignel)); a->script = edit_script_new(); edit_script_sub (&a->script, seq2->len); a->beg1 = p1->sepBefore + 2; a->beg2 = 1; a->end1 = p1->sepAfter; a->end2 = seq2->len; a->seq1 = seq1->v; a->seq2 = seq2->v; if (s < scoreThresh.s) a->s = scoreThresh.s; // so it won't be discarded else a->s = s; a->next = NULL; a->isTrivial = true; debugSnoopEditScripts_2; debugSnoopEditScripts_4; } else if ((partitionedTriviality) && (sp2->p != NULL)) // (seq2 is partitioned) { if (doHspsInBatches) suicidef ("internal error, attempt to add trivial self-alignment with batches in play"); freeSlot = anchors->len; for (partIx=0 ; partIxlen ; partIx++) { p1 = &sp1->p[partIx]; p2 = &sp2->p[partIx]; s = score_identical_partition (seq1, seq2, p1, p2, scoring); mp = msp[freeSlot++]; mp->pos1 = p1->sepBefore + 1; mp->pos2 = p2->sepBefore + 1; mp->end1 = p1->sepAfter - 1; mp->end2 = p2->sepAfter - 1; mp->leftAlign1 = mp->leftAlign2 = mp->rightAlign1 = mp->rightAlign2 = NULL; mp->leftSeg1 = mp->leftSeg2 = mp->rightSeg1 = mp->rightSeg2 = NULL; mp->firstSeg = NULL; save_seg (mp, mp->pos1, mp->pos2, mp->end1, mp->end2); insert_align(mp, &orderBegInc, &orderEndDec); mp->lastSeg = mp->firstSeg; mp->firstSeg->prevSeg = mp->lastSeg->nextSeg = NULL; a = mp->align = malloc_or_die ("gapped_extend", sizeof(alignel)); a->script = edit_script_new(); edit_script_sub (&a->script, p1->sepAfter - (p1->sepBefore+1)); a->beg1 = p1->sepBefore + 2; a->beg2 = p2->sepBefore + 2; a->end1 = p1->sepAfter; a->end2 = p2->sepAfter; a->seq1 = seq1->v; a->seq2 = seq2->v; if (s < scoreThresh.s) a->s = scoreThresh.s; // so it won't be discarded else a->s = s; a->next = NULL; a->isTrivial = true; debugSnoopEditScripts_3; debugSnoopEditScripts_4; } } // process each batch of anchors head = last = NULL; // (initialize the list of qualifying gapped alignments) for (batIx=0 ; batIxlen ; batIx++) { startSegIx = segBatches->batch[batIx].start; endSegIx = segBatches->batch[batIx].end; debugSnoopBatches_1; if (doHspsInBatches) // empty the list of bounding orderBegInc = orderEndDec = NULL; // .. alignments, for this batch // convert each anchor in this batch to a gapped extension, processing // the anchors from high score to low (they've previously been sorted // into that order) for (i=startSegIx ; ileftAlign1; io.rightAlign = mp->rightAlign1; io.leftSeg = mp->leftSeg1; io.rightSeg = mp->rightSeg1; // find the closest alignments ending before or starting after this // anchor debugSnoopAnchorToGapped_2; io.anchor1 = mp->pos1; io.anchor2 = mp->pos2; io.hspId = mp->hspId; dbg_timing_gapped_extend_sub (debugClockAboveBelow); get_above_below (&io, orderBegInc, orderEndDec); dbg_timing_gapped_extend_add (debugClockAboveBelow); // if either sequence is partitioned, figure out the limits of the // partition containing this anchor if (segBatches->batch[batIx].part != NULL) { p1 = segBatches->batch[batIx].part; goto set_limits1; } else if (sp1->p != NULL) { p1 = lookup_partition (seq1, io.anchor1); set_limits1: io.low1 = p1->sepBefore + 1; io.high1 = p1->sepAfter; } if (sp2->p != NULL) { p2 = lookup_partition (seq2, io.anchor2); io.low2 = p2->sepBefore + 1; io.high2 = p2->sepAfter; } // if we have a chore, further restrict the limits to the chore if (seq2->choresFile != NULL) { interval tInt = seq2->chore.targetInterval; interval qInt = seq2->chore.queryInterval; if (tInt.s > io.low1) io.low1 = tInt.s; if (tInt.e < io.high1) io.high1 = tInt.e; if (qInt.s > io.low2) io.low2 = qInt.s; if (qInt.e < io.high2) io.high2 = qInt.e; } // extend this anchor into a gapped alignment, in both directions dbg_timing_gapped_extend_sub (debugClockYdropAlign); ydrop_align (&io); dbg_timing_gapped_extend_add (debugClockYdropAlign); debugSnoopAnchorToGapped_3; mp->align = format_alignment (&io, mp); mp->pos1 = io.start1; mp->pos2 = io.start2; mp->end1 = io.stop1; mp->end2 = io.stop2; if (mp->firstSeg == NULL) // (the gapped alignment is empty, { // .. so skip it) debugSnoopEditScripts_5; continue; } debugSnoopEditScripts_4; // record the alignment's tail and detach the circular pointer mp->lastSeg = mp->firstSeg->prevSeg; mp->firstSeg->prevSeg = mp->lastSeg->nextSeg = NULL; // if this alignment doesn't meet the score threshold, discard it // now; otherwise, save it and use it to bound subsequent gapped // extensions; note that in blastz, the alignment was *always* // used as a bound (regardless of whether it met score threshold), // and we provide that functionality if allBounds is true (in which // case the low-scoring alignments are discarded later) if ((!allBounds) && (mp->align->s < scoreThresh.s)) { debugGappedExtendDbgShowIdentity_1; debugSnoopBlocks_2; free_align_list (mp->align); for (bp=mp->firstSeg ; bp!=NULL ; bp=bq) { bq = bp->nextSeg; free_if_valid ("gapped_extend seg", bp); } continue; } // record the horizontal bounding alignments/segments of this // anchor's gapped extension, and insert it into the vertically // ordered alignment lists dbg_timing_gapped_extend_sub (debugClockLeftRight); align_left_right (orderBegInc, mp); insert_align (mp, &orderBegInc, &orderEndDec); dbg_timing_gapped_extend_add (debugClockLeftRight); // if we have a limit on the number of paired bases we'll // accept, check for that now; if we've exceeded the limit, // we just stop processing HSPs if (maxPairedBases > 0) { newPairedBases = count_paired_bases (mp); pairedBases += newPairedBases; //fprintf (stderr, "paired bases: %s of %s\n", // commatize(newPairedBases), // commatize(pairedBases)); if (pairedBases > maxPairedBases) { if (overlyPairedWarn) warn_for_paired_bases_limit (seq2, maxPairedBases, overlyPairedKeep); if (!overlyPairedKeep) goto discard_alignments; break; // exit the HSP loop } } debugGappedExtendVerbosity_2; } // link the high scoring alignments together into a list; discard the // trivial self-alignment if it isn't wanted, and discard any alignments // that don't meet the score threshold // note that unless allBounds is true (blastz compatibility), there will // be no low-scoring alignments to discard here // also note that if we have more than one batch, the self-alignment // won't be here, so we don't have to worry that we'll discard it too // early (but since we like to worry, we do check for that case) debugSnoopBlocks_3 for (mp=orderBegInc; mp!=NULL ; mp=mpNext) { mpNext = mp->next; if (mp->align->s < scoreThresh.s) { debugGappedExtendDbgShowIdentity_2; goto free_up_all; } if ((inhibitTrivial) && (mp->align->isTrivial)) goto free_up_all; else if (delayedCheckForTrivial) { aliseg* seg = mp->firstSeg; char* name1, *name2; if (mp->lastSeg != seg) goto not_trivial; if (seg->type != diagSeg) goto not_trivial; if (doHspsInBatches) suicidef ("internal error, attempt to discard trivial self-alignment with batches in play"); if (sp1->p == NULL) { name1 = seq1->header; len1 = seq1->trueLen; if ((name1 != NULL) && (name1[0] == '>')) name1 = skip_whitespace(name1+1); } else if (segBatches->batch[batIx].part != NULL) { p1 = segBatches->batch[batIx].part; goto set_name1; } else { p1 = lookup_partition (seq1, mp->pos1); set_name1: name1 = &sp1->pool[p1->header]; len1 = p1->trueLen; } if (sp2->p == NULL) { name2 = seq2->header; len2 = seq2->trueLen; if ((name2 != NULL) && (name2[0] == '>')) name2 = skip_whitespace(name2+1); } else { p2 = lookup_partition (seq2, mp->pos2); name2 = &sp2->pool[p2->header]; len2 = p2->trueLen; } debugTriviality_1; if (len1 != len2) { debugTriviality_2; goto not_trivial; } if (mp->end1+1 - mp->pos1 != len1) { debugTriviality_3; goto not_trivial; } if (strcmp (name1, name2) != 0) { debugTriviality_4; goto not_trivial; } for (pos1=mp->pos1,pos2=mp->pos2 ; pos1<=mp->end1 ; pos1++,pos2++) { if (seq1->v[pos1] != seq2->v[pos2]) { debugTriviality_5; goto not_trivial; } } debugTriviality_6; goto free_up_all; } not_trivial: if (head == NULL) head = last = mp->align; else last = last->next = mp->align; debugGappedExtendDbgShowIdentity_3; debugSnoopBlocks_4; goto free_up_segments; // discard an alignment block free_up_all: debugSnoopBlocks_5; free_align_list (mp->align); // (fall thru to free_up_segments) // discard an alignment block's segments free_up_segments: for (bp=mp->firstSeg ; bp!=NULL ; bp=bq) { bq = bp->nextSeg; free_if_valid ("gapped_extend seg", bp); } } } debugSnoopAnchorToGapped_4; free_if_valid ("gapped_extend msp[]", msp); return head; ////////// // failure exit ////////// // discard all alignments we found discard_alignments: debugSnoopBlocks_3b for (mp=orderBegInc; mp!=NULL ; mp=mpNext) { mpNext = mp->next; // discard an alignment block and its segments debugSnoopBlocks_5; free_align_list (mp->align); for (bp=mp->firstSeg ; bp!=NULL ; bp=bq) { bq = bp->nextSeg; free_if_valid ("gapped_extend seg", bp); } } debugSnoopAnchorToGapped_5; free_if_valid ("gapped_extend msp[]", msp); return NULL; } //---------- // // batched_segments-- // Separate a list of segments in batches relative to a partitioned sequence. // Conceptually we produce one batch per partition. But if a partition // contains no segment, it is just left out of the batch table. // //---------- // // Arguments: // segtable* anchors: The anchors. Note that these well be reordered // .. by this routine, so that every batch is // .. sorted by descreasing score. // seqpartition* sp: The partitioning to separate by. This should // .. relate to sequence 1. A special value of // .. NULL indicates that partitioning is not // .. wanted, so we should treat the whole segment // .. table as a single batch. // // Returns: // A pointer to the batch table. This is allocated data, which the caller // must eventually dispose of by calling free_segment_batches(). // //---------- static sbtable* _segBatches = NULL; static sbtable* batched_segments (segtable* anchors, seqpartition* sp) { u32 entriesNeeded; size_t bytesNeeded; partition* part; segment* anc; unspos pEnd; u32 batIx, partIx, ancIx; u32 startIx, endIx; // allocate the batch table, or resize it if it's not big enough // nota bene: (as of this writing), we don't expect realloc to ever occur, // because the sequence that we're partitioning with is the same // one through the entire program; that may change in the future if (sp == NULL) entriesNeeded = 1; else entriesNeeded = sp->len; bytesNeeded = sbtable_bytes (entriesNeeded); if (_segBatches == NULL) { _segBatches = (sbtable*) malloc_or_die ("batched_segments", bytesNeeded); _segBatches->size = entriesNeeded; } else if (entriesNeeded > _segBatches->size) { _segBatches = (sbtable*) realloc_or_die ("batched_segments", _segBatches, bytesNeeded); _segBatches->size = entriesNeeded; } // if we don't have a partitioned sequence, create one batch convering the // entire segment list if (sp == NULL) { _segBatches->batch[0].part = NULL; _segBatches->batch[0].start = 0; _segBatches->batch[0].end = anchors->len; _segBatches->len = 1; sort_segments (anchors, qSegmentsByDecreasingScore); } // otherwise, create batches, dividing the segments by partition; note // that we assume each segment is contained within a single partition else { sort_segments (anchors, qSegmentsByPos1); batIx = 0; ancIx = 0; anc = &anchors->seg[ancIx]; for (partIx=0 ; partIxlen ; partIx++) { if (ancIx >= anchors->len) break; part = &sp->p[partIx]; pEnd = part->sepAfter; if (pEnd < anc->pos1 + anc->length) continue; _segBatches->batch[batIx].part = part; _segBatches->batch[batIx].start = startIx = ancIx++; anc++; while ((ancIx < anchors->len) && (pEnd >= anc->pos1 + anc->length)) { ancIx++; anc++; } _segBatches->batch[batIx].end = endIx = ancIx; batIx++; sort_some_segments (anchors, startIx, endIx, qSegmentsByDecreasingScore); } _segBatches->len = batIx; } return _segBatches; } //---------- // // free_segment_batches-- // Dispose of memory alloceted by batched_segments(). // //---------- // // Arguments: // (none) // // Returns: // (nothing) // //---------- void free_segment_batches (void) { free_if_valid ("free_segment_batches", _segBatches); _segBatches = NULL; } //---------- // // init_from_anchors-- // Initialize a gapped alignments set from a set of anchors. // //---------- // // Arguments: // segtable* anchors: The anchors, which we will expand into gapped // .. alignment records, intitally as single // .. points. // u32 numExtraSlots: The number of extra slots to add (see note 1 // .. below) // // Returns: // A pointer to newly allocated data (see description below); failures result // in program fatality. The caller must eventually dispose of the table, with // a call to free(). // // Data Description: // The allocated data is an array of pointers to galign blocks. It is // allocated as a single block from the heap (which is why it can be freed wih // a simple call to free). The first N locations form the array, and the data // following that contains the galign blocks. // //---------- // // Notes: // // (1) We allocate extra slots in the gapped alignment table, to allow a // trivial self-alignments to be added later. If the target sequence is // *not* partitioned, there will be at most one trivial alignment and we // add only one extra slot. The same is true if the query is *not* // partitioned. If both target and query are partitioned, there can be // a trivial alignment for each partition in whichever has fewer // partitions, so we add that many extra slots. We push this determination // back to the caller, who must provide the value for numExtraSlots. // // (2) We only copy pos1/pos2 from the anchors. We ignore the length and // don't set end1/end2. The expectation is that the anchors are single // points without length-- upstream code should have reduced them to // single points, and initial downstream code only uses pos1/pos2. // //---------- static galign** init_from_anchors (segtable* anchors, u32 numExtraSlots) { int numAnchors, numSlots, ix; u32 bytesNeeded, bytesPointers, bytesBlocks; galign** snakes, *snake; segment* anchor; numAnchors = anchors->len; numSlots = numAnchors + numExtraSlots; // allocate; note that we allocate extra slots as per note (1) above // $$$ Why is zalloc used here instead of malloc? If the table is large, // $$$ .. zalloc will spec significant time zeroing it, but why is that // $$$ .. necessary? This sets the anchor's length to zero, but would it // $$$ .. be more efficient to set that explicitely, so that we don't waste // $$$ .. time clearing out the extra slots which are often not used? bytesPointers = round_up_16 (numSlots * sizeof(galign*)); if ((u32max-16) / sizeof(galign) < (u32) numSlots) goto overflow1; bytesBlocks = round_up_16 (numSlots * sizeof(galign)); bytesNeeded = bytesPointers; bytesNeeded += bytesBlocks; if (bytesNeeded < bytesBlocks) goto overflow2; snakes = (galign**) zalloc_or_die ("init_from_anchors", bytesNeeded); gapped_extend_count_stat (zallocCallsA); gapped_extend_add_stat (zallocTotalA, bytesNeeded); // hook up pointers; this links the array of N pointers to the actual // blocks snakes[0] = (galign*) (((char*) snakes) + bytesPointers); for (ix=1 ; ixseg[ix]; snake->pos1 = anchor->pos1; snake->pos2 = anchor->pos2; snake->hspId = anchor->hspId; } return snakes; // failure exits #define suggestions " consider raising scoring threshold (--hspthresh or --exact)" \ " or breaking your target sequence into smaller pieces" overflow1: suicidef ("in init_from_anchors(), structure size would exceed 2^32" " (%u + %u*%u)\n" suggestions, bytesPointers, sizeof(galign), numSlots); return NULL; // (doesn't get here) overflow2: suicidef ("in init_from_anchors(), structure size exceeds 2^32" " (%u + %u)\n" suggestions, bytesPointers, bytesBlocks); return NULL; // (doesn't get here) } //---------- // // identical_sequences-- // Determine if two sequences are identical for practical purposes. // //---------- // // Arguments: // seq* seq1: One sequence // seq* seq2: The other sequence. // scoreset* scoring: The scoring scheme to use. This may be NULL if the // .. caller isn't interested in the score. // score* s: Place to return the score. This may be NULL. // // Returns: // true if the sequences are identical, but allowing differences in upper/lower // case; false otherwise. // //---------- // // Notes: // // (1) (this note also applies to score_identical_partition and // score_identical_partition_of) // When computing the score below we have to consider the possibility // of overflow. For example, if matches are worth 100, then a sequence // of length 21.5M will exceed 2^31 and overflow (assuming scores are // 32-bit ints). We prevent this by checking whether adding a positve // score will overflow. However, for negative match scores (e.g. N vs N), // we simply subtract. Thus we are assuming that that we won't have so // many negative scores that we underflow. // //---------- static int identical_sequences (seq* seq1, seq* seq2, scoreset* scoring, score* _s) { seqpartition* sp1 = &seq1->partition; seqpartition* sp2 = &seq2->partition; u8* a, *b; u8 aNuc, bNuc; score s, sub; unspos ix; if ((sp1->p != NULL) || (sp2->p != NULL)) // one or the return false; // .. other is partitioned if ((seq1->fileType == seq_type_qdna) != (seq2->fileType == seq_type_qdna)) return false; if (seq1->len != seq2->len) return false; if (seq1->revCompFlags != seq2->revCompFlags) return false; a = seq1->v; b = seq2->v; s = 0; for (ix=0 ; ixlen ; ix++) { aNuc = (u8) dna_toupper(a[ix]); bNuc = (u8) dna_toupper(b[ix]); if (aNuc != bNuc) return false; if (scoring == NULL) continue; sub = scoring->sub[aNuc][bNuc]; // (see note above about overflow) if (s == bestPossibleScore) ; else if ((sub <= 0) || (s < bestPossibleScore - sub)) s += sub; else s = bestPossibleScore; } if (_s != NULL) *(_s) = s; return true; } //---------- // // identical_partitioned_sequences-- // Determine if two partitioned sequences are identical for practical purposes. // //---------- // // Arguments: // seq* seq1: One sequence. // seq* seq2: The other sequence. // // Returns: // true if the sequences are identical across their partitions, but allowing // differences in upper/lower case; false otherwise. // //---------- static int identical_partitioned_sequences (seq* seq1, seq* seq2) { seqpartition* sp1 = &seq1->partition; seqpartition* sp2 = &seq2->partition; partition* p1, *p2; u32 partIx; u8* a, *b; u8 aNuc, bNuc; unspos len1, len2, ix; if ((sp1->p == NULL) || (sp2->p == NULL)) // one or the other return false; // .. is *not* partitioned if ((seq1->fileType == seq_type_qdna) != (seq2->fileType == seq_type_qdna)) return false; if (sp1->len != sp2->len) return false; if (seq1->revCompFlags != seq2->revCompFlags) return false; for (partIx=0 ; partIxlen ; partIx++) { p1 = &sp1->p[partIx]; p2 = &sp2->p[partIx]; len1 = p1->sepAfter - (p1->sepBefore+1); len2 = p2->sepAfter - (p2->sepBefore+1); if (len1 != len2) return false; a = seq1->v + p1->sepBefore+1; b = seq2->v + p2->sepBefore+1; for (ix=0 ; ixpartition; seqpartition* sp2 = &seq2->partition; partition* p1; u32 partIx; u8* a, *b; u8 aNuc, bNuc; unspos len1, len2, ix; int isMatch; if (sp1->p == NULL) return -1; // seq1 is *not* partitioned if (sp2->p != NULL) return -1; // seq2 is partitioned if ((seq1->fileType == seq_type_qdna) != (seq2->fileType == seq_type_qdna)) return -1; if (seq1->revCompFlags != seq2->revCompFlags) return -1; // if this is a previously unseen partitioning, prescan to find the min and // max lengths #ifdef cache_partition_lengths if (sp1 != cachedSp1) { cachedSp1 = sp1; partIx = 0; p1 = &sp1->p[partIx]; len1 = p1->sepAfter - (p1->sepBefore+1); cachedMinLen = cachedMaxLen = len1; for (partIx=1 ; partIxlen ; partIx++) { p1 = &sp1->p[partIx]; len1 = p1->sepAfter - (p1->sepBefore+1); if (len1 < cachedMinLen) cachedMinLen = len1; else if (len1 > cachedMaxLen) cachedMaxLen = len1; } //fprintf (stderr, "partition lengths: " unsposFmt ".." unsposFmt "\n", // cachedMinLen, cachedMaxLen); } #endif // cache_partition_lengths // if the second sequence is shorter or longer than all partitions, there // can be no match len2 = seq2->len; #ifdef cache_partition_lengths if (sp1 == cachedSp1) { if ((len2 < cachedMinLen) || (len2 > cachedMaxLen)) { // fprintf (stderr, "length out of range: " unsposFmt "\n", len2); return -1; } } #endif // cache_partition_lengths // otherwise, scan all partitions for a match for (partIx=0 ; partIxlen ; partIx++) { p1 = &sp1->p[partIx]; len1 = p1->sepAfter - (p1->sepBefore+1); if (len1 != len2) continue; a = seq1->v + p1->sepBefore+1; b = seq2->v; isMatch = true; for (ix=0 ; ixsepAfter - (p1->sepBefore+1); a = seq1->v + p1->sepBefore+1; b = seq2->v + p2->sepBefore+1; s = 0; for (ix=0 ; ixsub[aNuc][bNuc]; // (see note above about overflow) if (s == bestPossibleScore) ; else if ((sub <= 0) || (s < bestPossibleScore - sub)) s += sub; else s = bestPossibleScore; } return s; } //---------- // // score_identical_partition_of-- // Compute the score of a partition versus a sequence, which are assumed to // have indentical length. // //---------- // // Arguments: // seq* seq1: The sequence containing p1. // seq* seq2: The other sequence (not partitioned). // partition* p1: The partition. // scoreset* scoring: The scoring scheme to use. // // Returns: // the score. // //---------- // // Notes: // // (1) We assume, WITHOUT CHECKING, that the partition and sequence have the // same length. // // (2) (see note 1 in identical_sequences, regarding overflow) // //---------- static score score_identical_partition_of (seq* seq1, seq* seq2, partition* p1, scoreset* scoring) { u8* a, *b; u8 aNuc, bNuc; score s, sub; unspos len, ix; len = p1->sepAfter - (p1->sepBefore+1); a = seq1->v + p1->sepBefore+1; b = seq2->v; s = 0; for (ix=0 ; ixsub[aNuc][bNuc]; // (see note above about overflow) if (s == bestPossibleScore) ; else if ((sub <= 0) || (s < bestPossibleScore - sub)) s += sub; else s = bestPossibleScore; } return s; } //---------- // // new_traceback-- // Allocate some traceback data. // //---------- // // Arguments: // u32 size: The maximum number of bytes to allocate. // // Returns: // A pointer to a newly allocated traceback data; failures result in program // fatality. The caller must eventually dispose of the table, with a call to // free_traceback(). // //---------- tback* new_traceback (u32 size) { tback* tb; int cells; // sanity check if (size < sizeof(tback)) suicidef ("in new_traceback(), size can't be %u", size); // allocate cells = 1 + (size - sizeof(tback)) / sizeof(tb->space[0]); tb = (tback*) malloc_or_die ("new_traceback", size); // initialize tb->size = cells; return tb; } //---------- // // free_traceback-- // Deallocate traceback data. // //---------- // // Arguments: // tback* tb: The traceback data to dispose of. // // Returns: // (nothing) // //---------- void free_traceback (tback* tb) { free_if_valid ("free_traceback", tb); } //---------- // // ydrop_align-- // The Y-drop variant of dynamic programming is applied to create a gapped // alignment by extending in both directions from an "anchor point". // //---------- // // Arguments: // alignio* io: The collection of input arguments, and a place to store // .. results. The input arguments describe the alignment // .. problem to be solved. See the definition of the type // .. alignio for more details. // // Returns: // nothing; actual result values are in io: s, start1, start2, stop1, stop2 // and script. // //---------- static void lop_initial_indels (alignio* io); static void lop_final_indels (alignio* io); #ifndef snoopSubprobs #define snoopSubprobsA_V ; #define snoopSubprobsA_1 ; #define snoopSubprobsA_2 ; #define snoopSubprobsA_3 ; #endif // not snoopSubprobs #ifdef snoopSubprobs #define snoopSubprobsA_V \ partition* part; \ char* name1, *name2 #define snoopSubprobsA_1 \ if (snoopSubprobsSeq1->partition.p == NULL) \ { \ name1 = (snoopSubprobsSeq1->useFullNames)? snoopSubprobsSeq1->header \ : snoopSubprobsSeq1->shortHeader; \ if ((name1 == NULL) || (name1[0] == 0)) name1 = "seq1"; \ } \ else \ { \ part = lookup_partition (snoopSubprobsSeq1, anchor1); \ name1 = &snoopSubprobsSeq1->partition.pool[part->header]; \ } \ \ if (snoopSubprobsSeq2->partition.p == NULL) \ { \ name2 = (snoopSubprobsSeq2->useFullNames)? snoopSubprobsSeq2->header \ : snoopSubprobsSeq2->shortHeader; \ if ((name2 == NULL) || (name2[0] == 0)) name2 = "seq2"; \ } \ else \ { \ part = lookup_partition (snoopSubprobsSeq2, anchor2); \ name2 = &snoopSubprobsSeq2->partition.pool[part->header]; \ } \ \ fprintf (stderr, "align:\tbck\t%s\t" unsposFmt "\t" unsposFmt \ "\t%s\t" unsposFmt "\t" unsposFmt, \ name1, io->len1-anchor1-2, (anchor1+1)-io->low1, \ name2, io->len2-anchor2-2, (anchor2+1)-io->low2); #define snoopSubprobsA_2 \ if (snoopSubprobsSeq1->partition.p == NULL) \ { \ name1 = (snoopSubprobsSeq1->useFullNames)? snoopSubprobsSeq1->header \ : snoopSubprobsSeq1->shortHeader; \ if ((name1 == NULL) || (name1[0] == 0)) name1 = "seq1"; \ } \ else \ { \ part = lookup_partition (snoopSubprobsSeq1, anchor1); \ name1 = &snoopSubprobsSeq1->partition.pool[part->header]; \ } \ \ if (snoopSubprobsSeq2->partition.p == NULL) \ { \ name2 = (snoopSubprobsSeq2->useFullNames)? snoopSubprobsSeq2->header \ : snoopSubprobsSeq2->shortHeader; \ if ((name2 == NULL) || (name2[0] == 0)) name2 = "seq2"; \ } \ else \ { \ part = lookup_partition (snoopSubprobsSeq2, anchor2); \ name2 = &snoopSubprobsSeq2->partition.pool[part->header]; \ } \ \ fprintf (stderr, "\nalign:\tfwd\t%s\t" unsposFmt "\t" unsposFmt \ "\t%s\t" unsposFmt "\t" unsposFmt, \ name1, anchor1, io->high1-(anchor1+1), \ name2, anchor2, io->high2-(anchor2+1)); #define snoopSubprobsA_3 \ fprintf (stderr, "\n"); #endif // snoopSubprobs static void ydrop_align (alignio* io) { unspos anchor1, anchor2; unspos end1, end2; score scoreLeft, scoreRight; editscript* script, *scriptRight; u32 op; #ifdef snoopAlgorithm int snoop = true; #endif // snoopAlgorithm snoopSubprobsA_V; if (io == NULL) suicide ("ydrop_align() called with NULL pointer."); #ifdef snoopAlignioInput dump_alignio_input (stderr, io); #endif // snoopAlignioInput gapped_extend_count_stat (numAnchorsExtended); anchor1 = io->anchor1; anchor2 = io->anchor2; #ifdef snoopAlgorithm #if ((defined debugPosA1) && (defined debugPosB1)) snoop = ((anchor1 == debugPosA1) && (anchor2 == debugPosA2)) || ((anchor1 == debugPosB1) && (anchor2 == debugPosB2)); #elif (defined debugPosA1) snoop = (anchor1 == debugPosA1) && (anchor2 == debugPosA2); #elif (defined debugPosB1) snoop = (anchor1 == debugPosB1) && (anchor2 == debugPosB2); #endif // debugPosA1,debugPosB1 #endif // snoopAlgorithm #ifdef snoopAlgorithm if (snoop) { fprintf (stderr, "gapE = " scoreFmtSimple "\n", io->scoring->gapExtend); fprintf (stderr, "gapOE = " scoreFmtSimple "\n", io->scoring->gapOpen + io->scoring->gapExtend); fprintf (stderr, "yDrop = " scoreFmtSimple "\n", io->yDrop); dump_score_set (stderr, io->scoring, (u8*)"ACGTacgtNn", (u8*)"ACGTacgtNn"); } #endif // snoopAlgorithm snoopSubprobsA_1; script = edit_script_new(); dbg_timing_gapped_extend_sub (debugClockYdropOneSidedAlign); scoreLeft = ydrop_one_sided_align (io, /*reverse*/ true, io->rev1 + io->len1 - anchor1 - 2, io->rev2 + io->len2 - anchor2 - 2, (anchor1+1) - io->low1, (anchor2+1) - io->low2, io->trimToPeak, &script, &end1, &end2); dbg_timing_gapped_extend_add (debugClockYdropOneSidedAlign); io->start1 = anchor1 + 1 - end1; io->start2 = anchor2 + 1 - end2; snoopSubprobsA_2; scriptRight = edit_script_new(); dbg_timing_gapped_extend_sub (debugClockYdropOneSidedAlign); scoreRight = ydrop_one_sided_align (io, /*reverse*/ false, io->seq1 + anchor1, io->seq2 + anchor2, io->high1 - (anchor1+1), io->high2 - (anchor2+1), io->trimToPeak, &scriptRight, &end1, &end2); dbg_timing_gapped_extend_add (debugClockYdropOneSidedAlign); io->stop1 = anchor1 + end1; io->stop2 = anchor2 + end2; snoopSubprobsA_3; #ifdef snoopAlgorithm if (snoop) { fprintf (stderr, "==== left edit script:\n"); dump_edit_script (stderr, script); fprintf (stderr, "==== right edit script:\n"); dump_edit_script (stderr, scriptRight); } #endif // snoopAlgorithm edit_script_reverse (scriptRight); edit_script_append (&script, scriptRight); free_if_valid ("ydrop_align, edit script", scriptRight); io->s = scoreRight + scoreLeft; io->script = script; #ifdef snoopAlgorithm if (snoop) { fprintf (stderr, "==== combined edit script:\n"); dump_edit_script (stderr, script); fprintf (stderr, "\n"); } #endif // snoopAlgorithm // in rare cases, the DP algorithm can produce an alignment that starts or // ends with a gap; remove those here if (io->script->len != 0) { op = edit_op_operation (io->script->op[0]); if (op != editopSub) lop_initial_indels (io); // nota bene: lop_initial_indels will not set len to zero op = edit_op_operation (io->script->op[io->script->len-1]); if (op != editopSub) lop_final_indels (io); } #if ((defined(snoopAlignioOutput)) || (defined(snoopEditScripts))) dump_alignio_output (stderr, io); #endif // snoopAlignioOutput OR snoopEditScripts } // lop_initial_indels-- static void lop_initial_indels (alignio* io) { unspos pos1, pos2; u32 opIx, op, rpt, numIndelOps; // scan past all indel ops at the start pos1 = io->start1; pos2 = io->start2; for (opIx=0 ; opIxscript->len ; opIx++) { op = edit_op_operation (io->script->op[opIx]); rpt = edit_op_repeat (io->script->op[opIx]); if (op == editopSub) break; else if (op == editopIns) pos2 += rpt; else if (op == editopDel) pos1 += rpt; } numIndelOps = opIx; // special case if the alignment was nothing but indels; assign it a // phenominally bad score so that something downstream from us will discard // it if (numIndelOps == io->script->len) { io->s = worstPossibleScore; return; } // modify the starting position, copy the remainder of the edit script // forward, and rescore the alignment io->start1 = pos1; io->start2 = pos2; io->script->len -= numIndelOps; for (opIx=0 ; opIxscript->len ; opIx++) io->script->op[opIx] = io->script->op[opIx+numIndelOps]; io->s = score_alignment (io->scoring, io->seq1, io->start1, io->seq2, io->start2, io->script); } // lop_final_indels-- static void lop_final_indels (alignio* io) { unspos pos1, pos2; u32 opIx, op, rpt; // scan past all indel ops at the end pos1 = io->stop1; pos2 = io->stop2; for (opIx=io->script->len ; opIx>0 ; ) { opIx--; op = edit_op_operation (io->script->op[opIx]); rpt = edit_op_repeat (io->script->op[opIx]); if (op == editopSub) { opIx++; break; } else if (op == editopIns) pos2 -= rpt; else if (op == editopDel) pos1 -= rpt; } // special case if the alignment was nothing but indels; assign it a // phenominally bad score so that something downstream from us will discard // it if (opIx == 0) { io->s = worstPossibleScore; return; } // modify the ending position, truncate the edit script, and rescore the // alignment io->stop1 = pos1; io->stop2 = pos2; io->script->len = opIx; io->s = score_alignment (io->scoring, io->seq1, io->start1, io->seq2, io->start2, io->script); } //---------- // // ydrop_one_sided_align-- // Find an gapped alignment in one direction from a starting point. // //---------- // // Arguments: // alignio* io: Other input arguments for the alignment problem // .. being solved. Of interest to this routine are // .. the scoring parameters (scoring and yDrop) and // .. the rightSeg/leftSeg stuff, which is used to // .. track sweep row bounds caused by previous // .. alignments. // int reversed: true => search to the lower-left // false => search to the upper-right // u8* A: Query sequence (conceptually vertical) // .. (indexed by 1..M, see note below) // u8* B: Subject sequence (conceptually horizontal) // .. (indexed by 1..N, see note below) // unspos M: Length of query sequence. // unspos N: Length of subject sequence. // int trimToPeak: true => always trim end of alignment to score peak. // false => don't trim if extension runs into end of // .. sequence (see note 12 below). // editscript** script: Place to return the resulting edit script. Note // .. that this is in traceback order, with the // .. alignment's starting point corresponding to the // .. final operation in the edt script. // unspos* end1: Place to return the alignment end, in A. // unspos* end2: Place to return the alignment end, in B. // // Returns: // The score of the alignment. // //---------- // // Notes: // - A and B are indexed starting from 1, so the pointers must point to one // character before the start of the string. That character is ignored. // // - Whether reversed is true or false, the DP matrix is conceptually scanned // from bottom-up, and from left-to-right (see diagram below). What changes // is the equivalence between DP row,col and sequence row,col. // // ------ right -----> // +----+------------------+ // row M |A[M]| | ^ // ... |... | | | // row 1 |A[1]| DP Matrix | up // row 0 | | | | // | +------------------+ | // | B[1] ... B[N]| // +-----------------------+ // col col ... col // 0 1 N // //---------- // // Implementation Notes: // // (1) We manage the external bounds of the current sweep row by keeping track // of the alignment and the specific segments that constrain the feasible // region on the left (leftAlign,leftSeg) and right (rightAlign,rightSeg), // if they exist. These are initially set to the constraints of the anchor // point, and information is maintained to allow efficient updating as the // sweep row advances. // // For a given row, the values L and R bracket the legal DP columns to // compute. Neither column L nor R is a valid column in that row, being // part of a previous alignment or off the left or right edge of the // matrix. Similarly, LY and RY bracket the intersection of the y-drop // feasible region and the left/right bounds. // // (2) The new alignment should not be allowed to intersect any previously // found alignments, so we have to keep track of "active segments"-- gap- // free segments of earlier alignments that intersect the sweep row within // the feasible region. These are kept in a list that supports insertion // and deletion when the sweep row reaches either end of an earlier // alignment. // // As we advance the sweep row, update_active_segs identifies all columns // (within the feasible range) that intersect an existing alignment. It // marks the corresponding DP cells as 'masked' (the mask indicator value // is the current row number, so that we don't have to erase them for the // next sweep row). When the sweep encounters a masked cell, it refuses // to step to it, and sets all three scores to -infinity so that no step // will be taken from it. // // (3) Memory for traceback cells is provided by the caller (in io->tb). That // memory is carved into pieces matching the computed rows/slices of the DP // matrix. This routine allocates an array of indexes to the start of // each row (with the index offset to the left so that columns index the // correct cells). // // The memory for this index array is allocated here, and lives across // calls as a static variable. We allocate it in multiples of 512K in the // hope that this improves performance on some platforms (specifically, to // avoid expensive mmap/munmap activity in linux with gnu malloc); it is // not clear whether this helps. // // We use an array of indexes, rather than an array of pointers, to reduce // the memory requirement on 64-bit machines. Pointers on such machines // would require 8 bytes per entry, whereas indexes require only 4 bytes. // // The BLASTZ version of this function allocated a (potentially) very large // pointer array. It allocated at least M+1 entries, one for each possible // row. For long sequences such as human chromosome 1 (~ 250 MB), this // required 1GB of memory if the anchor was close to either end of the // sequence. Typical alignments are very unlikely to use that many rows, // especially if B is much shorter than A. So in this implementation, the // pointer array is initially allocated as 1% of the the sequence length, // but not more than the length of B, and rounded up to a multiple of 512K // bytes. It is "sort of" expanded at a rate of 6% whenever the DP scan // requires more rows (in practice, it just allocates the next chunk of // 512K bytes). We don't want to have to reallocate very often, because // realloc may have to copy the entire contents to the new buffer (even if // we don't need it to, as when we are starting a new alignment; too bad // malloc/realloc doesn't provide a routine that lets us say "I need more, // but I don't need you to preserve the contents). // // (4) We only use DP memory corresponding to one row/slice of the matrix. At // any time, we need information about cells in one row and the row above // it, but we only need both rows at the scan point. We use three local // variables to represent the extra cell there. // // One pointer (dp) scans the cells for the lower row (the row being read), // another (dq) scans the upper row. These pointers are the same when // the two rows have the same left bound. When the upper row has a // different bound (which is necessarily further to the right), dq will // lag behind dp. // // (5) The inner loop scans the current row from left to the right, starting // with column LY and ending before column RY. Note that column LY is an // invalid column (it is column zero, or it is on a segment, or it is // beyond ydrop); this corresponds to the traditional full DP algorithm // computing column zero even though column one is the first character of // B. // // The following are true for each pass through the inner loop (see note // (4) for information about dp and dq): // // We enter the loop each time with // dp points to cell to read DP[*][col] // dq points to cell to write DP[*][col] // C[row-1][col] in dp->CC // C[row ][col] in c (proposed, not final) // D[row ][col] in dp->DD // I[row ][col] in i // A[row ][* ] reflected in sub[*] // B[col+1] in *b // And during the loop, we will write // C[row ][col ] to dq->CC, and check/update bestScore // C[row ][col+1] to c (proposed) // I[row ][col+1] to i // D[row+1][col ] to dq->DD // link into C[row][col] to *tbp // // It is important to realize that the stored values in one cell are for // the C node in one row and the D node for the row above it. // // (6) At the beginning of each row, we grab a pointer (sub) to the row of the // substitution scores matrix for the character A[row]. At the end of each // pass through the inner loop, we look up the score for the character // B[col+1]. Thus we are computing the value of a substitution into // node C[row,col+1]. // // (7) The trace-back information, i.e., to trace backwards along an optimal // path, is packed into one byte per dynamic-programming grid point. At // each grid point there are three nodes: C which is entered by one // diagonal edge, D which is entered by two vertical edges (from a D node // and a C node), and I which is entered by two horizontal edges (from an // I node and a C node). See ref [1] for further details. // // The right-most two bits hold cFromC (0), cFromI (1) or cFromD (2), // telling how this cell's C node gets the maximum score from an entering // edge. // // iExtend (4) and dExtend (8) relate to edges *exiting* this cell to the I // or D node in the next column or row, respectively. iExtend is set iff // it is better to take a horizontal edge from this cell's I node (a gap // extend) than a horizontal edge from the C node (a gap open). dExtend // has a similar meaning for vertical edges. // // (8) The loop to compute the first row of the DP matrix is limited by col<=N, // but instead this ought to be col=-yDrop, will keep us from // going very far (with typical values, this condition will stop us after // around 300 columns). // // (9) When we stop computing one row at its right bound, we may still need // to prolong computation in the row to support an overhang on the row // above. We do so by allowing a series of insertions. // // The following are true for each pass through the row prolongation loop: // dq points to cell to write DP[*][col] // I[row][col] in i // // (10) We use the loop body below left even though the one below right seems // more mathematically correct. // // dq->CC = c; dq->CC = c; // (dq++)->DD = c - gapOE; c -= gapE; // c -= gapE; (dq++)->DD = c; // *(tbp++) = cFromI; *(tbp++) = cFromI; // // However, using the one on the right caused problems aligning // // A=GTTTTTTTTTTTTTTTCTT // B=TTTTTTTTTTTTTTTTTTCTT // // It preferred the alignment on the left to the one on the right, because // it failed to charge a second gap open penalty. // // ---GTTTTTTTTTTTTTTTCTT --GTTTTTTTTTTTTTTTCTT // TTT-TTTTTTTTTTTTTTTCTT TTTTTTTTTTTTTTTTTTCTT // // (11) Blastz had a potential problem with the feasibility bounds LY and RY. // update_LR_bounds occasionally will create the condition RY < LY. This // means that no columns of this dp matrix row will be computed, because // the col=RY will halt the routine. However, the calculation // of tbNeeded is suspect, since it could conceivably be negative. I don't // believe it ever goes negative in practice (because I don't think RY-LY // is ever less than -2), but for safety sake we clip it at zero. // // (12) Correction for short read "mismatch shadow". The original intent of // this routine is that the two sequences would be long (e.g. chromosomes). // In that context, it made sense to trim any negatively-scoring suffix // from the end of an alignment. However, if one of the sequences is a // short read (e.g. 50 bases), we are usually interested in aligning the // entire read if we can. If we don't, we will be creating a bias against // reporting mismatches near the end of the reads. // // Suppose scores are +1 for a match and -7 for a mismatch. Then any // mismatch within 7 bases of the end will create a negatively-scoring // suffix, even if all the bases between it and the end are matches. If // this is trimmed, that mismatch is never reported as part of an // alignment. // // The argument trimToPeak is added to control this behavior. If // trimToPeak is true, we *always* trim the alignment end back to the // maximum scoring location. If trimToPeak is false, we don't trim if we // happen to encounter the end of the sequence. // // (13) The BLASTZ loop to compute the first row of dp matrix was equivalent to // this: // // dq = dynProg->p; // dq->CC = 0; // c = (dq++)->DD = -gapOE; // // for (col=1 ; (col<=N)&&(c>=-yDrop) ; col++) // { // dq->CC = c; // (dq++)->DD = c - gapOE; // c -= gapE; // } // // The (c>=-yDrop) part of this test is not correct, because it is testing // the value of the D entry in the cell, instead of the C entry. // // (14) Bounding segment calculations are initially done w.r.t. the full DP // matrix with forward orientation of the sequences. When we are doing a // reversed alignment, L and R (the left and right bounds) must be // swapped. This swap appeared to be performed wrong in BLASTZ, and has // been changed, as shown in this table: // // left right | BLASTZ | LASTZ | // bound? bound? | L' R' | L' R' | // ---------------+---------+-----------+ // no no | L R | L R | // no yes | -R R | -R+1 N+1 | // yes no | L -L | 0 -L-1 | // yes yes | -R -L | -R+1 -L-1 | // ---------------+---------+-----------+ // //---------- //=== macros for ydrop_one_sided_align === #define cFromC 0 // (c bit is no bits) see note (7) #define cFromI 1 #define cFromD 2 #define iExtend 4 #define dExtend 8 #define cidBits (cFromC | cFromI | cFromD) #define op_string(op) ((op==cFromC)? "SUB" : (op==cFromI)? "INS" : (op==cFromD)? "DEL" : "???") #define prune \ c = dp->CC + sub[*(b++)]; /* propose C[row][col+1] */ \ if (col == LY) /* (we're at first col in sweep) */ \ LY++; \ else \ { /* 'set' I[row][col+1] */ \ i = dq->DD = dq->CC = negInf; /* .. and set D[row+1][col] */ \ dq++; /* .. and set C[row+1][col] */ \ } \ dp++; \ *(tbp++) = 0; //=== stuff for link_to_string === #define include_link_to_string \ static char* link_to_string (int link) \ { \ static char str[100]; \ char* s = str; \ *(s++) = '{'; \ if ((link & cidBits) == 0) { *(s++) = 'c'; *(s++) = ','; } \ if ((link & cFromI) != 0) { *(s++) = 'i'; *(s++) = ','; } \ if ((link & cFromD) != 0) { *(s++) = 'd'; *(s++) = ','; } \ if ((link & iExtend) != 0) { *(s++) = 'i'; *(s++) = 'i'; *(s++) = ','; } \ if ((link & dExtend) != 0) { *(s++) = 'd'; *(s++) = 'd'; *(s++) = ','; } \ if (s > str+1) s--; \ *(s++) = '}'; \ *s = 0; \ return str; \ } //=== stuff for snoopAlgorithm === #ifndef snoopAlgorithm #define snoopAlgorithm_1 ; #define snoopAlgorithm_2 ; #define snoopAlgorithm_3 ; #define snoopAlgorithm_4A ; #define snoopAlgorithm_4B ; #define snoopAlgorithm_4C ; #define snoopAlgorithm_5 ; #define snoopAlgorithm_5A ; #define snoopAlgorithm_5B ; #define snoopAlgorithm_6 ; #define snoopAlgorithm_7A ; #define snoopAlgorithm_7B ; #define snoopAlgorithm_7C ; #define snoopAlgorithm_7D ; #define snoopAlgorithm_7E ; #define snoopAlgorithm_8 ; #endif // not snoopAlgorithm #ifdef snoopAlgorithm static char ccSnoop[100]; static char ddSnoop[100]; static char iiSnoop[100]; static char* relative_to_infinity (score s) { static char str[100]; if (s < negInf) sprintf (str, "-inf" scoreFmtSimple, s-negInf); else if (s == negInf) sprintf (str, "-inf"); else if (s <= negInf+2000) sprintf (str, "-inf+" scoreFmtSimple, s-negInf); else sprintf (str, scoreFmtSimple, s); return str; } include_link_to_string #define snoopAlgorithm_1 \ if (snoop) \ { \ char A25[26], B25[26]; \ \ strncpy (A25, (char*) A+1, sizeof(A25)); A25[sizeof(A25)-1] = 0; \ strncpy (B25, (char*) B+1, sizeof(B25)); B25[sizeof(B25)-1] = 0; \ \ fprintf (stderr, "\nydrop_one_sided_align(" unsposSlashFmt ") %s" \ " A=%s (" unsposFmt ")" \ " B=%s (" unsposFmt ")\n", \ io->anchor1, io->anchor2, \ (reversed)?" reversed":"forward", \ A25, M, B25, N); \ } #define snoopAlgorithm_2 \ if (snoop) \ { \ strcpy (ccSnoop, relative_to_infinity ((dq-1)->CC)); \ strcpy (ddSnoop, relative_to_infinity ((dq-1)->DD)); \ fprintf (stderr,"init edge, dynProg=%08lX dpLen=%08X\n", \ (long) dynProg, dynProg->len); \ fprintf (stderr, "\n[%3u,%3u]" \ " [%d].cc=%s [%d].DD=%s" \ " dq=%d\n", \ 0, 0, 0, ccSnoop, 0, ddSnoop, 0); \ } #define snoopAlgorithm_3 \ if (snoop) \ { \ strcpy (ccSnoop, relative_to_infinity ((dq-1)->CC)); \ strcpy (ddSnoop, relative_to_infinity ((dq-1)->DD)); \ fprintf (stderr, "[%3u,%3u]" \ " [" unsposFmt "].cc=%s" \ " [" unsposFmt "].DD=%s" \ " dq=%d\n", \ 0, col, col, ccSnoop, col, ddSnoop, col); \ } #define snoopAlgorithm_4A \ if (snoop) \ fprintf (stderr, "\nrow " unsposFmt \ " LY=" unsposFmt \ " RY=" unsposFmt, \ row, LY, RY); #define snoopAlgorithm_4B \ if (snoop) \ fprintf (stderr, " -> L=" sgnposFmt \ " R=" sgnposFmt \ " LY=" unsposFmt \ " RY=" unsposFmt, \ L, R, LY, RY); #define snoopAlgorithm_4C \ if (snoop) \ fprintf (stderr, " tbNeeded=%d dp=%ld dq=%ld\n", \ tbNeeded, (long) (dp-dynProg->p), (long) (dq-dynProg->p)); #define snoopAlgorithm_5 \ if (snoop) \ { \ strcpy (ccSnoop, relative_to_infinity (c)); \ strcpy (ddSnoop, relative_to_infinity (dp->DD)); \ strcpy (iiSnoop, relative_to_infinity (i)); \ fprintf (stderr, "[%3u,%3u] %c%c B[%ld]=%c" \ " c=%s d=%s i=%s" \ " dp=%ld dq=%ld", \ row, col, A[row], B[col], (long) (b-B), *b, \ ccSnoop, ddSnoop, iiSnoop, \ (long) (dp-dynProg->p), (long) (dq-dynProg->p)); \ } #define snoopAlgorithm_5A \ if (snoop) \ { \ strcpy (ccSnoop, relative_to_infinity (bestScore)); \ fprintf (stderr, " **new best score %s**", ccSnoop); \ } #define snoopAlgorithm_5B \ if (snoop) \ fprintf (stderr, "\n"); #define snoopAlgorithm_6 \ if (snoop) \ { \ strcpy (ccSnoop, relative_to_infinity ((dq-1)->CC)); \ strcpy (ddSnoop, relative_to_infinity ((dq-1)->DD)); \ fprintf (stderr, " [%lu].CC=%s [%lu].DD=%s link=%s\n", \ (long) ((dq-1)-dynProg->p), ccSnoop, \ (long) ((dq-1)-dynProg->p), ddSnoop, \ link_to_string (link)); \ } #define snoopAlgorithm_7A \ if (snoop) \ { \ if ((rightSeg != NULL) && (R > 0)) \ fprintf (stderr, "NN <- " sgnposFmt " (from R)\n", R-1); \ else \ fprintf (stderr, "NN <- " sgnposFmt " (from N)\n", (sgnpos)N); \ } #define snoopAlgorithm_7B \ if (snoop) \ fprintf (stderr, "RY <- " unsposFmt \ " (hit ydrop prior to RY)\n", \ RY); #define snoopAlgorithm_7C \ if (snoop) \ prevRY = RY; #define snoopAlgorithm_7D \ if ((snoop) && (RY != prevRY)) \ fprintf (stderr, "RY <- " unsposFmt \ " (feasible overhang)\n", \ RY); #define snoopAlgorithm_7E \ if (snoop) \ fprintf (stderr, "RY <- " unsposFmt \ " (room at right boundary)\n", \ RY); #define snoopAlgorithm_8 \ if (snoop) \ fprintf (stderr, "alignment ends at [%3u,%3u]\n", *_end1, *_end2); #endif // snoopAlgorithm #ifdef snoopAlgorithmTrap static unspos trapRow = 497; static unspos trapCol = (unspos) -1; #endif // snoopAlgorithmTrap //=== stuff for snoopTraceback === #ifndef snoopTraceback #define snoopTraceback_1 ; #define snoopTraceback_2 ; #define snoopTraceback_3 ; #endif // not snoopTraceback #ifdef snoopTraceback #define snoopTraceback_1 \ fprintf (stderr, "(tbp=%08X) tb[%3d,%3d] <- %-2d %-9s\n", \ (u32) tbp-1, \ row, col, link, link_to_string(link)); #define snoopTraceback_2 \ fprintf (stderr, "(tbp=%08X) tb[%3d,%3d] is %-2d %-9s" \ " prevOp=%s --> op=%d", \ (u32) &tb->space[tbRow[row] + col], \ row, col, link, link_to_string(link), \ op_string(prevOp), op); #define snoopTraceback_3 \ fprintf (stderr, " -> %s -> row %d, col %d\n", \ op_string(op), row, col); #ifndef snoopAlgorithm include_link_to_string #endif //snoopAlgorithm #endif // snoopTraceback //=== stuff for snoopSubprobs === #ifndef snoopSubprobs #define snoopSubprobsB_1 ; #define snoopSubprobsB_2 ; #define snoopSubprobsB_3 ; #endif // not snoopSubprobs #ifdef snoopSubprobs #define snoopSubprobsB_1 \ if (leftSeg == NULL) \ fprintf (stderr, "\n leftSeg=(none)"); \ else \ fprintf (stderr, "\n leftSeg=%s" \ " " unsposSlashFmt \ " " unsposSlashFmt, \ (leftSeg->type==diagSeg)? "diag" \ : (leftSeg->type==horzSeg)? "horz" \ : (leftSeg->type==vertSeg)? "vert" \ : "????", \ leftSeg->b1, leftSeg->b2, \ leftSeg->e1, leftSeg->e2); \ if (rightSeg == NULL) \ fprintf (stderr, "\n rightSeg=(none)"); \ else \ fprintf (stderr, "\n rightSeg=%s" \ " " unsposSlashFmt \ " " unsposSlashFmt, \ (rightSeg->type==diagSeg)? "diag" \ : (rightSeg->type==horzSeg)? "horz" \ : (rightSeg->type==vertSeg)? "vert" \ : "????", \ rightSeg->b1, rightSeg->b2, \ rightSeg->e1, rightSeg->e2); \ if (!reversed) \ fprintf (stderr, "\n L=" sgnposFmt " R=" sgnposFmt, \ L, R); \ else \ { \ sgnpos tempL = L; \ sgnpos tempR = R; \ sgnpos tempT = 0; \ if ((leftSeg == NULL) && (rightSeg != NULL)) { tempL = -R+1; tempR = N+1; } \ else if ((leftSeg != NULL) && (rightSeg == NULL)) { tempR = -L-1; tempL = 0; } \ else if ((leftSeg != NULL) && (rightSeg != NULL)) { tempT = -L-1; tempL = -R+1; tempR = tempT; } \ fprintf (stderr, "\n L=" sgnposFmt " R=" sgnposFmt \ " -> L=" sgnposFmt " R=" sgnposFmt, \ L, R, tempL, tempR); \ } \ fprintf (stderr, "\n ----"); #define snoopSubprobsB_2 \ fprintf (stderr, "\n row " unsposFmt \ " -> L=" sgnposFmt \ " R=" sgnposFmt \ " LY=" unsposFmt \ " RY=" unsposFmt, \ row, L, R, LY, RY); #ifndef collect_stats #define snoopSubprobsB_3 \ fprintf (stderr, "\n"); \ fprintf (stderr, "\trows=" unsposFmt, row); #endif // not collect_stats #ifdef collect_stats #define snoopSubprobsB_3 \ fprintf (stderr, "\n"); \ fprintf (stderr, "\trows=" unsposFmt, row); \ fprintf (stderr, "\tmaxDpRows=%s", commatize(gappedExtendStats.maxDpRows)); #endif // collect_stats #endif // snoopSubprobs //=== stuff for snoopBounds === #ifndef snoopBounds #define debugSnoopBounds_1 ; #define debugSnoopBounds_2 ; #define debugSnoopBounds_3 ; #endif // not snoopBounds #ifdef snoopBounds static char* segmentTypeName[3] = {"diag", "horz", "vert"}; #define debugSnoopBounds_1 \ { \ fprintf (stderr, "bounds: leftSeg(%s)=(" unsposFmt ")/" unsposFmt \ " anchor=(" unsposFmt ")/" unsposFmt "\n", \ segmentTypeName[(u8)leftSeg->type], \ leftSeg->b1, leftSeg->b2, \ anchor1, anchor2); \ if (leftSeg->type == diagSeg) \ fprintf (stderr, "bounds: L=(" unsposFmt "-" unsposFmt ")" \ "-(" unsposFmt "-" unsposFmt ")" \ "=" sgnposFmt "\n", \ leftSeg->b2, anchor2, leftSeg->b1, anchor1, L); \ else \ fprintf (stderr, "bounds: L=(" unsposFmt "-" unsposFmt ")" \ "=" sgnposFmt "\n", \ leftSeg->b2, anchor2, L); \ } #define debugSnoopBounds_2 \ { \ fprintf (stderr, "bounds: rightSeg(%s)=(" unsposFmt ")/" unsposFmt \ " anchor=(" unsposFmt ")/" unsposFmt "\n", \ segmentTypeName[(u8)rightSeg->type], \ rightSeg->b1, rightSeg->b2, \ anchor1, anchor2); \ if (rightSeg->type == diagSeg) \ fprintf (stderr, "bounds: R=(" unsposFmt "-" unsposFmt ")" \ "-(" unsposFmt "-" unsposFmt ")" \ "=" sgnposFmt "\n", \ rightSeg->b2, anchor2, rightSeg->b1, anchor1, R); \ else \ fprintf (stderr, "bounds: R=(" unsposFmt "-" unsposFmt ")" \ "=" sgnposFmt "\n", \ rightSeg->b2, anchor2, R); \ } #define debugSnoopBounds_3 \ fprintf (stderr, "bounds: swapped to L=" sgnposFmt " R=" sgnposFmt "\n", \ L, R); #endif // snoopBounds //=== traceback row memory (see note (3)) === // // minTbRowsNeeded is the inverse of the round_up result in tbrow_needed #define minTbRowsNeeded (((((512*1024)/sizeof(u32))-1)*16)/17) static u32* tbRow = NULL; // memory to track row positions in static u32 tbRowLen = 0; // .. traceback array static void tbrow_needed (u32 rowsNeeded); static void tbrow_needed (u32 rowsNeeded) { size_t needed; if (rowsNeeded <= tbRowLen) return; needed = round_up(sizeof(u32)*(rowsNeeded+1+rowsNeeded/16), 512*1024); tbRow = realloc_or_die ("ydrop_one_sided_align tbRow", tbRow, needed); tbRowLen = needed / sizeof(u32); } void free_traceback_rows (void) { free_if_valid ("free_traceback_rows", tbRow); tbRow = NULL; tbRowLen = 0; } //=== ydrop_one_sided_align === static score ydrop_one_sided_align (alignio* io, int reversed, // true => search to the lower-left u8* A, // vertical sequence text u8* B, // horizontal sequence text unspos M, // vertical sequence length unspos N, // horizontal sequence length int trimToPeak, editscript** script, unspos* _end1, unspos* _end2) { tback* tb; // tb is space provided to record traceback; u8* tbp; // .. tbp is the current traceback cell int tbLen; // .. and steps linearly from tbp->space; int tbNeeded; // .. tbRow[r] is the conceptual start of // .. the traceback cells for row r; it is // .. indexed as tbRow[r][c] for c=LY..R', // .. where R' is the last cell used in // .. row r unspos anchor1, anchor2; // anchor positions in A and B scorerow* allSub; // substitution scores matrix score* sub; // substitution scores vector, current row score gapOE, gapE; // gap penalties score yDrop; // score drop threshold int yDropTail; // length of shortest score fall >= yDrop // nota bene: row, col, leftCol, L, R, // .. LY, RY, prevLY, NN, npCol, end1 // .. and end2 are all relative to the DP // .. matrix unspos row, col = 0; // current DP grid point unspos leftCol; // (copy of left column, for stats) sgnpos L, R; // external column limits for current row unspos LY, RY; // actual column limits for current row; // .. ("Y" is for y-drop, not cartesian Y) unspos prevLY; // left column limit for previous row sgnpos NN; // truncated right side bound, current row unspos npCol; // last non-pruned cell in current row unspos end1, end2; // end of optimal alignment int endIsBoundary; // true => report boundaryScore instead of // .. bestScore score bestScore; // score of best alignment seen so far score boundaryScore; // score of best alignment seen so far that // .. ends at the end of either sequence dpMatrix* dynProg; // DP cells dpCell* dp; // scans previous row of dp matrix dpCell* dq; // scans current row of dp matrix u8* b; // scans horizontal sequence score c, d, i; // running scores for DP cells score cOpen, cNext, cTemp;// scratch values for cell scores activeseg* active; // list of segments that intersect the // .. sweep row within the feasible region galign* alignList; galign* rightAlign, *leftAlign; aliseg* rightSeg, *leftSeg; u8 link = 0; // traceback link u8 op, prevOp; // edit operations #if ((defined snoopAlgorithm) || (defined snoopAlgorithmTrap)) int snoop = true; #ifdef snoopAlgorithm unspos prevRY = 0; #endif // snoopAlgorithm #endif // snoopAlgorithm or snoopAlgorithmTrap // sanity check; if either sequence is empty there's no alignment if ((N <= 0) || (M <= 0)) { *(_end1) = *(_end2) = 0; return 0; } #if ((defined snoopAnchors) || ((defined snoopAlgorithm) && (defined debugPosA1))) fprintf (stderr,"[hspId=" u64Fmt "] ydrop_one_sided_align(" unsposSlashFmt ") %s\n", io->hspId, io->anchor1, io->anchor2, (reversed)?" reversed":""); #endif // snoopAnchors #ifdef snoopAlgorithm #if ((defined debugPosA1) && (defined debugPosB1)) snoop = ((io->anchor1 == debugPosA1) && (io->anchor2 == debugPosA2)) || ((io->anchor1 == debugPosB1) && (io->anchor2 == debugPosB2)); #elif (defined debugPosA1) snoop = (io->anchor1 == debugPosA1) && (io->anchor2 == debugPosA2); #elif (defined debugPosB1) snoop = (io->anchor1 == debugPosB1) && (io->anchor2 == debugPosB2); #endif // debugPosA1,debugPosB1 #endif // snoopAlgorithm gapped_extend_count_stat (numExtensions); dbg_timing_count_stat (numExtensions); snoopAlgorithm_1; // extract scoring constants allSub = io->scoring->sub; gapE = io->scoring->gapExtend; gapOE = io->scoring->gapOpen + gapE; // (any new gap gets both penalties) yDrop = io->yDrop; tb = io->tb; tbLen = tb->size; if (gapE != 0) yDropTail = (yDrop/gapE) + 6; else { // when gapE is zero, the above results would be infinite; but we can // limit yDropTail to the distance from the length of the sequence; // this can increase the amount of memory needed; the solution here is // not completely sufficient; "truncating alignment" reports are still // likely int maxYDropTail = 500*1000; if (N < (unsigned int) maxYDropTail) yDropTail = N+1; else yDropTail = maxYDropTail; } // determine initial left and right constraints L = 0; R = N+1; // (in blastz this was R=N) anchor1 = io->anchor1; anchor2 = io->anchor2; leftSeg = io->leftSeg; if (leftSeg != NULL) { L = signed_difference (leftSeg->b2, anchor2); if (leftSeg->type == diagSeg) L -= signed_difference (leftSeg->b1, anchor1); debugSnoopBounds_1; } rightSeg = io->rightSeg; if (rightSeg != NULL) { R = signed_difference (rightSeg->b2, anchor2); if (rightSeg->type == diagSeg) R -= signed_difference (rightSeg->b1, anchor1); debugSnoopBounds_2; } snoopSubprobsB_1; // if we're doing a reversed alignment we need to swap the L-R bounds (see // note (14)) if (reversed) { sgnpos temp = 0; // (placate compiler) if ((leftSeg == NULL) && (rightSeg != NULL)) { L = -R+1; R = N+1; } else if ((leftSeg != NULL) && (rightSeg == NULL)) { R = -L-1; L = 0; } else if ((leftSeg != NULL) && (rightSeg != NULL)) { temp = -L-1; L = -R+1; R = temp; } debugSnoopBounds_3; } active = NULL; rightAlign = io->rightAlign; leftAlign = io->leftAlign; alignList = (!reversed)? io->aboveList : io->belowList; // make sure we have a reasonable number of traceback rows to start with // (see note (3)) tbrow_needed (minTbRowsNeeded); tbRow[0] = 0; tbp = tb->space; ////////// // compute first row of dp matrix ////////// // make sure we have enough traceback and dp space for the first row tbNeeded = yDropTail; if (tbNeeded > tbLen) suicide ("not enough space in trace_back array"); dynProg = zalloc_or_die ("ydrop_one_sided_align dynProg", sizeof(dpMatrix)); dp_ready (dynProg, tbNeeded); gapped_extend_count_stat (zallocCallsB); gapped_extend_add_stat (zallocTotalB, sizeof(dpMatrix)); // compute first row of dp matrix (see notes (8), (10) and (13)) dq = dynProg->p; dq->CC = cTemp = 0; // set C[0][0] c = (dq++)->DD = -gapOE; // set D[1][0] *(tbp++) = 0; snoopAlgorithm_2; for (col=1 ; (col<=N)&&(cTemp>=-yDrop) ; col++) { dq->CC = cTemp = c; // set C[0][col] (dq++)->DD = c - gapOE; // set D[1][col] c -= gapE; *(tbp++) = cFromI; snoopAlgorithm_3; } gapped_extend_add_stat (dpCellsVisited, col); LY = 0; RY = col; // (1 column beyond the feasible region) ////////// // compute additional rows of DP matrix ////////// end1 = end2 = 0; bestScore = 0; boundaryScore = negInf; endIsBoundary = false; for (row=1; row<=M ; row++) { #ifdef snoopAlgorithm if (snoop) cTemp = 0; // (place to set a breakpoint) #endif // snoopAlgorithm snoopAlgorithm_4A; #ifdef snoopAlgorithmTrap if ((snoop) && (row == trapRow) && (trapCol == (unspos) -1)) cTemp = 0; #endif // update sweep row bounds, active segments, masking prevLY = LY; update_LR_bounds (reversed, &rightSeg, &leftSeg, &rightAlign, &leftAlign, row, anchor1, anchor2, &L, &R, &LY, &RY); update_active_segs (reversed, &active, &alignList, dynProg->p-prevLY, row, anchor1, anchor2, LY, RY); snoopAlgorithm_4B; snoopSubprobsB_2; // make sure we have enough traceback and dp space for this row (see // note (3)) tbrow_needed (row+1); gapped_extend_max_stat (maxDpRows, row+1); if (RY < LY) RY = LY; // (see note 11) tbNeeded = RY - LY + yDropTail; if ((tbp - tb->space) + tbNeeded >= tbLen) { if (gapped_extend_inhibitTruncationReport) goto dp_finished; if (!reversed) fprintf (stderr, "truncating alignment ending at (" unsposCommaFmt ");", end1 + anchor1 + 1, end2 + anchor2 + 1); else fprintf (stderr, "truncating alignment starting at (" unsposCommaFmt ");", anchor1 + 2 - end1, anchor2 + 2 - end2); fprintf(stderr, " anchor at (" unsposCommaFmt ")\n", anchor1, anchor2); if (!haveReportedTruncation) { haveReportedTruncation = true; fprintf(stderr, "truncation can be reduced by using --allocate:traceback to increase traceback memory\n"); } goto dp_finished; } tbRow[row] = (tbp - tb->space) - LY; // set up DP pointers for this sweep row (see note (5)) dp_ready (dynProg, tbNeeded); // make sure we have enough DP space dq = dynProg->p; // dq cells start at col == LY dp = dq + LY - prevLY; // dp cells start at col == prevLY snoopAlgorithm_4C; // compute DP values for all bounded columns in this row (see note (4)) sub = allSub[A[row]]; col = leftCol = LY; b = B + col + 1; // (b scans horizontal sequence, one column ahead) npCol = col; // npCol records the last non-pruned position i = negInf; // 'set' I[row][col] c = negInf; // propose C[row][col] for ( ; (colDD; // get D[row][col] // at this point d, i and c contain the DP values for the cell at // [row][col]; the latter is the value for reaching C from previous // grid points, but since C also has edges entering it from the // current D and I nodes, we might improve it // nota bene: when we *can* improve c, we make an arbitrary choice // to prefer deletion to insertion (when i and d are equal) // nota bene 2: all paths through this series of ifs assign a value // to link if ((active != NULL) && (dp->mask == row)) { snoopAlgorithm_5; prune; snoopAlgorithm_5B; continue; } if ((d > c) || (i > c)) // === we CAN improve C === { // nota bene: both iExtend and dExtend are set here because // the value of the C and I (or C and D) are equal, so traceback // may as well take a gap extend into this cell if (d >= i) { c = d; link = cFromD | iExtend | dExtend; } else { c = i; link = cFromI | iExtend | dExtend; } snoopAlgorithm_5; if (c < bestScore - yDrop) { prune; snoopAlgorithm_5B; continue; } #ifndef allowBackToBackGaps // not allowing back-to-back gaps, so we don't need to consider // opening a gap here i -= gapE; // 'set' I[row][col+1] dq->DD = d - gapE; // set D[row+1][col] #else // back-to-back gaps are allowed, so we must consider gap opens cOpen = c - gapOE; d -= gapE; // set D[row+1][col] if (cOpen > d) { dq->DD = cOpen; link &= ~dExtend; } else dq->DD = d; i -= gapE; // 'set' I[row][col+1] if (cOpen > i) { i = cOpen; link &= ~iExtend; } #endif // allowBackToBackGaps } else // === we CANNOT improve C === { snoopAlgorithm_5; if (c < bestScore - yDrop) { prune; snoopAlgorithm_5B; continue; } if (c >= bestScore) { bestScore = c; end1 = row; end2 = col; endIsBoundary = false; snoopAlgorithm_5A; } if ((!trimToPeak) && (c >= boundaryScore) && ((row == M) || (col == N))) { boundaryScore = c; end1 = row; end2 = col; endIsBoundary = true; } cOpen = c - gapOE; d -= gapE; // set D[row+1][col] if (cOpen > d) { dq->DD = cOpen; link = cFromC; } else { dq->DD = d; link = cFromC | dExtend; } i -= gapE; // 'set' I[row][col+1] if (cOpen > i) i = cOpen; else link |= iExtend; } npCol = col; // save as last non-pruned position // save C for this column, and compute proposed C for the next // column (see note (6)) cNext = (dp++)->CC+sub[*(b++)]; // propose C[row][col+1] (dq++)->CC = c; // set C[row][col] c = cNext; *(tbp++) = link; // set link into C[row][col] snoopAlgorithm_6; //snoopTraceback_1; } gapped_extend_add_stat (dpCellsVisited, col-leftCol); // if the feasible region is empty, we're done if (LY >= RY) goto dp_finished; // finish up this row, by either moving the right bound left or // prolonging the row to support an overhang on the row above snoopAlgorithm_7A; NN = ((rightSeg != NULL) && (R > 0))? (R-1) : ((sgnpos) N); if (RY > npCol+1) // we hit ydrop prior to RY { RY = npCol+1; snoopAlgorithm_7B; } else { // the current row reached its right bound, but the row above may // still have a feasible overhang so we prolong this row with // insertions (see note (9)) snoopAlgorithm_7C; while ((i >= bestScore - yDrop) && (((sgnpos)RY) <= NN)) { if (((u32)(dq - dynProg->p)) >= dynProg->len) suicidef("(in ydrop_one_sided_align:%d, dq-dynProg->p==%d, dynProg->len=" unsposFmt ")", __LINE__, dq - dynProg->p, dynProg->len); dq->CC = i; // set C[row][col] (dq++)->DD = i - gapOE; // set D[row+1][col] i -= gapE; // 'set' I[row][col+1] *(tbp++) = cFromI; RY++; } snoopAlgorithm_7D; } // terminate the cell at the right boundary, so that nothing will // step from it (termination occurs if the if is false and we thus // *fail* to increment RY) if (((sgnpos)RY) <= NN) { if (((u32)(dq - dynProg->p)) >= dynProg->len) suicidef("(in ydrop_one_sided_align:%d, dq-dynProg->p==%d, dynProg->len=" unsposFmt ")", __LINE__, dq - dynProg->p, dynProg->len); dq->DD = dq->CC = negInf; // set D[row+1][col] RY++; // .. and set C[row][col] snoopAlgorithm_7E; } } dp_finished: snoopSubprobsB_3; *(_end1) = row = end1; *(_end2) = col = end2; snoopAlgorithm_8; ////////// // traceback the alignment to create the edit script ////////// #ifdef snoopAlgorithm if (snoop) cTemp = 0; // (place to set a breakpoint) #endif // snoopAlgorithm for (prevOp=0 ; (row>=1) || (col>0) ; prevOp=op) { link = tb->space[tbRow[row] + col]; op = link & cidBits; if ((prevOp == cFromI) && ((link & iExtend) != 0)) op = cFromI; if ((prevOp == cFromD) && ((link & dExtend) != 0)) op = cFromD; snoopTraceback_2 if (op == cFromI) { col--; edit_script_ins (script, 1); } else if (op == cFromD) { row--; edit_script_del (script, 1); } else { row--; col--; edit_script_sub (script, 1); } snoopTraceback_3 } filter_active_segs (&active, 2); // (disposes of everything in the list) free_if_valid ("ydrop_one_sided_align dynProg->p", dynProg->p); free_if_valid ("ydrop_one_sided_align dynProg", dynProg); if (endIsBoundary) return boundaryScore; else return bestScore; } //---------- // // dp_ready-- // Allocate the dynamic-programming "sweep row" array, and ensure that the // first n elements exist and have been initialized. // //---------- // // Arguments: // dpMatrix* dynProg: The current sweep array. // unspos needed: The number of elements required. // // Returns: // nothing; the contents of dynProg are (potentially) modified by adding // enough empty cells to satisfy the number needed // //---------- static void dp_ready (dpMatrix* dynProg, unspos needed) { u32 oldLen; int added; if (dynProg == NULL) suicide ("dp_ready called with NULL pointer"); if (needed > 0xFFFFFFFF) suicidef ("in dp_ready, number of DP cells needed exceeds limit\n" " " unsposFmt " cells requested", needed); // make sure that there's room for n cells if (dynProg->p == NULL) { oldLen = 0; dynProg->len = needed + 1000; dynProg->p = malloc_or_die ("dp_ready", dynProg->len * sizeof(dpCell)); } else { oldLen = dynProg->len; if (needed > oldLen) { dynProg->len = needed + dynProg->len/16 + 1000; dynProg->p = realloc_or_die ("dp_ready", dynProg->p, dynProg->len * sizeof(dpCell)); } } // 0..oldLen-1 are already in use, but zero the rest to clear their mask // fields added = dynProg->len - oldLen; if (added > 0) memset (dynProg->p + oldLen, 0, sizeof(dynProg->p[0])*added); gapped_extend_max_stat (maxDpColumns, needed); } //---------- // // msp_left_right-- // Given MSP m, we search obi, the list of alignments ordered by beginning // point (in vertical/seq1), to find the alignments and their contained // segments (gap-free pieces) that are the closest to the left and right of // the MSP's anchor-point. // //---------- // // Arguments: // galign* obi: The list of (previously computed) alignments, ordered by // .. increasing beginning point. This may be NULL. // galign* m: The MSP to check against the list. // // Returns: // false if the anchor-point is in an already-computed alignment; true // otherwise; actual result values (the segment pointers) are deposited in // the MSP. // //---------- static int msp_left_right (galign* obi, galign* m) { unspos pos1 = m->pos1; unspos pos2 = m->pos2; unspos right, left; galign* mRight, *mLeft; aliseg* bRight, *bLeft; aliseg* bp; sgnpos x; right = left = seqposInfinity; mRight = mLeft = NULL; bRight = bLeft = NULL; // process all alignments in the obi that overlap (along vertical axis) // m's anchor point for ( ; (obi!=NULL)&&(obi->pos1<=pos1) ; obi=obi->next) { if (obi->end1 < pos1) continue; // invariant: obi->pos1 <= pos1 <= obi->end1 for (bp=obi->firstSeg ; bp!=NULL ; bp=bp->nextSeg) { if (bp->e1 >= pos1) break; } if (bp == NULL) continue; // invariant: bp is the first segment (in this alignment) that // intersects the line y=pos1 if (bp->type == horzSeg) suicide ("msp_left_right: cannot be horizontal"); if (bp->type == diagSeg) x = signed_difference (bp->b2, pos2) // x is how far to the + signed_difference (pos1, bp->b1); // .. right of m the else // (bp->type == vertSeg) // .. segment is, along x = signed_difference (bp->b2, pos2); // .. the line y=pos1 if (x == 0) return false; if ((x > 0) && ((unspos) x < right)) { // a new closest segment to the right right = (unspos) x; mRight = obi; bRight = bp; } else if ((x < 0) && ((unspos) -x < left)) { // a new closest segment to the left left = (unspos) -x; mLeft = obi; bLeft = bp; } } m->rightAlign1 = m->rightAlign2 = mRight; m->rightSeg1 = m->rightSeg2 = bRight; m->leftAlign1 = m->leftAlign2 = mLeft; m->leftSeg1 = m->leftSeg2 = bLeft; return true; } //---------- // // get_above_below-- // In preparation for extending the anchor-point of an MSP, find the closest // alignment ending below the anchor-point and the closest alignment starting // above the anchor-point. // //---------- // // Arguments: // alignio* io: io->anchor1 is the anchor's vertical/seq1 position; // this routine fills in io->belowList and io->aboveList // galign* obi: The list of (previously computed) alignments, ordered by // .. increasing beginning point. This may be NULL. // galign* oed: The list of (previously computed) alignments, ordered by // .. decreasing ending point. This may be NULL. // // Returns: // nothing; values are written to io->belowList and io->aboveList. // //---------- static void get_above_below (alignio* io, galign* obi, galign* oed) { unspos pos1 = io->anchor1; galign* mp; for (mp=oed ; mp!=NULL ; mp=mp->prev) { if (mp->end1 < pos1) break; } io->belowList = mp; for (mp=obi ; mp!=NULL ; mp=mp->next) { if (mp->pos1 > pos1) break; } io->aboveList = mp; } //---------- // // align_left_right-- // Given a new alignment, determine the alignments and their segments that are // the closest in either horizontal direction at both ends of the alignment. // //---------- // // Arguments: // galign* obi: The list of (previously computed) alignments, ordered by // .. increasing beginning point. This may be NULL. // galign* m: The alignment to check. // // Returns: // (nothing) // //---------- static void align_left_right (galign* obi, galign* m) { unspos pos1 = m->pos1, pos2 = m->pos2; unspos end1 = m->end1, end2 = m->end2; unspos rightOfBottom, rightOfTop, leftOfBottom, leftOfTop; galign* mRightOfBottom, *mRightOfTop, *mLeftOfBottom, *mLeftOfTop; aliseg* bRightOfBottom, *bRightOfTop, *bLeftOfBottom, *bLeftOfTop; aliseg* bp; sgnpos x; rightOfBottom = rightOfTop = leftOfBottom = leftOfTop = seqposInfinity; mRightOfBottom = mLeftOfBottom = mRightOfTop = mLeftOfTop = NULL; bRightOfBottom = bRightOfTop = bLeftOfBottom = bLeftOfTop = NULL; for ( ; obi!=NULL ; obi=obi->next) { if ((obi->pos1 > end1) || (obi->end1 < pos1)) continue; // invariant: obi->pos1 <= end1 and obi->end1 >= pos1 // meaning: obi overlaps m along vertical/seq1 for (bp=obi->firstSeg ; bp!=NULL ; bp=bp->nextSeg) { if ((bp->type != horzSeg) && (bp->e1 >= pos1)) break; } // invariant: bp is either NULL or the first non-horizontal segment (in // this alignment) that intersects the line y=pos1 if ((bp != NULL) && (bp->b1 <= pos1)) { if (bp->type == diagSeg) x = signed_difference (bp->b2, pos2) // x is how far to the + signed_difference (pos1, bp->b1); // .. right of m the else // (bp->type == vertSeg) // .. segment is, along x = signed_difference (bp->b2, pos2); // .. the line y=pos1 if ((x > 0) && ((unspos) x < rightOfBottom)) { rightOfBottom = (unspos) x; mRightOfBottom = obi; bRightOfBottom = bp; } else if ((x < 0) && ((unspos) -x < leftOfBottom)) { leftOfBottom = (unspos) -x; mLeftOfBottom = obi; bLeftOfBottom = bp; } } for ( ; bp!=NULL ; bp=bp->nextSeg) { if (bp->type != horzSeg && bp->e1 >= end1) break; } if ((bp != NULL) && (bp->type != horzSeg) && (bp->e1 >= end1)) { if (bp->type == diagSeg) x = signed_difference (bp->b2, end2) // x is how far to the + signed_difference (end1, bp->b1); // .. right of m the else // (bp->type == vertSeg) // .. segment is, along x = signed_difference (bp->b2, end2); // .. the line y=end1 if ((x > 0) && ((unspos) x < rightOfTop)) { rightOfTop = (unspos) x; mRightOfTop = obi; bRightOfTop = bp; } else if ((x < 0) && ((unspos) -x < leftOfTop)) { leftOfTop = (unspos) -x; mLeftOfTop = obi; bLeftOfTop = bp; } } } m->rightAlign1 = mRightOfBottom; m->rightSeg1 = bRightOfBottom; m->rightAlign2 = mRightOfTop; m->rightSeg2 = bRightOfTop; m->leftAlign1 = mLeftOfBottom; m->leftSeg1 = bLeftOfBottom; m->leftAlign2 = mLeftOfTop; m->leftSeg2 = bLeftOfTop; } //---------- // // insert_align-- // Insert a new alignment into the two lists of alignments, one ordered by // increasing beginning point (in vertical/seq1) and the other ordered by // decreasing end point. // //---------- // // Arguments: // galign* m: The MSP that was used to anchor the alignment. // galign** obi: Pointer to the list of (previously computed) alignments, // .. ordered by increasing beginning point. The list may // .. be empty (*obi == NULL). This is updated by this // .. function. // galign** oed: Pointer to the list of (previously computed) alignments, // .. ordered by decreasing ending point. The list may // .. be empty (*oed == NULL). This is updated by this // .. function. // // Returns: // (nothing) // //---------- //=== stuff for snoopBlocks === #ifndef snoopBlocks #define debugSnoopBlocks_6 ; #endif // not snoopBlocks #ifdef snoopBlocks #define debugSnoopBlocks_6 \ fprintf (stderr, "adding alignment block [%8p]" \ " b " unsposFmt " " unsposFmt \ " e " unsposFmt " " unsposFmt \ " s " scoreFmtSimple "\n", \ m, m->pos1, m->pos2, m->end1, m->end2, \ (m->align == NULL)? ((score) 0) : m->align->s); #endif // snoopBlocks static void insert_align (galign* m, galign** _obi, galign** _oed) { galign* mp; // (mp scans list, galign* mq; // .. mq is predecessor in scan direction) galign* obi = *_obi; galign* oed = *_oed; if (m->firstSeg == NULL) suicide ("insert_align: null first segment"); debugSnoopBlocks_6; for (mq=NULL,mp=obi ; mp!=NULL ; mq=mp,mp=mp->next) { if (mp->pos1 >= m->pos1) break; } if (mq != NULL) { mq->next = m; m->next = mp; } else { m->next = obi; obi = m; } for (mq=NULL,mp=oed ; mp!=NULL ; mq=mp,mp=mp->prev) { if (mp->end1 <= m->end1) break; } if (mq != NULL) { mq->prev = m; m->prev = mp; } else { m->prev = oed; oed = m; } *(_obi) = obi; *(_oed) = oed; } //---------- // // update_LR_bounds-- // As we move to a new DP row, update LY and RY, the actual column limits // for dynamic programming. // //---------- // // Arguments: // int reversed: true => the DP row advances downward // false => it advances upward // aliseg** rightSeg: Current constraining segments; these may be // aliseg** leftSeg: .. updated by this function. // galign** rightAlign: Alignments containing those segments; these may be // galign** leftAlign: .. updated by this function. // unspos row: The sweep row. // unspos anchor1: The position at which the alignment // unspos anchor2: .. began. // sgnpos* L, (pointer to) most recent bounding constraints; // sgnpos* R: .. these may be updated by this function. // unspos* LY: (pointer to) most recent Y-drop constraints; these // unspos* RY: .. may lie strictly inside L and R; these may be // .. updated by this function. // // Returns: // (nothing) // //---------- #if ((!defined snoopAlgorithm)&&(!defined snoopSubprobs)) #define snoopLR_forward_left_1 ; #define snoopLR_forward_left_2 ; #define snoopLR_forward_left_3 ; #define snoopLR_forward_right_1 ; #define snoopLR_forward_right_2 ; #define snoopLR_forward_right_3 ; #define snoopLR_reverse_left_1 ; #define snoopLR_reverse_left_2 ; #define snoopLR_reverse_left_3 ; #define snoopLR_reverse_right_1 ; #define snoopLR_reverse_right_2 ; #define snoopLR_reverse_right_3 ; #endif // not snoopAlgorithm and not snoopSubprobs #if (defined snoopAlgorithm) #define snoopLR_forward_left_1 \ if (snoop) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ if ((*leftSeg)->type == diagSeg) \ fprintf (stderr, "update_LR: L <- " sgnposFmt \ " (on same diag seg)\n", \ L); \ else \ fprintf (stderr, "update_LR: L <- " sgnposFmt \ " (on same vert seg)\n", \ L); \ } #define snoopLR_forward_left_2 \ if ((snoop) && (*leftSeg != NULL)) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ fprintf (stderr, "update_LR: L <- " sgnposFmt \ " (on new seg " unsposSlashFmt " " unsposSlashFmt \ ")\n", \ L, \ (*leftSeg)->b1, (*leftSeg)->b2, \ (*leftSeg)->e1, (*leftSeg)->e2); \ } #define snoopLR_forward_left_3 \ if ((snoop) && (L > (sgnpos) LY)) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ fprintf (stderr, "update_LR: LY <- " unsposFmt \ "\n", \ (unspos) max ((sgnpos) LY, L)); \ } #define snoopLR_forward_right_1 \ if (snoop) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ if ((*rightSeg)->type == diagSeg) \ fprintf (stderr, "update_LR: L <- " sgnposFmt \ " (on same diag seg)\n", \ R); \ else \ fprintf (stderr, "update_LR: R <- " sgnposFmt \ " (on same vert seg)\n", \ R); \ } #define snoopLR_forward_right_2 \ if ((snoop) && (*rightSeg != NULL)) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ fprintf (stderr, "update_LR: R <- " sgnposFmt \ " (on new seg " unsposSlashFmt " " unsposSlashFmt \ ")\n", \ R, \ (*rightSeg)->b1, (*rightSeg)->b2, \ (*rightSeg)->e1, (*rightSeg)->e2); \ } #define snoopLR_forward_right_3 \ if ((snoop) && (R < (sgnpos) RY)) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ fprintf (stderr, "update_LR: RY <- " unsposFmt \ "\n", \ special_min (RY, R)); \ } #define snoopLR_reverse_left_1 \ if (snoop) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ if ((*rightSeg)->type == diagSeg) \ fprintf (stderr, "update_LR: L <- " sgnposFmt \ " (on same diag seg)\n", \ L); \ else \ fprintf (stderr, "update_LR: L <- " sgnposFmt \ " (on same vert seg)\n", \ L); \ } #define snoopLR_reverse_left_2 \ if ((snoop) && (*rightSeg != NULL)) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ fprintf (stderr, "update_LR: L <- " sgnposFmt \ " (on new seg " unsposSlashFmt " " unsposSlashFmt \ ")\n", \ L, \ (*rightSeg)->b1, (*rightSeg)->b2, \ (*rightSeg)->e1, (*rightSeg)->e2); \ } #define snoopLR_reverse_left_3 \ if ((snoop) && (L > (sgnpos) LY)) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ fprintf (stderr, "update_LR: LY <- " unsposFmt \ "\n", \ (unspos) max ((sgnpos) LY, L)); \ } #define snoopLR_reverse_right_1 \ if (snoop) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ if ((*leftSeg)->type == diagSeg) \ fprintf (stderr, "update_LR: R <- " sgnposFmt \ " (on same diag seg)\n", \ R); \ else \ fprintf (stderr, "update_LR: R <- " sgnposFmt \ " (on same vert seg)\n", \ R); \ } #define snoopLR_reverse_right_2 \ if ((snoop) && (*leftSeg != NULL)) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ fprintf (stderr, "update_LR: R <- " sgnposFmt \ " (on new seg " unsposSlashFmt " " unsposSlashFmt \ ")\n", \ R, \ (*leftSeg)->b1, (*leftSeg)->b2, \ (*leftSeg)->e1, (*leftSeg)->e2); \ } #define snoopLR_reverse_right_3 \ if ((snoop) && (R < (sgnpos) RY)) \ { \ if (!havePrinted) \ { fprintf (stderr, "\n"); havePrinted = true; } \ fprintf (stderr, "update_LR: RY <- " unsposFmt \ "\n", \ special_min (RY, R)); \ } #endif // snoopAlgorithm #if ((!defined snoopAlgorithm)&&(defined snoopSubprobs)) #define snoopLR_forward_left_1 \ if (true) \ { \ if ((*leftSeg)->type == diagSeg) \ fprintf (stderr, "\n update_LR: L <- " sgnposFmt \ " (on same diag seg)", \ L); \ else \ fprintf (stderr, "\n update_LR: L <- " sgnposFmt \ " (on same vert seg)", \ L); \ } #define snoopLR_forward_left_2 \ if (*leftSeg != NULL) \ { \ fprintf (stderr, "\n update_LR: L <- " sgnposFmt \ " (on new seg " unsposSlashFmt " " unsposSlashFmt \ ")", \ L, \ (*leftSeg)->b1, (*leftSeg)->b2, \ (*leftSeg)->e1, (*leftSeg)->e2); \ } #define snoopLR_forward_left_3 \ if (L > (sgnpos) LY) \ { \ fprintf (stderr, "\n update_LR: LY <- " unsposFmt, \ (unspos) max ((sgnpos) LY, L)); \ } #define snoopLR_forward_right_1 \ if (true) \ { \ if ((*rightSeg)->type == diagSeg) \ fprintf (stderr, "\n update_LR: L <- " sgnposFmt \ " (on same diag seg)", \ R); \ else \ fprintf (stderr, "\n update_LR: R <- " sgnposFmt \ " (on same vert seg)", \ R); \ } #define snoopLR_forward_right_2 \ if (*rightSeg != NULL) \ { \ fprintf (stderr, "\n update_LR: R <- " sgnposFmt \ " (on new seg " unsposSlashFmt " " unsposSlashFmt \ ")", \ R, \ (*rightSeg)->b1, (*rightSeg)->b2, \ (*rightSeg)->e1, (*rightSeg)->e2); \ } #define snoopLR_forward_right_3 \ if (R < (sgnpos) RY) \ { \ fprintf (stderr, "\n update_LR: RY <- " unsposFmt, \ special_min (RY, R)); \ } #define snoopLR_reverse_left_1 \ if (true) \ { \ if ((*rightSeg)->type == diagSeg) \ fprintf (stderr, "\n update_LR: L <- " sgnposFmt \ " (on same diag seg)", \ L); \ else \ fprintf (stderr, "\n update_LR: L <- " sgnposFmt \ " (on same vert seg)", \ L); \ } #define snoopLR_reverse_left_2 \ if (*rightSeg != NULL) \ { \ fprintf (stderr, "\n update_LR: L <- " sgnposFmt \ " (on new seg " unsposSlashFmt " " unsposSlashFmt \ ")", \ L, \ (*rightSeg)->b1, (*rightSeg)->b2, \ (*rightSeg)->e1, (*rightSeg)->e2); \ } #define snoopLR_reverse_left_3 \ if (L > (sgnpos) LY) \ { \ fprintf (stderr, "\n update_LR: LY <- " unsposFmt, \ (unspos) max ((sgnpos) LY, L)); \ } #define snoopLR_reverse_right_1 \ if (true) \ { \ if ((*leftSeg)->type == diagSeg) \ fprintf (stderr, "\n update_LR: R <- " sgnposFmt \ " (on same diag seg)", \ R); \ else \ fprintf (stderr, "\n update_LR: R <- " sgnposFmt \ " (on same vert seg)", \ R); \ } #define snoopLR_reverse_right_2 \ if (*leftSeg != NULL) \ { \ fprintf (stderr, "\n update_LR: R <- " sgnposFmt \ " (on new seg " unsposSlashFmt " " unsposSlashFmt \ ")", \ R, \ (*leftSeg)->b1, (*leftSeg)->b2, \ (*leftSeg)->e1, (*leftSeg)->e2); \ } #define snoopLR_reverse_right_3 \ if (R < (sgnpos) RY) \ { \ fprintf (stderr, "\n update_LR: RY <- " unsposFmt, \ special_min (RY, R)); \ } #endif // not snoopAlgorithm and snoopSubprobs static unspos special_min(unspos RY, sgnpos R); static unspos special_min(unspos RY, sgnpos R) { // ideally we would return min(RY,R); but R is signed, so if R is negative // we need to return zero (since we are to return an unsigned); if R is // positive then we can return min(RY,R), but we may as well compute it // ourselves (instead of calling min); note that the cast of R to unsigned // only occurs when R is positive, so we needn't worry about it causing // overflaw (and yes, I mean overflaw not overflow) if (R <= 0) return 0; else if ((unspos) R < RY) return R; else return RY; } static void update_LR_bounds (int reversed, aliseg** rightSeg, aliseg** leftSeg, galign** rightAlign, galign** leftAlign, unspos row, unspos anchor1, unspos anchor2, sgnpos* _L, sgnpos* _R, unspos* _LY, unspos* _RY) { sgnpos L = *_L; sgnpos R = *_R; unspos LY = *_LY; unspos RY = *_RY; #ifdef snoopAlgorithm int snoop = true; int havePrinted = false; #endif // snoopAlgorithm dbg_timing_gapped_extend_sub (debugClockUpdateLrBounds); #ifdef snoopAlgorithm #if ((defined debugPosA1) && (defined debugPosB1)) snoop = ((anchor1 == debugPosA1) && (anchor2 == debugPosA2)) || ((anchor1 == debugPosB1) && (anchor2 == debugPosB2)); #elif (defined debugPosA1) snoop = (anchor1 == debugPosA1) && (anchor2 == debugPosA2); #elif (defined debugPosB1) snoop = (anchor1 == debugPosB1) && (anchor2 == debugPosB2); #endif // debugPosA1,debugPosB1 #endif // snoopAlgorithm // handle forward search (DP row advances upward) if (!reversed) { if (*leftSeg != NULL) { if ((*leftSeg)->e1 >= row + anchor1) // stay on same segment { if ((*leftSeg)->type == diagSeg) L++; snoopLR_forward_left_1; } else // move to next segment { L = next_sweep_seg (/*right*/ false, leftSeg, leftAlign, row, anchor1, anchor2) + 1; snoopLR_forward_left_2; } } if (*leftSeg != NULL) // (next_sweep_seg may have changed *leftSeg) { snoopLR_forward_left_3; LY = (unspos) max ((sgnpos) LY, L); } if (*rightSeg != NULL) { if ((*rightSeg)->e1 >= row + anchor1) // stay on same segment { if ((*rightSeg)->type == diagSeg) R++; snoopLR_forward_right_1; } else // move to next segment { R = next_sweep_seg (/*right*/ true, rightSeg, rightAlign, row, anchor1, anchor2) - 1; snoopLR_forward_right_2; } } if (*rightSeg != NULL) // (next_sweep_seg may have changed *rightSeg) { snoopLR_forward_right_3; RY = special_min (RY, R); } } // handle reversed search (DP row advances downward) else { if (*rightSeg != NULL) { if ((*rightSeg)->b1 <= anchor1 - row) // stay on same segment { if ((*rightSeg)->type == diagSeg) L++; snoopLR_reverse_left_1; } else // move to next segment { L = prev_sweep_seg (/*right*/ true, rightSeg, rightAlign, row, anchor1, anchor2) + 1; snoopLR_reverse_left_2; } } if (*rightSeg != NULL) // (prev_sweep_seg may have changed *rightSeg) { snoopLR_reverse_left_3; LY = (unspos) max ((sgnpos) LY, L); } if (*leftSeg != NULL) { if ((*leftSeg)->b1 <= anchor1 - row) // stay on same segment { if ((*leftSeg)->type == diagSeg) R++; snoopLR_reverse_right_1; } else // move to next segment { R = prev_sweep_seg (/*right*/ false, leftSeg, leftAlign, row, anchor1, anchor2) - 1; snoopLR_reverse_right_2; } } if (*leftSeg != NULL) // (prev_sweep_seg may have changed *leftSeg) { snoopLR_reverse_right_3; RY = special_min (RY, R); } } *_L = L; *_R = R; *_LY = LY; *_RY = RY; dbg_timing_gapped_extend_add (debugClockUpdateLrBounds); } //---------- // // next_sweep_seg-- // Move to the next leftSeg or rightSeg in a forward (upward) sweep. // prev_sweep_seg-- // Move to the previous leftSeg or rightSeg in a reverse (downward) sweep. // //---------- // // Arguments: // int lookRight: true => move to the next alignment to the right // .. when we run off the end of an alignment // false => move to next alignment to left instead // aliseg** bp: Current constraining segment; may be updated by // .. this function. // galign** mp: Alignment containing that segment; may be updated // .. by this function. // unspos row: The sweep row. // unspos anchor1: The position at which the alignment // unspos anchor2: .. began. // // Returns: // The column position of the next (or previous) leftSeg or rightSeg if there // is one; otherwise, the column position of the next (or previous) leftAlign // or rightAlign. If there is also no such alignment, zero is returned. // //---------- static sgnpos next_sweep_seg (int lookRight, aliseg** bp, galign** mp, unspos row, unspos anchor1, unspos anchor2) { sgnpos col; dbg_timing_gapped_extend_sub (debugClockNextSweepSeg); // move to the next segment, if there is one *bp = (*bp)->nextSeg; if (*bp != NULL) { if (((*bp)->type == horzSeg) && ((*bp = (*bp)->nextSeg) == NULL)) suicide ("Last alignment segment was horizontal"); col = signed_difference ((*bp)->b2, anchor2); dbg_timing_gapped_extend_add (debugClockNextSweepSeg); return col; } // we've run off the end (top) of an alignment; move to the next one if (lookRight) { *bp = (*mp)->rightSeg2; *mp = (*mp)->rightAlign2; } else { *bp = (*mp)->leftSeg2; *mp = (*mp)->leftAlign2; } if (*bp == NULL) { dbg_timing_gapped_extend_add (debugClockNextSweepSeg); return 0; // no constraint; there was no "next alignment" } // figure out the column where the start of the new segment intersects the // line y=row if ((*bp)->type == diagSeg) col = (sgnpos) row + signed_difference ((*bp)->b2, anchor2) - signed_difference ((*bp)->b1, anchor1); else // we jumped to a vertical segment col = signed_difference ((*bp)->b2, anchor2); dbg_timing_gapped_extend_add (debugClockNextSweepSeg); return col; } static sgnpos prev_sweep_seg (int lookRight, aliseg** bp, galign** mp, unspos row, unspos anchor1, unspos anchor2) { sgnpos col; dbg_timing_gapped_extend_sub (debugClockPrevSweepSeg); // move to the previous segment, if there is one *bp = (*bp)->prevSeg; if (*bp != NULL) { if (((*bp)->type == horzSeg) && ((*bp = (*bp)->prevSeg) == NULL)) suicide ("First alignment segment was horizontal"); col = signed_difference (anchor2, (*bp)->e2); dbg_timing_gapped_extend_add (debugClockPrevSweepSeg); return col; } // we've run off the front (bottom) of an alignment; move to the previous // one if (lookRight) { *bp = (*mp)->rightSeg1; *mp = (*mp)->rightAlign1; } else { *bp = (*mp)->leftSeg1; *mp = (*mp)->leftAlign1; } if (*bp == NULL) { dbg_timing_gapped_extend_add (debugClockPrevSweepSeg); return 0; // no constraint; there was no "previous alignment" } // figure out the column where the end of the new segment intersects the // line y=row if ((*bp)->type == diagSeg) col = (sgnpos) row + signed_difference (anchor2, (*bp)->e2) - signed_difference (anchor1, (*bp)->e1); else // we jumped to a vertical segment col = signed_difference (anchor2, (*bp)->e2); dbg_timing_gapped_extend_add (debugClockPrevSweepSeg); return col; } //---------- // // update_active_segs-- // As we move to a new sweep row, update the list of segments that intersect // the sweep row within the feasible region. // //---------- // // Arguments: // int reversed: true => the DP row advances downward // false => it advances upward // activeseg** active: The list of active segments. This may be updated // .. by this function. // galign** alignList: Alignments in advance of the sweep row. This may // .. be updated by this function. // dpCell* dp: First DP cell in (conceptual) sweep row. This is // .. indexed from LY to RY, inclusive. // unspos row: The sweep row. // unspos anchor1: The position at which the alignment // unspos anchor2: .. began. // unspos LY, RY: Current Y-drop constraints. // // Returns: // (nothing) // //---------- static aliseg* next_seg (aliseg* bp, int reversed) { return (reversed)? bp->prevSeg : bp->nextSeg; } static void update_active_segs (int reversed, activeseg** _active, galign** _alignList, dpCell* dp, unspos row, unspos anchor1, unspos anchor2, unspos LY, unspos RY) { activeseg* active = *_active; galign* alignList = *_alignList; activeseg* act; dbg_timing_gapped_extend_sub (debugClockUpdateActiveSegs); // process currently active segments (those that intersect the sweep row) for (act=active ; act!=NULL ; act=act->next) { if (act->type == horzSeg) suicide ("Impossible horizontal segment."); if (act->lastRow >= row) { // sweep row still intersects this segment if (act->type == diagSeg) act->x++; if ((act->x >= LY) && (act->x <= RY)) dp[act->x].mask = row; } else if ((act->seg = next_seg(act->seg, reversed)) != NULL) { // sweep row intersects the next segment of this alignment; // move to the next segment and mask its intial DP cells build_active_seg (reversed, act, dp, row, anchor1, anchor2, LY, RY); if (act->type == horzSeg) { act->seg = next_seg (act->seg, reversed); build_active_seg (reversed, act, dp, row, anchor1, anchor2, LY, RY); } } else { // sweep row has passed the end of this alignment act->filter = 1; // (mark it for deletion) } } // add any other alignments the sweep row now intersects, adding the // first segment (from the appropriate end) to the active list and // masking its intial DP cells if (!reversed) { while ((alignList!=NULL) && (alignList->pos1 - anchor1 == row)) { active = add_new_active (reversed, active, alignList, dp, row, anchor1, anchor2, LY, RY); alignList = alignList->next; } } else { while ((alignList != NULL) && (anchor1 - alignList->end1 == row)) { active = add_new_active (reversed, active, alignList, dp, row, anchor1, anchor2, LY, RY); alignList = alignList->prev; } } filter_active_segs (&active, 0); // delete blocks whose filter is not 0 *_active = active; *_alignList = alignList; dbg_timing_gapped_extend_add (debugClockUpdateActiveSegs); } //---------- // // build_active_seg-- // Create an active segment record for a given segment, and mask any DP // cells that it intersects. // //---------- // // Arguments: // int reversed: true => the DP row advances downward // false => it advances upward // activeseg* act: The active segment record, which already contains // .. the proper segment, but nothing else. // dpCell* dp: First DP cell in (conceptual) sweep row. This is // .. indexed from LY to RY, inclusive. // unspos row: The sweep row. // unspos anchor1: The position at which the alignment // unspos anchor2: .. began. // unspos LY, RY: Current Y-drop constraints. // // Returns: // (nothing) // //---------- static void build_active_seg (int reversed, activeseg* act, dpCell* dp, unspos row, unspos anchor1, unspos anchor2, unspos LY, unspos RY) { unspos horzEnd, iMin, iMax, i; act->type = act->seg->type; // nota bene: the following assigns to act->x and act->lastRow always // result in non-neagtive values if (!reversed) { act->x = act->seg->b2 - anchor2; act->lastRow = act->seg->e1 - anchor1; } else { act->x = anchor2 - act->seg->e2; act->lastRow = anchor1 - act->seg->b1; } if (act->type != horzSeg) { if ((act->x >= LY) && (act->x <= RY)) dp[act->x].mask = row; } else { horzEnd = (!reversed)? act->seg->e2 - anchor2 : anchor2 - act->seg->b2; iMin = max (LY, act->x); iMax = min (RY, horzEnd); for (i=iMin ; i<=iMax ; i++) dp[i].mask = row; } } //---------- // // add_new_active-- // Create a new active segment record and add it to the list, containing // the given alignment's terminal segment. // //---------- // // Arguments: // int reversed: true => the DP row advances downward // false => it advances upward // activeseg* active: The list of active segments. Upon return, the // .. caller should assign this function's return // .. value to this. // galign* alignList: Alignments in advance of the sweep row. This may // .. be updated by this function. // dpCell* dp: First DP cell in (conceptual) sweep row. This is // .. indexed from LY to RY, inclusive. // unspos row: The sweep row. // unspos anchor1: The position at which the alignment // unspos anchor2: .. began. // unspos LY, RY: Current Y-drop constraints. // // Returns: // Pointer to the head of the active segment list, which may be the newly // created node. // //---------- static activeseg* add_new_active (int reversed, activeseg* active, galign* alignList, dpCell* dp, unspos row, unspos anchor1, unspos anchor2, unspos LY, unspos RY) { activeseg* act = malloc_or_die ("add_new_active", sizeof(activeseg)); act->filter = 0; if (!reversed) act->seg = alignList->firstSeg; else act->seg = alignList->lastSeg; act->next = active; build_active_seg (reversed, act, dp, row, anchor1, anchor2, LY, RY); return act; } //---------- // // filter_active_segs-- // Remove (dispose of) active segments with filter values NOT equal to some // specified value. // //---------- // // Arguments: // activeseg** active: List of segments. // int filter: The filter value of segments to KEEP. The active // .. records for all other segments are removed from the // .. list and disposed of. // // Returns: // (nothing) // //---------- static void filter_active_segs (activeseg** active, int filter) { activeseg* prevAct, *act; dbg_timing_gapped_extend_sub (debugClockFilterActiveSegs); for (prevAct=NULL,act=(*active); act!=NULL ; ) { if (act->filter == filter) { prevAct = act; act = act->next; } else if (prevAct != NULL) { prevAct->next = act->next; free_if_valid ("filter_active_segs", act); act = prevAct->next; } else { *active = act->next; free_if_valid ("filter_active_segs", act); act = *active; } } dbg_timing_gapped_extend_add (debugClockFilterActiveSegs); } //---------- // // format_alignment-- // Process the edit script for a newly computed alignment by (1) storing a // linked-list form with the seeding MSP and (2) augmenting the edit script // into a form suitable for returning to the calling program. // //---------- // // Arguments: // alignio* io: The alignment to format. // galign* m: The MSP that was used to anchor the alignment. // // Returns: // Pointer to the newly allocated alignment description. // //---------- static alignel* format_alignment (alignio* io, galign* m) { unspos beg1, end1, beg2, end2; unspos height, width, i, j, startI, startJ, run; u32 opIx; editscript* script; u8* seq1, *seq2; alignel* a; beg1 = io->start1 + 1; end1 = io->stop1 + 1; beg2 = io->start2 + 1; end2 = io->stop2 + 1; script = io->script; seq1 = io->seq1; seq2 = io->seq2; height = end1 - beg1 + 1; width = end2 - beg2 + 1; opIx = 0; for (i=j=0 ; (iscript = script; a->beg1 = beg1; a->beg2 = beg2; a->end1 = end1; a->end2 = end2; a->seq1 = seq1; a->seq2 = seq2; a->s = io->s; a->next = NULL; a->isTrivial = false; a->hspId = m->hspId; return a; } //---------- // // save_seg-- // Add a gap-free segment to a seeding MSP, inserting a vertical or horizontal // piece before it if appropriate. // //---------- // // Arguments: // galign* m: The MSP to add the segment to. // unspos b1, b2: The segment's starting position. // unspos e1, e2: The segment's ending position. // // Returns: // (nothing) // //---------- static void insert_seg_to_tail (galign* mp, aliseg* bp); static void save_seg (galign* m, unspos b1, unspos b2, unspos e1, unspos e2) { aliseg* bp = malloc_or_die ("save_seg bp", sizeof(aliseg)); aliseg* bq; bp->b1 = b1; bp->b2 = b2; bp->e1 = e1; bp->e2 = e2; bp->type = diagSeg; // if the alignment is empty, create it with this as the first segment if (m->firstSeg == NULL) { m->firstSeg = bp->prevSeg = bp->nextSeg = bp; return; } // otherwise, we insert it at the tail, with a preceding vertical or // horizontal segment; we assume the previous tail was a diagonal segment // (since they are the only type ever added to the tail); further, we // assume that the previous tail ends on either the y=b1-1 or x=b2-1 line // (but not both) bq = malloc_or_die ("save_seg bq", sizeof(aliseg)); bq->type = ((b1 == m->firstSeg->prevSeg->e1+1)? horzSeg : vertSeg); bq->b1 = m->firstSeg->prevSeg->e1 + 1; bq->b2 = m->firstSeg->prevSeg->e2 + 1; bq->e1 = b1 - 1; bq->e2 = b2 - 1; insert_seg_to_tail (m, bq); insert_seg_to_tail (m, bp); } static void insert_seg_to_tail (galign* mp, aliseg* bp) { bp->prevSeg = mp->firstSeg->prevSeg; bp->nextSeg = mp->firstSeg; mp->firstSeg->prevSeg->nextSeg = bp; mp->firstSeg->prevSeg = bp; } //---------- // [[-- a seed hit reporter function --]] // // gappily_extend_hsps-- // Perform a gapped extension of a seed hit or HSP. // // Arguments and Return value: (see seed_search.h) // //---------- u32 gappily_extend_hsps (void* _info, unspos pos1, unspos pos2, unspos length, score s) { hitrepgappily* info = (hitrepgappily*) _info; seq* seq1 = info->seq1; seq* seq2 = info->seq2; seqpartition* sp1 = &seq1->partition; seqpartition* sp2 = &seq2->partition; partition* p1, *p2; unspos peak; alignio io; galign mp; aliseg* bp, *bq; u32 hIx, h; u32 returnVal; //fprintf (stderr, "working on segment " unsposSlashFmt " " unsposFmt "\n", // pos1-length, pos2-length, length); if (gapped_extend_dbgShowHsps) { fprintf (stderr, "\n"); dump_aligned_nucleotides (stderr, seq1, pos1-length, seq2, pos2-length, length); } // (move pos1/pos2 from end of segment to start) pos1 -= length; pos2 -= length; // $$$ move this test to set_up_hit_processor() if (info->scoreThresh.t != 'S') suicidef ("gappily_extend_hsps can't handle score threshold %s", score_thresh_to_string (&info->scoreThresh)); // reduce the HSP to a single point peak = segment_peak (seq1->v+pos1, seq2->v+pos2, length, info->scoring); pos1 += peak; pos2 += peak; //length = 0; // (unnecessary) #ifdef debugHspImmediate fprintf (stderr, "hsp: " unsposSlashFmt " -> " unsposSlashFmt "\n", pos1-peak, pos2-peak, pos1, pos2); #endif if (gapped_extend_dbgShowAnchors) { segtable st; segment* seg = &st.seg[0]; st.size = 1; st.len = 1; st.haveScores = true; st.coverageLimit = 0; st.coverage = length; st.lowScore = s; seg->pos1 = pos1; seg->pos2 = pos2; seg->length = 0; seg->s = s; seg->id = seq2->revCompFlags; seg->scoreCov = length; seg->filter = false; write_segments (stderr, &st, seq1, seq2, false, gapped_extend_dbgShowAnchorsHowOften); } // build the alignio record for ydrop_align() io.seq1 = seq1->v; io.seq2 = seq2->v; io.rev1 = info->rev1; io.rev2 = info->rev2; io.low1 = 0; io.len1 = io.high1 = seq1->len; io.low2 = 0; io.len2 = io.high2 = seq2->len; io.scoring = info->scoring; io.yDrop = info->yDrop; io.trimToPeak = info->trimToPeak; io.anchor1 = pos1; io.anchor2 = pos2; if (sp1->p != NULL) { p1 = lookup_partition (seq1, io.anchor1); io.low1 = p1->sepBefore + 1; io.high1 = p1->sepAfter; } if (sp2->p != NULL) { p2 = lookup_partition (seq2, io.anchor2); io.low2 = p2->sepBefore + 1; io.high2 = p2->sepAfter; } if (info->traceback == NULL) suicide ("gappily_extend_hsps was given a NULL traceback pointer."); io.tb = info->traceback; // perform alignment io.leftAlign = NULL; // (convince ydrop_align() that there are no io.rightAlign = NULL; // .. neighboring/bounding alignments io.leftSeg = NULL; // .. to worry about) io.rightSeg = NULL; io.aboveList = NULL; io.belowList = NULL; ydrop_align (&io); // (find the gapped alignment) #ifdef debugHspImmediate fprintf (stderr, " gappily: " unsposSlashFmt " " unsposSlashFmt " " scoreFmt "\n", io.start1, io.start2, io.stop1, io.stop2, io.s); #endif if (io.s < info->scoreThresh.s) { #ifdef debugHspImmediate fprintf (stderr, " gappily: (fails score thresh, " scoreFmt "<" scoreFmt ")\n", io.s, info->scoreThresh.s); #endif free_if_valid ("gappily_extend_hsps io.script", io.script); mp.firstSeg = NULL; mp.align = NULL; goto return_zero; // (the alignment score is too low) } mp.firstSeg = NULL; // (convert the alignment to a linked list) mp.align = format_alignment (&io, &mp); mp.pos1 = io.start1; mp.pos2 = io.start2; mp.end1 = io.stop1; mp.end2 = io.stop2; if (mp.firstSeg == NULL) goto return_zero; // (the alignment is empty) mp.lastSeg = mp.firstSeg->prevSeg; // (record the alignment's tail and mp.firstSeg->prevSeg // .. detach the circular pointer) = mp.lastSeg->nextSeg = NULL; // subject the alignment to the user's guantlet of filtering settings if ((info->minIdentity > 0) || (info->maxIdentity < 1)) { mp.align = filter_aligns_by_identity (seq1, seq2, mp.align, info->minIdentity, info->maxIdentity); if (mp.align == NULL) goto return_zero; } if ((info->minCoverage > 0) || (info->maxCoverage < 1)) { mp.align = filter_aligns_by_coverage (seq1, seq2, mp.align, info->minCoverage, info->maxCoverage); if (mp.align == NULL) goto return_zero; } if ((info->minContinuity > 0) || (info->maxContinuity < 1)) { mp.align = filter_aligns_by_continuity (mp.align, info->minContinuity, info->maxContinuity); if (mp.align == NULL) goto return_zero; } if (info->minMatchCount > 0) { mp.align = filter_aligns_by_match_count (seq1, seq2, mp.align, info->minMatchCount); if (mp.align == NULL) goto return_zero; } if (info->maxMismatchCount >= 0) { mp.align = filter_aligns_by_mismatch_count (seq1, seq2, mp.align, info->maxMismatchCount); if (mp.align == NULL) goto return_zero; } if (info->maxSeparateGapsCount >= 0) { mp.align = filter_aligns_by_num_gaps (mp.align, info->maxSeparateGapsCount); if (mp.align == NULL) goto return_zero; } if (info->maxGapColumnsCount >= 0) { mp.align = filter_aligns_by_num_gap_columns (mp.align, info->maxGapColumnsCount); if (mp.align == NULL) goto return_zero; } // if we're to prevent duplicates, compute the alignment's hash and make // sure we haven't already seen it if (info->alignmentHashes != NULL) { // if the list is already *past* full, we don't need to check whether // this alignment is a duplicate, we can just reject it if (info->alignmentHashesSeen > info->alignmentHashesSize) goto return_zero; // compute the hash value and check if we've seen it before; note that // if the list is *just* full, we still need to check, and only reject // if it is truly a duplicate; this allows the calling routine to // properly detect when the number of alignments exceeds the limit h = alignment_hash (mp.align->beg1, mp.align->end1, seq1->revCompFlags, mp.align->beg2, mp.align->end2, seq2->revCompFlags, /*script*/ NULL); for (hIx=0 ; hIxalignmentHashesSeen ; hIx++) { if (hIx >= info->alignmentHashesSize) break; if (info->alignmentHashes[hIx] == h) goto return_zero; } hIx = info->alignmentHashesSeen++; if (hIx >= info->alignmentHashesSize) goto return_one; info->alignmentHashes[hIx] = h; } // the alignment has satisfied all the user's criteria-- print it if (info->deGapifyOutput) print_align_list_segments (mp.align); else print_align_list (mp.align); return_one: free_align_list (mp.align); mp.align = NULL; returnVal = 1; goto cleanup; return_zero: if (mp.align != NULL) { free_align_list (mp.align); mp.align = NULL; } returnVal = 0; goto cleanup; cleanup: if ((mp.align != NULL) && (mp.align->script != NULL)) free_if_valid ("gappily_extend_hsps mp.script", mp.align->script); for (bp=mp.firstSeg ; bp!=NULL ; bp=bq) { bq = bp->nextSeg; free_if_valid ("gappily_extend_hsps seg", bp); } return returnVal; } //---------- // // dump_alignio_input, dump_alignio_output-- // Dump the contents of ydrop_align's io record. // //---------- #ifdef snoopAlignioInput static void dump_alignio_input (FILE* f, alignio* io) { fprintf (f, "=====\n"); fprintf (f, "seq1 = %8p\n", io->seq1); fprintf (f, "rev1 = %8p\n", io->rev1); fprintf (f, "len1 = " unsposFmt "\n", io->len1); fprintf (f, "low1 = " unsposFmt "\n", io->low1); fprintf (f, "high1 = " unsposFmt "\n", io->high1); fprintf (f, "anchor1 = " unsposFmt "\n", io->anchor1); fprintf (f, "seq2 = %8p\n", io->seq2); fprintf (f, "rev2 = %8p\n", io->rev2); fprintf (f, "len2 = " unsposFmt "\n", io->len2); fprintf (f, "low2 = " unsposFmt "\n", io->low2); fprintf (f, "high2 = " unsposFmt "\n", io->high2); fprintf (f, "anchor2 = " unsposFmt "\n", io->anchor2); fprintf (f, "scoring = %8p\n", io->scoring); fprintf (f, "yDrop = " scoreFmt "\n", io->yDrop); fprintf (f, "trim = %s\n", (io->trimToPeak)? "yes" : "no"); fprintf (f, "tb = %8p\n", io->tb); fprintf (f, "leftAlign = %8p\n", io->leftAlign); fprintf (f, "rightAlign = %8p\n", io->rightAlign); fprintf (f, "leftSeg = %8p\n", io->leftSeg); fprintf (f, "rightSeg = %8p\n", io->rightSeg); fprintf (f, "aboveList = %8p\n", io->aboveList); fprintf (f, "belowList = %8p\n", io->belowList); } #endif // snoopAlignioInput #if ((defined(snoopAlignioOutput)) || (defined(snoopEditScripts))) static void dump_alignio_output (FILE* f, alignio* io) { fprintf (f, "s(core) = " scoreFmt "\n", io->s); fprintf (f, " start1 = " unsposFmt "\n", io->start1); fprintf (f, " stop1 = " unsposFmt "\n", io->stop1); fprintf (f, " start2 = " unsposFmt "\n", io->start2); fprintf (f, " stop2 = " unsposFmt "\n", io->stop2); // fprintf (f, " script = %p\n", io->script); } #endif // snoopAlignioOutput OR snoopEditScripts //---------- // // score_alignment-- // Determine the score of gapped alignment in two subsequences. // // Note that this is generally used only after an alignment has been modified. // The gapped_extend() produces the same alignment score as the alignment is // found. // //---------- // // Arguments: // scoreset* scoring: The scoring scheme to use. // u8* seq1: The first sequence. // unspos pos1: The subsequence start position in seq1 (origin-0). // u8* seq2: The second sequence. // unspos pos2: The subsequence start position in seq2 (origin-0). // editscript* s: The script describing the alignment. // // Returns: // The alignment's score. // //---------- score score_alignment (scoreset* scoring, u8* seq1, unspos pos1, u8* seq2, unspos pos2, editscript* script) { u8* s1 = seq1 + pos1; u8* s2 = seq2 + pos2; u8* stop; u32 opIx; editop op; u32 rpt; score similarity = 0; for (opIx=0 ; opIxlen ; opIx++) { // score s = similarity; op = script->op[opIx]; rpt = edit_op_repeat(op); if (rpt == 0) continue; op = edit_op_operation(op); switch (op) { case editopSub: stop = s1 + rpt; while (s1 < stop) similarity += scoring->sub[*(s1++)][*(s2++)]; //fprintf (stderr, "match " unsposFmt " -> " scoreFmt "\n", rpt, similarity - s); break; case editopIns: similarity -= scoring->gapOpen + (rpt * scoring->gapExtend); s2 += rpt; //fprintf (stderr, "insert " unsposFmt " -> " scoreFmt "\n", rpt, similarity - s); break; case editopDel: similarity -= scoring->gapOpen + (rpt * scoring->gapExtend); s1 += rpt; //fprintf (stderr, "delete " unsposFmt " -> " scoreFmt "\n", rpt, similarity - s); break; } } return similarity; } //---------- // // count_paired_bases-- // Count the number of paired bases in an alignment. // //---------- // // Arguments: // galign* mp: The MSP that was used to anchor the alignment. // // Returns: // The number of paired bases. // //---------- static u64 count_paired_bases (galign* mp) { aliseg* bp; u64 pairedBases; pairedBases = 0; for (bp=mp->firstSeg ; bp!=NULL ; bp=bp->nextSeg) { if (bp->type == diagSeg) pairedBases += bp->e1+1 - bp->b1; } return pairedBases; } //---------- // // warn_for_paired_bases_limit-- // Tell the user that this query exceeded the limit for paired bases. // //---------- // // Arguments: // seq* seq2: The query sequence we were aligning. // u64 maxPairedBases: (same meaning as for gapped_extend) // int overlyPairedKeep: (same meaning as for gapped_extend) // // Returns: // (nothing) // //---------- static void warn_for_paired_bases_limit (seq* seq2, u64 maxPairedBases, int overlyPairedKeep) { static int firstReport = true; seqpartition* sp2 = &seq2->partition; char* name2; char strand; if (sp2->p != NULL) name2 = "seq2"; // (seq2 is partitioned) else if (seq2->useFullNames) name2 = seq2->header; else name2 = seq2->shortHeader; strand = ((seq2->revCompFlags & rcf_rev) == 0)? '+' : '-'; fprintf (stderr, "WARNING. Query %s (%c strand) contains more than %s paired bases.\n", name2, strand, commatize(maxPairedBases)); if (firstReport) { if (overlyPairedKeep) fprintf (stderr, "Any gapped alignments already found for this query/strand are reported but the\n" "query/strand is not processed further.\n"); else fprintf (stderr, "All gapped alignments for this query/strand are discarded and the query/strand\n" "is not processed further.\n"); firstReport = false; } } //---------- // // gapped_extend_zero_stats-- // Clear the statistics for this module. // //---------- // // Arguments: // (none) // // Returns: // (nothing) // //---------- void gapped_extend_zero_stats (void) { dbg_timing_set_stat (numExtensions, 0); #ifdef collect_stats // set 'em en masse to zero memset (&gappedExtendStats, 0, sizeof(gappedExtendStats)); // set any values that might be floating point to zero (fp bit pattern for // zero may not be all-bits-zero) gapped_extend_set_stat (totalPeakScore, 0); #endif // collect_stats } //---------- // // gapped_extend_show_stats-- // Show the statistics that have been collected for this module. // //---------- // // Arguments: // FILE* f: The file to print the stats to. // // Returns: // (nothing) // //---------- void gapped_extend_show_stats (arg_dont_complain(FILE* f)) { dbg_timing_report_stat (numExtensions, "gapped extensions"); #ifdef collect_stats if (f == NULL) return; fprintf (f, " number of anchors: %s\n", commatize(gappedExtendStats.numAnchors)); fprintf (f, " anchors >= %2d bp: %s\n", anchorPeakLen, commatize(gappedExtendStats.numPeaks)); if (gappedExtendStats.numPeaks > 0) fprintf (f, "average peak score: %.1f\n", ((float)gappedExtendStats.totalPeakScore) / gappedExtendStats.numPeaks); fprintf (f, " anchors extended: %s\n", commatize(gappedExtendStats.numAnchorsExtended)); fprintf (f, " gapped extensions: %s\n", commatize(gappedExtendStats.numExtensions)); fprintf (f, " DP cells visited: %s\n", commatize(gappedExtendStats.dpCellsVisited)); if (gappedExtendStats.numExtensions > 0) fprintf (f, "DP cells/extension: %s\n", commatize((2*gappedExtendStats.dpCellsVisited+gappedExtendStats.numExtensions)/(2*gappedExtendStats.numExtensions))); fprintf (f, " max DP rows: %s\n", commatize(gappedExtendStats.maxDpRows)); fprintf (f, " max DP columns: %s\n", commatize(gappedExtendStats.maxDpColumns)); if (gappedExtendStats.zallocCallsA != 0) fprintf (f, " zalloc calls A: %s (%s bytes per)\n", commatize(gappedExtendStats.zallocCallsA), commatize((gappedExtendStats.zallocTotalA + gappedExtendStats.zallocCallsA/2) / gappedExtendStats.zallocCallsA)); fprintf (f, " zalloc total A: %s\n", commatize(gappedExtendStats.zallocTotalA)); if (gappedExtendStats.zallocCallsB != 0) fprintf (f, " zalloc calls B: %s (%s bytes per)\n", commatize(gappedExtendStats.zallocCallsB), commatize((gappedExtendStats.zallocTotalB + gappedExtendStats.zallocCallsB/2) / gappedExtendStats.zallocCallsB)); fprintf (f, " zalloc total B: %s\n", commatize(gappedExtendStats.zallocTotalB)); fprintf (f, "-------------------\n"); #endif // collect_stats } void gapped_extend_generic_stats (arg_dont_complain(FILE* f), arg_dont_complain(void (*func) (FILE*, const char*, ...))) { #ifdef collect_stats if (f == NULL) return; (*func) (f, "num_anchors=%" PRId64 "\n", gappedExtendStats.numAnchors); #endif // collect_stats }