From 39afb6da60f0e8913f9c5da5193be2a12277914a Mon Sep 17 00:00:00 2001
From: Sergey Koren <sergek@umd.edu>
Date: Thu, 14 Nov 2019 17:12:31 -0500
Subject: [PATCH] Use ratio when computing the unitig positions within contigs
 too

---
 src/gfa/alignGFA.C | 22 ++++++++++++++++++----
 1 file changed, 18 insertions(+), 4 deletions(-)

diff --git a/src/gfa/alignGFA.C b/src/gfa/alignGFA.C
index 24a99ff69..6549d318b 100644
--- a/src/gfa/alignGFA.C
+++ b/src/gfa/alignGFA.C
@@ -476,19 +476,23 @@ checkRecord_align(char *label,
 bool
 checkRecord(bedRecord   *record,
             sequences   &ctgs,
+            sequences   &ctgs_orig,
             sequences   &utgs,
+            sequences   &utgs_orig,
             bool         beVerbose,
             bool         UNUSED(doPlot)) {
 
   char   *Aseq = ctgs[record->_Aid].seq;
   char   *Bseq = utgs[record->_Bid].seq, *Brev = NULL;
 
-  int32  Abgn  = record->_bgn;
-  int32  Aend  = record->_end;
-
   int32  Alen = ctgs[record->_Aid].len;
   int32  Blen = utgs[record->_Bid].len;
 
+  double ratio = max((double)Alen/ctgs_orig[record->_Aid].len, (double)Blen/utgs_orig[record->_Bid].len);
+
+  int32  Abgn  = (int32)(ceil(record->_bgn*ratio));
+  int32  Aend  = (int32)(ceil(record->_end*ratio));
+
   bool   success    = true;
   int32  alignScore = 0;
 
@@ -688,11 +692,21 @@ processBED(char   *tigName,
 
   bedFile   *bed  = new bedFile(inBED);
 
+  fprintf(stderr, "-- Loading sequences from tigStore '%s' version %u.\n", tigName, tigVers-1);
+
+  sequences *utgs_origp = new sequences(tigName, tigVers-1);
+  sequences &utgs_orig  = *utgs_origp;
+
   fprintf(stderr, "-- Loading sequences from tigStore '%s' version %u.\n", tigName, tigVers);
 
   sequences *utgsp = new sequences(tigName, tigVers);
   sequences &utgs  = *utgsp;
 
+  fprintf(stderr, "-- Loading sequences from tigStore '%s' version %u.\n", seqName, seqVers-1);
+
+  sequences *ctgs_origp = new sequences(seqName, seqVers-1);
+  sequences &ctgs_orig  = *ctgs_origp;
+
   fprintf(stderr, "-- Loading sequences from tigStore '%s' version %u.\n", seqName, seqVers);
 
   sequences *ctgsp = new sequences(seqName, seqVers);
@@ -713,7 +727,7 @@ processBED(char   *tigName,
   for (uint32 ii=0; ii<iiLimit; ii++) {
     bedRecord *record = bed->_records[ii];
 
-    if (checkRecord(record, ctgs, utgs, (verbosity > 0), false)) {
+    if (checkRecord(record, ctgs, ctgs_orig, utgs, utgs_orig, (verbosity > 0), false)) {
       pass++;
     } else {
       delete bed->_records[ii];