UpdateTestsCatch

Author: telatin

.github/workflows/build.yml

@@ -40,7 +40,7 @@ jobs:
       uses: actions/cache@v1
       with:
         path: $HOME/htslib
-        key: ${{ runner.os }}-htslib-1.10
+        key: ${{ runner.os }}-htslib-1.12
 
     # Install Dependencies
     - name: Install dependencies (Linux)
@@ -73,10 +73,10 @@ jobs:
       with:
         version: ${{ matrix.nimversion }}
 
-    - uses: actions-rs/toolchain@v1
-      with:
-        toolchain: stable
-    - uses: actions-rs/cargo@v1
+#    - uses: actions-rs/toolchain@v1
+#      with:
+#        toolchain: stable
+#    - uses: actions-rs/cargo@v1
 
 
     # Build and Test
@@ -87,16 +87,14 @@ jobs:
         ls -l bin/
         
         bin/bamtocov --version
-        bin/bamtocounts --version
+        
         bin/bamtocov input/mini.bam
-        echo "--"
-        cat input/regions.bed
-        echo "--"
+
+        bin/bamtocounts --version
         bin/bamtocounts input/regions.bed input/mini.bam --debug --header
         # TESTS
-        # bash tests/all.sh
-        tests/bin/shpec tests/shpec/bamtocov.sh
-        tests/unit/bamtocov-base.sh
+        #bash tests/all.sh
+        
         
     - name: Build and Test Linux executable
       if: runner.os == 'Linux'
@@ -108,7 +106,8 @@ jobs:
         bin/bamtocounts --version
         bin/bamtocov input/mini.bam
         bin/bamtocounts input/regions.bed input/mini.bam --debug --header
-        echo SKIPPING bash tests/all.sh
+        #bash tests/all.sh
         tests/bin/shpec tests/shpec/bamtocov.sh
         tests/unit/bamtocov-base.sh
-
+        echo "WHERE: $HOME/htslib" 
+        cat $HOME/htslib/version.h 

bin/covtotarget

@@ -1,8 +1,10 @@
 # in htslib: cd 
 # nim c -a -d:static -d:release 
 
+
 -d:release
 --opt:speed
+
 @if static:
   passC:"-static"
   passl:"-static"

bin/gff2bed

@@ -95,7 +95,7 @@ proc alignments_count(table: var OrderedTable[string, stranded_counts], bam:Bam,
       stderr.writeLine("[alignments_count] Got chrom: ", chrom, " tot=", len(regions))
     for aln in bam.query(chrom):
       if debug:
-          stderr.writeLine("[alignments_count]  Got aln ", aln.start)       
+          stderr.writeLine("[alignments_count]  Got aln ", aln.start)    
       if not regions.contains(chrom) or regions[chrom].len == 0:
         continue
 
@@ -107,24 +107,32 @@ proc alignments_count(table: var OrderedTable[string, stranded_counts], bam:Bam,
       if (aln.flag and eflag) != 0: continue
 
       let readAsInterval = (aln.tid, pos_t(aln.start), pos_t(aln.stop), aln.flag.reverse)
-      if debug:
-          stderr.writeLine("[alignments_count]\tGot aln: ", readAsInterval)        
-      for interval in intersections(readAsInterval, regions, target_idx):
-        # Returns: genomic_interval_t[tuple[l1: T, l2: string]
+      var
+        s: string
+        c = 0
+      try:     
+        for interval in intersections(readAsInterval, regions, target_idx):
+          # Returns: genomic_interval_t[tuple[l1: T, l2: string]]
+          s = $interval
+          c = c + 1
+          #let feature : target_feature = (chrom: "", cid: interval.chrom.int, start: interval.start.int, stop: interval.stop.int, feature: interval.label.l2)
+          #if interval.label.l2 notin table:
+          #  stderr.writeLine("[alignments_count] Warning: unknown feature: ", interval.label.l2)
+          #  table[interval.label.l2] = (fwd: 0, rev: 0)
 
-        #let feature : target_feature = (chrom: "", cid: interval.chrom.int, start: interval.start.int, stop: interval.stop.int, feature: interval.label.l2)
-        if interval.label.l2 notin table:
-          stderr.writeLine("[alignments_count] Warning: unknown feature: ", interval.label.l2)
-          table[interval.label.l2] = (fwd: 0, rev: 0)
-        if debug:
-          stderr.writeLine("[alignments_count]\t\tfor ", interval.label.l2)
-        table[interval.label.l2].inc(aln.flag.reverse)
+          try:
+            table[interval.label.l2].inc(aln.flag.reverse)
+          except Exception as e:
+            stderr.writeLine("[alignments_count] Error key table: ", e.msg)
+          
+          if debug:
+            stderr.writeLine("[alignments_count]\t\tdone ", interval.label.l2)
+      except Exception:
+        # ⛔️ ERROR TODO FIXME ALOHA [index not in ...]
         if debug:
-          stderr.writeLine("[alignments_count]\t\tdone ", interval.label.l2)
-      if debug:
-        stderr.writeLine("[alignments_count]\t\tdone aln -----")
-  if debug:
-    stderr.writeLine("[alignments_count] Done")
+          stderr.writeLine("[alignments_count] intersections loop broken at chr=", chrom, " aln=", aln.qname, " last=", s, " c=", c)
+
+ 
 
 
 proc targetSort(x, y: target_feature): int =

nim.cfg

@@ -103,9 +103,39 @@ iterator intersections*[T](query: genomic_interval_t[T], target: target_t, idx:
     while not (query << intervals[i]) and i <= max_idx:
       yield intersection_both(query, intervals[i])
       i += 1
-  #else:
-  #  dbEcho("target seeking: chrom", chrom, "not in target")
-    
+
+iterator intersections2*[T](query: genomic_interval_t[T], target: target_t, idx: var target_index_t): genomic_interval_t[tuple[l1: T, l2: string]] =
+  #dbEcho("target intersection for query interval:", query, ", starting from:", idx)
+  let chrom = query.chrom
+  echo "IS2.1 ", chrom
+  if chrom in target:
+    let
+      intervals = target[chrom]
+      max_idx = len(intervals) - 1
+    # if the chrom changes, start from leftmost interval
+    if idx.chrom != chrom:
+      #dbEcho("target seeking: change chrom to:", idx)
+      idx = (chrom, 0)
+    # advance target interval until we reach the query interval
+    while intervals[idx.interval] << query and idx.interval < max_idx:
+      idx.interval += 1
+      echo "IS2.2 idx=", idx.interval
+      #dbEcho("target seeking: advance target index to:", idx.interval, "at target interval:", intervals[idx.interval])
+    # yield all intersections
+    var i = idx.interval
+    try:
+      discard intervals[i]
+    except Exception as e:
+      stderr.writeLine("INTERSECT2:", e.msg)
+    while not (query << intervals[i]) and i <= max_idx: 
+      echo "IS2.3 idx=", idx.interval
+      try:
+        let r= intersection_both(query, intervals[i])
+        yield r
+        i += 1
+      except Exception as e:
+        stderr.writeLine("-------- INTERSECT2:", e.msg)
+        
 proc intersects*[T](query: genomic_interval_t[T], target: target_t, idx: var target_index_t): bool =
   for i in intersections(query, target, idx):
     return to_bool(i)