"bufio"
"compress/gzip"
"encoding/gob"
+ "encoding/json"
"errors"
"flag"
"fmt"
+ "image/color/palette"
"io"
"net/http"
_ "net/http/pprof"
"git.arvados.org/arvados.git/sdk/go/arvados"
log "github.com/sirupsen/logrus"
+ "gonum.org/v1/plot"
+ "gonum.org/v1/plot/plotter"
+ "gonum.org/v1/plot/vg"
+ "gonum.org/v1/plot/vg/draw"
)
type importer struct {
skipOOO bool
outputTiles bool
includeNoCalls bool
+ outputStats string
encoder *gob.Encoder
}
flags.BoolVar(&cmd.skipOOO, "skip-ooo", false, "skip out-of-order tags")
flags.BoolVar(&cmd.outputTiles, "output-tiles", false, "include tile variant sequences in output file")
flags.BoolVar(&cmd.includeNoCalls, "include-no-calls", false, "treat tiles with no-calls as regular tiles")
+ flags.StringVar(&cmd.outputStats, "output-stats", "", "output stats to `file` (json)")
priority := flags.Int("priority", 500, "container request priority")
pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
loglevel := flags.String("loglevel", "info", "logging threshold (trace, debug, info, warn, error, fatal, or panic)")
fmt.Sprintf("-skip-ooo=%v", cmd.skipOOO),
fmt.Sprintf("-output-tiles=%v", cmd.outputTiles),
fmt.Sprintf("-include-no-calls=%v", cmd.includeNoCalls),
+ "-output-stats", "/mnt/output/stats.json",
"-tag-library", cmd.tagLibraryFile,
"-ref", cmd.refFile,
"-o", cmd.outputFile,
return 0
}
-func (cmd *importer) tileFasta(tilelib *tileLibrary, infile string) (tileSeq, error) {
+func (cmd *importer) tileFasta(tilelib *tileLibrary, infile string) (tileSeq, []importStats, error) {
var input io.ReadCloser
input, err := os.Open(infile)
if err != nil {
- return nil, err
+ return nil, nil, err
}
defer input.Close()
if strings.HasSuffix(infile, ".gz") {
input, err = gzip.NewReader(input)
if err != nil {
- return nil, err
+ return nil, nil, err
}
defer input.Close()
}
starttime := time.Now()
errs := make(chan error, 1)
todo := make(chan func() error, len(infiles)*2)
+ allstats := make([][]importStats, len(infiles)*2)
var encodeJobs sync.WaitGroup
- for _, infile := range infiles {
- infile := infile
+ for idx, infile := range infiles {
+ idx, infile := idx, infile
var phases sync.WaitGroup
phases.Add(2)
variants := make([][]tileVariantID, 2)
defer phases.Done()
log.Printf("%s starting", infile)
defer log.Printf("%s done", infile)
- tseqs, err := cmd.tileFasta(tilelib, infile)
+ tseqs, stats, err := cmd.tileFasta(tilelib, infile)
+ allstats[idx*2] = stats
var kept, dropped int
variants[0], kept, dropped = tseqs.Variants()
log.Printf("%s found %d unique tags plus %d repeats", infile, kept, dropped)
defer phases.Done()
log.Printf("%s starting", infile2)
defer log.Printf("%s done", infile2)
- tseqs, err := cmd.tileFasta(tilelib, infile2)
+ tseqs, stats, err := cmd.tileFasta(tilelib, infile2)
+ allstats[idx*2+1] = stats
var kept, dropped int
variants[1], kept, dropped = tseqs.Variants()
log.Printf("%s found %d unique tags plus %d repeats", infile2, kept, dropped)
+
return err
}
} else if fastaFilenameRe.MatchString(infile) {
defer phases.Done()
log.Printf("%s starting", infile)
defer log.Printf("%s done", infile)
- tseqs, err := cmd.tileFasta(tilelib, infile)
+ tseqs, stats, err := cmd.tileFasta(tilelib, infile)
+ allstats[idx*2] = stats
if err != nil {
return err
}
defer phases.Done()
log.Printf("%s phase %d starting", infile, phase+1)
defer log.Printf("%s phase %d done", infile, phase+1)
- tseqs, err := cmd.tileGVCF(tilelib, infile, phase)
+ tseqs, stats, err := cmd.tileGVCF(tilelib, infile, phase)
+ allstats[idx*2] = stats
var kept, dropped int
variants[phase], kept, dropped = tseqs.Variants()
log.Printf("%s phase %d found %d unique tags plus %d repeats", infile, phase+1, kept, dropped)
}
}
remain := len(todo) + int(atomic.LoadInt64(&running)) - 1
- ttl := time.Now().Sub(starttime) * time.Duration(remain) / time.Duration(cap(todo)-remain)
- eta := time.Now().Add(ttl)
- log.Printf("progress %d/%d, eta %v (%v)", cap(todo)-remain, cap(todo), eta, ttl)
+ if remain < cap(todo) {
+ ttl := time.Now().Sub(starttime) * time.Duration(remain) / time.Duration(cap(todo)-remain)
+ eta := time.Now().Add(ttl)
+ log.Printf("progress %d/%d, eta %v (%v)", cap(todo)-remain, cap(todo), eta, ttl)
+ }
}
}()
}
tileJobs.Wait()
encodeJobs.Wait()
+
go close(errs)
- return <-errs
+ err := <-errs
+ if err != nil {
+ return err
+ }
+
+ if cmd.outputStats != "" {
+ f, err := os.OpenFile(cmd.outputStats, os.O_CREATE|os.O_WRONLY, 0666)
+ if err != nil {
+ return err
+ }
+ var flatstats []importStats
+ for _, stats := range allstats {
+ flatstats = append(flatstats, stats...)
+ }
+ err = json.NewEncoder(f).Encode(flatstats)
+ if err != nil {
+ return err
+ }
+ }
+
+ return nil
}
-func (cmd *importer) tileGVCF(tilelib *tileLibrary, infile string, phase int) (tileseq tileSeq, err error) {
+func (cmd *importer) tileGVCF(tilelib *tileLibrary, infile string, phase int) (tileseq tileSeq, stats []importStats, err error) {
if cmd.refFile == "" {
err = errors.New("cannot import vcf: reference data (-ref) not specified")
return
return
}
defer consensus.Wait()
- tileseq, err = tilelib.TileFasta(fmt.Sprintf("%s phase %d", infile, phase+1), stdout)
+ tileseq, stats, err = tilelib.TileFasta(fmt.Sprintf("%s phase %d", infile, phase+1), stdout)
if err != nil {
return
}
}
return flat
}
+
+type importstatsplot struct{}
+
+func (cmd *importstatsplot) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
+ err := cmd.Plot(stdin, stdout)
+ if err != nil {
+ log.Errorf("%s", err)
+ return 1
+ }
+ return 0
+}
+
+func (cmd *importstatsplot) Plot(stdin io.Reader, stdout io.Writer) error {
+ var stats []importStats
+ err := json.NewDecoder(stdin).Decode(&stats)
+ if err != nil {
+ return err
+ }
+
+ p, err := plot.New()
+ if err != nil {
+ return err
+ }
+ p.Title.Text = "coverage preserved by import (excl X<0.65)"
+ p.X.Label.Text = "input base calls ÷ sequence length"
+ p.Y.Label.Text = "output base calls ÷ input base calls"
+ p.Add(plotter.NewGrid())
+
+ data := map[string]plotter.XYs{}
+ for _, stat := range stats {
+ data[stat.InputLabel] = append(data[stat.InputLabel], plotter.XY{
+ X: float64(stat.InputCoverage) / float64(stat.InputLength),
+ Y: float64(stat.TileCoverage) / float64(stat.InputCoverage),
+ })
+ }
+
+ nextInPalette := 0
+ for label, xys := range data {
+ s, err := plotter.NewScatter(xys)
+ if err != nil {
+ return err
+ }
+ s.GlyphStyle.Color = palette.Plan9[nextInPalette%len(palette.Plan9)]
+ s.GlyphStyle.Radius = vg.Millimeter / 2
+ s.GlyphStyle.Shape = draw.CrossGlyph{}
+ nextInPalette += 7
+ p.Add(s)
+ if false {
+ p.Legend.Add(label, s)
+ }
+ }
+ p.X.Min = 0.65
+ p.X.Max = 1
+
+ w, err := p.WriterTo(8*vg.Inch, 6*vg.Inch, "svg")
+ if err != nil {
+ return err
+ }
+ _, err = w.WriteTo(stdout)
+ return err
+}
projects / lightning.git / blobdiff