GTF support

Author: telatin

bin/covtotarget

@@ -114,20 +114,35 @@ def makeTarget(genomeSize, targetSize, numFeatures=1):
             target[start] = end
     return target
 
+def savetarget(target, file):
+    """
+    Save target in bed format
+    """
+    bed = open(file, "w")
+    chrname = "chromosome"
+    
+    for t in target:
+        name = f"chromosome:{t}-{target[t]}"
+        print(f"{chrname}\t{t}\t{target[t]}\t{name}", file=bed)
+    
+
+    
 
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser(description="Simulate a BAM file with long reads mapped against a hypothetical genome")
     parser.add_argument("-o", "--output", help="Output BAM file", required=True)
     parser.add_argument("-l", "--length", help="Length of the genome [default: %(default)s]", default="100M")
 
-    parser.add_argument("-n", "--num-reads", help="Number of reads [default: %(default)s]", default="100M")
+    parser.add_argument("-n", "--num-reads", help="Number of reads [default: %(default)s]", default="1M")
     parser.add_argument("-t", "--target-size", help="Bases in target [default: %(default)s]", default="10M")
     parser.add_argument("-f", "--target-features", help="Number of features [default: %(default)s]", type=int, default=10)
 
-    parser.add_argument("-m", "--min-len", help="Minimum read length [default: 1000]", type=int, default=1000)
-    parser.add_argument("-M", "--max-len", help="Minimum read length [default: 10000]", type=int,default=10000)
+    parser.add_argument("-m", "--min-len", help="Minimum read length [default: 1000]", type=int, default=50)
+    parser.add_argument("-M", "--max-len", help="Minimum read length [default: 10000]", type=int,default=300)
     parser.add_argument("-s", "--seed", help="Random seed [default: 42]", type=int,default=42)
+    parser.add_argument("--multiply", help="Multiply the number of reads by this number", type=int, default=1)
+    parser.add_argument("--randomcigar", help="Use random CIGAR strings", action="store_true")
     parser.add_argument("--progress", help="Print progress every INT reads [default: 10000]", type=int, default=10000)
 
 
@@ -144,14 +159,18 @@ def makeTarget(genomeSize, targetSize, numFeatures=1):
     bamHeader = makeHeader(genome)
     target = makeTarget(genomeSize, stringSizeToInt(opts.target_size), opts.target_features)
 
+    bedOutput = opts.output.replace("bam", "bed")
+    savetarget(target, bedOutput)
     totalReads = stringSizeToInt(opts.num_reads)
     generatedBases = 0
     n = 0
     seqPosLength = {}
 
-    while n < totalReads:
+    for n in track(range(totalReads), total=totalReads, description="Generating positions..."):
+    #while n < totalReads:
+        #n += 1
         readLength = int(opts.min_len + (opts.max_len - opts.min_len) * random.random())
-        n += 1
+        
         # Generate a random position being included in a random interval of the target
         intervalStart = random.choice(list(target.keys()))
         intervalEnd = target[intervalStart]
@@ -162,16 +181,15 @@ def makeTarget(genomeSize, targetSize, numFeatures=1):
         else:
             seqPosLength[pos] = [readLength]
 
-    
-    eprint(f"Generated positions for {totalReads} reads (from {len(seqPosLength)} positions)")
 
+    
+    
     # Sort seqPosLength keys ascending
     seqPosLength = {k: v for k, v in sorted(seqPosLength.items())}
-    eprint(f"Generating {totalReads} reads, {generatedBases} bp")
     # Generate reads
     with pysam.AlignmentFile(opts.output, "wb", header=bamHeader) as outf:
         n = 0
-        for index, pos in  track(enumerate(seqPosLength), total=len(seqPosLength), description="Writing BAM..."):
+        for index, pos in  track(enumerate(seqPosLength), total=len(seqPosLength), description="Writing BAM...         "):
 
 
             for length in seqPosLength[pos]:
@@ -184,9 +202,7 @@ def makeTarget(genomeSize, targetSize, numFeatures=1):
                     seqString   = "A" * length
                     try:
                         a = pysam.AlignedSegment()
-                        
-            
-                        a.query_sequence = "chr1"
+                        #a.query_sequence = "chr1"
                         a.flag = 0
                         a.reference_id = 0
                         a.reference_start = pos
@@ -208,6 +224,7 @@ def makeTarget(genomeSize, targetSize, numFeatures=1):
                             raise Exception("Read extends past end of genome")
                         outf.write(a)
 
+
                     except Exception as e:
                         eprint("ERROR:", sys.exc_info(), a)
                         eprint(f"{n} {pos} {seqPosLength[pos]}")

bin/gff2bed

@@ -210,6 +210,8 @@ template doAssert(condition: bool, message: string) = # FIXME is this already in
     stderr.writeLine("ERROR: ", message)
     quit(1)
 
+
+# COVERAGE FUNCTIONS #
 proc newCov(f = 0, r = 0): coverage_t =
   coverage_t(forward: f, reverse: r)
 
@@ -228,12 +230,24 @@ proc dec(c: var coverage_t, reverse=false) =
 proc tot(c: coverage_t): int =
   c.forward + c.reverse
 
+proc max(c1: coverage_t, c2: coverage_t): coverage_t =
+  newCov(max(c1.forward, c2.forward), max(c1.reverse, c2.reverse))
+proc min(c1: coverage_t, c2: coverage_t): coverage_t =
+  newCov(min(c1.forward, c2.forward), min(c1.reverse, c2.reverse))
+proc `+`(c1: coverage_t, c2: coverage_t): coverage_t =
+  newCov(c1.forward + c2.forward, c1.reverse + c2.reverse)
+
+proc `/`(c: coverage_t, by: float): tuple[forward: float, reverse: float] = 
+  (float(c.forward)/by, float(c.reverse)/by)
+proc `*`(c: coverage_t, by: int): coverage_t =
+  newCov(c.forward*by, c.reverse*by)
+
+
 proc topStop(q: HeapQueue): int64 =
   if not q[0].isNil:
     return q[0].stop
   return -1
 
-
 proc topReverse(q: HeapQueue): bool =
   if not q[0].isNil:
     return q[0].reverse
@@ -391,6 +405,24 @@ type
     chrom2str: TableRef[chrom_t, string]
     chrom2len: TableRef[chrom_t, pos_t]
 
+proc output_wig_span(span: genomic_interval_t[coverage_t], opts: output_option_t) =
+  let span_length = opts.span_length # FIXME the actual span can be less than span_length!
+  let value_str =
+    if opts.strand:
+      case opts.span_func:
+        of sf_max, sf_min: $span.label.forward & "\t" & $span.label.reverse
+        of sf_mean: 
+          let mean = span.label/float(span_length) 
+          $mean.forward & "\t" & $mean.reverse
+    else:
+      let tot = span.label.forward + span.label.reverse
+      case opts.span_func:
+        of sf_max, sf_min: $tot
+        of sf_mean: $(float(tot)/float(span_length))
+  echo $span.start & "\t" & value_str
+
+
+
 proc write_output(o: var output_t, i: genomic_interval_t[coverage_t]) =
   if o.current_span.chrom != i.chrom or i.start < i.stop: # skip empty intervals
     case o.opts.output_format:
@@ -412,16 +444,13 @@ proc write_output(o: var output_t, i: genomic_interval_t[coverage_t]) =
         if len(o.quantization_index2label) > 0:
           stderr.writeLine("wig output does not support quantized coverage")
           raise
-        if o.opts.strand:
-          stderr.writeLine("wig output does not support stranded coverage")
-          raise
+        #if o.opts.strand:
+        #  stderr.writeLine("wig output does not support stranded coverage")
+        #  raise
         let span_length = o.opts.span_length
         if o.current_span.chrom != i.chrom: # start new contig
           if o.current_span.chrom != -1 and o.current_span.start < o.chrom2len[o.current_span.chrom]: # output last possibly incomplete span from previous chrom
-            let span_value = case o.opts.span_func:
-              of sf_max, sf_min: $o.current_span.label.forward
-              of sf_mean: $(float(o.current_span.label.forward)/float(span_length)) # FIXME the actual span is less than span_length!
-            echo $o.current_span.start & "\t" & span_value
+            output_wig_span(o.current_span, o.opts)
           if i.chrom == -1:
             return
 
@@ -434,22 +463,18 @@ proc write_output(o: var output_t, i: genomic_interval_t[coverage_t]) =
         while o.current_span.start <= i.stop:
           let inter = intersection_first(o.current_span, i)
           if not is_empty(inter): # update the current span value
-            o.current_span.label.forward = case o.opts.span_func:
-              of sf_max: max(o.current_span.label.forward, i.label.forward)
-              of sf_min: min(o.current_span.label.forward, i.label.forward)
-              of sf_mean: o.current_span.label.forward + i.label.forward*int(len(inter))
+            o.current_span.label = case o.opts.span_func:
+              of sf_max: max(o.current_span.label, i.label)
+              of sf_min: min(o.current_span.label, i.label)
+              of sf_mean: o.current_span.label + i.label*int(len(inter))
           if inter.stop == o.current_span.stop: # span is concluded
-            # output span
-            let span_value = case o.opts.span_func:
-              of sf_max, sf_min: $o.current_span.label.forward
-              of sf_mean: $(float(o.current_span.label.forward)/float(span_length))
-            echo $o.current_span.start & "\t" & span_value
+            output_wig_span(o.current_span, o.opts)
             # next span
             o.current_span.start += span_length
             o.current_span.stop = o.current_span.start + span_length
-            o.current_span.label.forward = case o.opts.span_func
-              of sf_max, sf_mean: 0
-              of sf_min: high(int)
+            o.current_span.label = case o.opts.span_func
+              of sf_max, sf_mean: newCov(0, 0)
+              of sf_min: newCov(high(int), high(int))
           else: # span extends beyond the interval, we are done
             break
 

scripts/benchmarking/make-bam-targeted.py

@@ -42,7 +42,36 @@ proc tostring*(r: region_t, s:var string) {.inline.} =
     s.add(r.name & "\t")
   s.add($r.count)
 
+# Converts a GTF line to region object
+proc gtf_line_to_region*(line: string, gffField = "exon", gffSeparator = ";", gffIdentifier = "gene_id"): region_t =
+  #NC_001422.1     Prodigal:002006 gene    51      221     .       +       0       gene_id "nbis-gene-1"; ID "nbis-gene-1"; inference "ab initio prediction:Prodigal:002006"; locus_tag "PhiX_01"; product "hypothetical protein";
+  var
+   cse = line.strip().split('\t')
+
+  if len(cse) < 8:
+    stderr.write_line("[warning] skipping GTF line (fields not found):", line.strip())
+    return nil
+
+  # Skip non CDS fields (or user provided)
+  if cse[2] != gffField:
+    return nil
 
+  var
+    s = parse_int(cse[3])  - 1
+    e = parse_int(cse[4])
+    reg = region_t(chrom: cse[0], start: s, stop: e, count:0)
+  
+  # In the future, 8th field could be requireed [TODO]
+  if len(cse) == 9:
+    for gffAnnotPartRaw in cse[8].split(gffSeparator):
+      let gffAnnotPart = gffAnnotPartRaw.strip(chars = {'"', '\'', ' '})
+      if gffAnnotPart.startsWith(gffIdentifier):
+        try:
+          reg.name = gffAnnotPart.split("=")[1].strip(chars = {'"', '\'', ' '}) 
+        except:
+          reg.name = gffAnnotPart.split(" ")[1].strip(chars = {'"', '\'', ' '})
+        break
+  return reg
 
 # Converts a GFF line to region object
 proc gff_line_to_region*(line: string, gffField = "CDS", gffSeparator = ";", gffIdentifier = "ID"): region_t =
@@ -123,6 +152,30 @@ proc bed_to_table*(bed: string): TableRef[string, seq[region_t]] =
   return bed_regions
 
 
+proc gtf_to_table*(bed: string, gffField, gffSeparator, gffIdentifier: string): TableRef[string, seq[region_t]] =
+  var bed_regions = newTable[string, seq[region_t]]()
+  var hf = hts.hts_open(cstring(bed), "r")
+  var kstr: hts.kstring_t
+  kstr.l = 0
+  kstr.m = 0
+  kstr.s = nil
+  while hts_getline(hf, cint(10), addr kstr) > 0:
+    if ($kstr.s).startswith("##FASTA"):
+      break
+    if $kstr.s[0] == "#":
+      continue
+
+    var v = gtf_line_to_region($kstr.s, gffField, gffSeparator, gffIdentifier)
+    if v == nil: continue
+    discard bed_regions.hasKeyOrPut(v.chrom, new_seq[region_t]())
+    bed_regions[v.chrom].add(v)
+
+  # since it is read into mem, can also well sort.
+  for chrom, ivs in bed_regions.mpairs:
+    sort(ivs, proc (a, b: region_t): int = a.start - b.start)
+
+  hts.free(kstr.s)
+  return bed_regions
 
 proc gff_to_table*(bed: string, gffField, gffSeparator, gffIdentifier: string): TableRef[string, seq[region_t]] =
   var bed_regions = newTable[string, seq[region_t]]()