"pca-py": &pythonPCA{},
"plot": &pythonPlot{},
"diff-fasta": &diffFasta{},
- "stats": &stats{},
+ "stats": &statscmd{},
"merge": &merger{},
})
)
"bufio"
"compress/gzip"
"encoding/gob"
+ "encoding/json"
"errors"
"flag"
"fmt"
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)
}
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
}
wg.Add(1)
go func() {
defer wg.Done()
- code := (&stats{}).RunCommand("lightning stats", []string{"-local"}, statsin, statsout, os.Stderr)
+ code := (&statscmd{}).RunCommand("lightning stats", []string{"-local"}, statsin, statsout, os.Stderr)
c.Check(code, check.Equals, 0)
}()
wg.Wait()
c.Logf("len(merged) %d", merged.Len())
statsout := &bytes.Buffer{}
- code = (&stats{}).RunCommand("lightning stats", []string{"-local"}, bytes.NewReader(merged.Bytes()), statsout, os.Stderr)
+ code = (&statscmd{}).RunCommand("lightning stats", []string{"-local"}, bytes.NewReader(merged.Bytes()), statsout, os.Stderr)
c.Check(code, check.Equals, 0)
c.Check(statsout.Len() > 0, check.Equals, true)
c.Logf("%s", statsout.String())
log "github.com/sirupsen/logrus"
)
-type stats struct {
+type statscmd struct {
debugUnplaced bool
}
-func (cmd *stats) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
+func (cmd *statscmd) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
var err error
defer func() {
if err != nil {
return 0
}
-func (cmd *stats) doStats(input io.Reader, output io.Writer) error {
+func (cmd *statscmd) doStats(input io.Reader, output io.Writer) error {
var ret struct {
Genomes int
CalledBases []int64
return nil
}
-func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader) (tileSeq, error) {
+type importStats struct {
+ InputFile string
+ InputLabel string
+ InputLength int
+ InputCoverage int
+ TileCoverage int
+ PathLength int
+ DroppedOutOfOrderTiles int
+}
+
+func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader) (tileSeq, []importStats, error) {
ret := tileSeq{}
type jobT struct {
label string
totalFoundTags := 0
totalPathLen := 0
skippedSequences := 0
+ stats := make([]importStats, 0, len(todo))
for job := range todo {
if len(job.fasta) == 0 {
continue
ret[job.label] = pathcopy
log.Debugf("%s %s tiled with path len %d, skipped %d", filelabel, job.label, len(path), skipped)
- log.Infof("%s %s fasta in %d coverage in %d coverage out %d", filelabel, job.label, len(job.fasta), countBases(job.fasta), basesOut)
+ basesIn := countBases(job.fasta)
+ log.Infof("%s %s fasta in %d coverage in %d coverage out %d", filelabel, job.label, len(job.fasta), basesIn, basesOut)
+ stats = append(stats, importStats{
+ InputFile: filelabel,
+ InputLabel: job.label,
+ InputLength: len(job.fasta),
+ InputCoverage: basesIn,
+ TileCoverage: basesOut,
+ PathLength: len(path),
+ DroppedOutOfOrderTiles: skipped,
+ })
+
totalPathLen += len(path)
}
log.Printf("%s tiled with total path len %d in %d sequences (skipped %d sequences with '_' in name, skipped %d out-of-order tags)", filelabel, totalPathLen, len(ret), skippedSequences, totalFoundTags-totalPathLen)
- return ret, scanner.Err()
+ return ret, stats, scanner.Err()
}
func (tilelib *tileLibrary) Len() int {
func (s *tilelibSuite) TestSkipOOO(c *check.C) {
// tags appear in seq: 4, 0, 2 (but skipOOO is false)
tilelib := &tileLibrary{taglib: &s.taglib, skipOOO: false}
- tseq, err := tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
+ tseq, _, err := tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
s.tag[4]+
"ggggggggggggggggggggggg\n"+
s.tag[0]+
// tags appear in seq: 0, 1, 2 -> don't skip
tilelib = &tileLibrary{taglib: &s.taglib, skipOOO: true}
- tseq, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
+ tseq, _, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
s.tag[0]+
"cccccccccccccccccccc\n"+
s.tag[1]+
// tags appear in seq: 2, 3, 4 -> don't skip
tilelib = &tileLibrary{taglib: &s.taglib, skipOOO: true}
- tseq, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
+ tseq, _, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
s.tag[2]+
"cccccccccccccccccccc\n"+
s.tag[3]+
// tags appear in seq: 4, 0, 2 -> skip 4
tilelib = &tileLibrary{taglib: &s.taglib, skipOOO: true}
- tseq, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
+ tseq, _, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
s.tag[4]+
"cccccccccccccccccccc\n"+
s.tag[0]+
// tags appear in seq: 0, 2, 1 -> skip 2
tilelib = &tileLibrary{taglib: &s.taglib, skipOOO: true}
- tseq, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
+ tseq, _, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
s.tag[0]+
"cccccccccccccccccccc\n"+
s.tag[2]+
// tags appear in seq: 0, 1, 1, 2 -> skip second tag1
tilelib = &tileLibrary{taglib: &s.taglib, skipOOO: true}
- tseq, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
+ tseq, _, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
s.tag[0]+
"cccccccccccccccccccc\n"+
s.tag[1]+
// tags appear in seq: 0, 1, 3, 0, 4 -> skip second tag0
tilelib = &tileLibrary{taglib: &s.taglib, skipOOO: true}
- tseq, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
+ tseq, _, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
s.tag[0]+
"cccccccccccccccccccc\n"+
s.tag[1]+
// tags appear in seq: 0, 1, 3 -> don't skip
tilelib = &tileLibrary{taglib: &s.taglib, skipOOO: true}
- tseq, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
+ tseq, _, err = tilelib.TileFasta("test-label", bytes.NewBufferString(">test-seq\n"+
s.tag[0]+
"cccccccccccccccccccc\n"+
s.tag[1]+