Faster annotate.
authorTom Clegg <tom@tomclegg.ca>
Mon, 2 Nov 2020 05:55:42 +0000 (00:55 -0500)
committerTom Clegg <tom@tomclegg.ca>
Mon, 2 Nov 2020 05:55:42 +0000 (00:55 -0500)
Arvados-DCO-1.1-Signed-off-by: Tom Clegg <tom@tomclegg.ca>

annotate.go
throttle.go [new file with mode: 0644]

index 811aefc2ecf2325bdbf16f2e392436928bfe8d38..449e64bb3449906e561cb583b0fe77caa6a8dbca 100644 (file)
@@ -25,6 +25,7 @@ import (
 type annotatecmd struct {
        variantHash bool
        maxTileSize int
+       tag2tagid   map[string]tagID
 }
 
 func (cmd *annotatecmd) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
@@ -144,9 +145,9 @@ func (cmd *annotatecmd) exportTileDiffs(outw io.Writer, tilelib *tileLibrary) er
        for name := range tilelib.refseqs {
                refs = append(refs, name)
        }
-       tag2tagid := make(map[string]tagID, len(tagset))
+       cmd.tag2tagid = make(map[string]tagID, len(tagset))
        for tagid, tagseq := range tagset {
-               tag2tagid[string(tagseq)] = tagID(tagid)
+               cmd.tag2tagid[string(tagseq)] = tagID(tagid)
        }
        sort.Strings(refs)
        log.Infof("len(refs) %d", len(refs))
@@ -163,9 +164,13 @@ func (cmd *annotatecmd) exportTileDiffs(outw io.Writer, tilelib *tileLibrary) er
        }()
        defer close(outch)
 
-       limiter := make(chan bool, runtime.NumCPU()+1)
-       var diffwg sync.WaitGroup
-       defer diffwg.Wait()
+       nseqs := 0
+       for _, refcs := range tilelib.refseqs {
+               nseqs += len(refcs)
+       }
+
+       throttle := &throttle{Max: runtime.NumCPU() + nseqs + 1}
+       defer throttle.Wait()
 
        for _, refname := range refs {
                refname := refname
@@ -177,94 +182,100 @@ func (cmd *annotatecmd) exportTileDiffs(outw io.Writer, tilelib *tileLibrary) er
                sort.Strings(seqnames)
                for _, seqname := range seqnames {
                        seqname := seqname
-                       var refseq []byte
-                       // tilestart[123] is the index into refseq
-                       // where the tile for tag 123 was placed.
-                       tilestart := map[tagID]int{}
-                       tileend := map[tagID]int{}
-                       for _, libref := range refcs[seqname] {
-                               if libref.Variant < 1 {
-                                       return fmt.Errorf("reference %q seq %q uses variant zero at tag %d", refname, seqname, libref.Tag)
-                               }
-                               seq := tilelib.TileVariantSequence(libref)
-                               if len(seq) < taglen {
-                                       return fmt.Errorf("reference %q seq %q uses tile %d variant %d with sequence len %d < taglen %d", refname, seqname, libref.Tag, libref.Variant, len(seq), taglen)
-                               }
-                               overlap := taglen
-                               if len(refseq) == 0 {
-                                       overlap = 0
-                               }
-                               tilestart[libref.Tag] = len(refseq) - overlap
-                               refseq = append(refseq, seq[overlap:]...)
-                               tileend[libref.Tag] = len(refseq)
+                       throttle.Acquire()
+                       go func() {
+                               defer throttle.Release()
+                               throttle.Report(cmd.annotateSequence(throttle, outch, tilelib, taglen, refname, seqname, refcs[seqname]))
+                       }()
+               }
+       }
+       throttle.Wait()
+       return throttle.Err()
+}
+
+func (cmd *annotatecmd) annotateSequence(throttle *throttle, outch chan<- string, tilelib *tileLibrary, taglen int, refname, seqname string, reftiles []tileLibRef) error {
+       var refseq []byte
+       // tilestart[123] is the index into refseq
+       // where the tile for tag 123 was placed.
+       tilestart := map[tagID]int{}
+       tileend := map[tagID]int{}
+       for _, libref := range reftiles {
+               if libref.Variant < 1 {
+                       return fmt.Errorf("reference %q seq %q uses variant zero at tag %d", refname, seqname, libref.Tag)
+               }
+               seq := tilelib.TileVariantSequence(libref)
+               if len(seq) < taglen {
+                       return fmt.Errorf("reference %q seq %q uses tile %d variant %d with sequence len %d < taglen %d", refname, seqname, libref.Tag, libref.Variant, len(seq), taglen)
+               }
+               overlap := taglen
+               if len(refseq) == 0 {
+                       overlap = 0
+               }
+               tilestart[libref.Tag] = len(refseq) - overlap
+               refseq = append(refseq, seq[overlap:]...)
+               tileend[libref.Tag] = len(refseq)
+       }
+       log.Infof("seq %s len(refseq) %d len(tilestart) %d", seqname, len(refseq), len(tilestart))
+       for tag, tvs := range tilelib.variant {
+               tag := tagID(tag)
+               refstart, ok := tilestart[tag]
+               if !ok {
+                       // Tag didn't place on this
+                       // reference sequence. (It
+                       // might place on the same
+                       // chromosome in a genome
+                       // anyway, but we don't output
+                       // the annotations that would
+                       // result.)
+                       continue
+               }
+               for variant := 1; variant <= len(tvs); variant++ {
+                       variant, hash := tileVariantID(variant), tvs[variant-1]
+                       tileseq := tilelib.TileVariantSequence(tileLibRef{Tag: tag, Variant: variant})
+                       if len(tileseq) < taglen {
+                               return fmt.Errorf("tilevar %d,%d has sequence len %d < taglen %d", tag, variant, len(tileseq), taglen)
                        }
-                       log.Infof("seq %s len(refseq) %d len(tilestart) %d", seqname, len(refseq), len(tilestart))
-                       for tag, tvs := range tilelib.variant {
-                               tag := tagID(tag)
-                               refstart, ok := tilestart[tag]
-                               if !ok {
-                                       // Tag didn't place on this
-                                       // reference sequence. (It
-                                       // might place on the same
-                                       // chromosome in a genome
-                                       // anyway, but we don't output
-                                       // the annotations that would
-                                       // result.)
-                                       continue
-                               }
-                               for variant := 1; variant <= len(tvs); variant++ {
-                                       variant, hash := tileVariantID(variant), tvs[variant-1]
-                                       tileseq := tilelib.TileVariantSequence(tileLibRef{Tag: tag, Variant: variant})
-                                       if len(tileseq) < taglen {
-                                               return fmt.Errorf("tilevar %d,%d has sequence len %d < taglen %d", tag, variant, len(tileseq), taglen)
-                                       }
-                                       var refpart []byte
-                                       endtag := string(tileseq[len(tileseq)-taglen:])
-                                       if endtagid, ok := tag2tagid[endtag]; !ok {
-                                               // Tile variant doesn't end on a tag, so it can only place at the end of a chromosome.
-                                               refpart = refseq[refstart:]
-                                               log.Warnf("%x tilevar %d,%d endtag not in ref: %s", hash[:13], tag, variant, endtag)
-                                       } else if refendtagstart, ok := tilestart[endtagid]; !ok {
-                                               // Ref ends a chromsome with a (possibly very large) variant of this tile, but genomes with this tile don't.
-                                               // Give up. (TODO: something smarter)
-                                               log.Debugf("%x not annotating tilevar %d,%d because end tag %d is not in ref", hash[:13], tag, variant, endtagid)
-                                               continue
+                       var refpart []byte
+                       endtag := string(tileseq[len(tileseq)-taglen:])
+                       if endtagid, ok := cmd.tag2tagid[endtag]; !ok {
+                               // Tile variant doesn't end on a tag, so it can only place at the end of a chromosome.
+                               refpart = refseq[refstart:]
+                               log.Warnf("%x tilevar %d,%d endtag not in ref: %s", hash[:13], tag, variant, endtag)
+                       } else if refendtagstart, ok := tilestart[endtagid]; !ok {
+                               // Ref ends a chromsome with a (possibly very large) variant of this tile, but genomes with this tile don't.
+                               // Give up. (TODO: something smarter)
+                               log.Debugf("%x not annotating tilevar %d,%d because end tag %d is not in ref", hash[:13], tag, variant, endtagid)
+                               continue
+                       } else {
+                               // Non-terminal tile vs. non-terminal reference.
+                               refpart = refseq[refstart : refendtagstart+taglen]
+                               log.Tracef("\n%x tilevar %d,%d endtag %s endtagid %d refendtagstart %d", hash[:13], tag, variant, endtag, endtagid, refendtagstart)
+                       }
+                       if len(refpart) > cmd.maxTileSize {
+                               log.Warnf("%x tilevar %d,%d skipping long diff, ref %s seq %s pos %d ref len %d", hash[:13], tag, variant, refname, seqname, refstart, len(refpart))
+                               continue
+                       }
+                       if len(tileseq) > cmd.maxTileSize {
+                               log.Warnf("%x tilevar %d,%d skipping long diff, ref %s seq %s variant len %d", hash[:13], tag, variant, refname, seqname, len(tileseq))
+                               continue
+                       }
+                       // log.Printf("\n%x @ refstart %d \n< %s\n> %s\n", tv.Blake2b, refstart, refpart, tileseq)
+
+                       throttle.Acquire()
+                       go func() {
+                               defer throttle.Release()
+                               diffs, _ := hgvs.Diff(strings.ToUpper(string(refpart)), strings.ToUpper(string(tileseq)), 0)
+                               for _, diff := range diffs {
+                                       diff.Position += refstart
+                                       var varid string
+                                       if cmd.variantHash {
+                                               varid = fmt.Sprintf("%x", hash)[:13]
                                        } else {
-                                               // Non-terminal tile vs. non-terminal reference.
-                                               refpart = refseq[refstart : refendtagstart+taglen]
-                                               log.Tracef("\n%x tilevar %d,%d endtag %s endtagid %d refendtagstart %d", hash[:13], tag, variant, endtag, endtagid, refendtagstart)
-                                       }
-                                       if len(refpart) > cmd.maxTileSize {
-                                               log.Warnf("%x tilevar %d,%d skipping long diff, ref %s seq %s pos %d ref len %d", hash[:13], tag, variant, refname, seqname, refstart, len(refpart))
-                                               continue
+                                               varid = fmt.Sprintf("%d", variant)
                                        }
-                                       if len(tileseq) > cmd.maxTileSize {
-                                               log.Warnf("%x tilevar %d,%d skipping long diff, ref %s seq %s variant len %d", hash[:13], tag, variant, refname, seqname, len(tileseq))
-                                               continue
-                                       }
-                                       // log.Printf("\n%x @ refstart %d \n< %s\n> %s\n", tv.Blake2b, refstart, refpart, tileseq)
-
-                                       diffwg.Add(1)
-                                       limiter <- true
-                                       go func() {
-                                               defer func() {
-                                                       <-limiter
-                                                       diffwg.Done()
-                                               }()
-                                               diffs, _ := hgvs.Diff(strings.ToUpper(string(refpart)), strings.ToUpper(string(tileseq)), 0)
-                                               for _, diff := range diffs {
-                                                       diff.Position += refstart
-                                                       var varid string
-                                                       if cmd.variantHash {
-                                                               varid = fmt.Sprintf("%x", hash)[:13]
-                                                       } else {
-                                                               varid = fmt.Sprintf("%d", variant)
-                                                       }
-                                                       outch <- fmt.Sprintf("%d\t%s\t%s\t%s:g.%s\n", tag, varid, refname, seqname, diff.String())
-                                               }
-                                       }()
+                                       outch <- fmt.Sprintf("%d\t%s\t%s\t%s:g.%s\n", tag, varid, refname, seqname, diff.String())
                                }
-                       }
+                       }()
                }
        }
        return nil
diff --git a/throttle.go b/throttle.go
new file mode 100644 (file)
index 0000000..a841100
--- /dev/null
@@ -0,0 +1,42 @@
+package main
+
+import (
+       "sync"
+       "sync/atomic"
+)
+
+type throttle struct {
+       Max       int
+       wg        sync.WaitGroup
+       ch        chan bool
+       err       atomic.Value
+       setupOnce sync.Once
+       errorOnce sync.Once
+}
+
+func (t *throttle) Acquire() {
+       t.setupOnce.Do(func() { t.ch = make(chan bool, t.Max) })
+       t.wg.Add(1)
+       t.ch <- true
+}
+
+func (t *throttle) Release() {
+       t.wg.Done()
+       <-t.ch
+}
+
+func (t *throttle) Report(err error) {
+       if err != nil {
+               t.errorOnce.Do(func() { t.err.Store(err) })
+       }
+}
+
+func (t *throttle) Err() error {
+       err, _ := t.err.Load().(error)
+       return err
+}
+
+func (t *throttle) Wait() error {
+       t.wg.Wait()
+       return t.Err()
+}