X-Git-Url: https://git.arvados.org/lightning.git/blobdiff_plain/176a9ff1a681d7bf252e0c445c67b361dbe2d80d..8fc5d19688d99d1ebacb9671b85ba30bbf3c3d35:/slicenumpy.go diff --git a/slicenumpy.go b/slicenumpy.go index dad23420e8..5aa97ac5a7 100644 --- a/slicenumpy.go +++ b/slicenumpy.go @@ -13,6 +13,7 @@ import ( "net/http" _ "net/http/pprof" "os" + "regexp" "runtime" "sort" "strings" @@ -21,12 +22,13 @@ import ( "git.arvados.org/arvados.git/sdk/go/arvados" "github.com/arvados/lightning/hgvs" "github.com/kshedden/gonpy" - "github.com/sirupsen/logrus" log "github.com/sirupsen/logrus" + "golang.org/x/crypto/blake2b" ) type sliceNumpy struct { - filter filter + filter filter + threads int } func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int { @@ -47,6 +49,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s ref := flags.String("ref", "", "reference name (if blank, choose last one that appears in input)") regionsFilename := flags.String("regions", "", "only output columns/annotations that intersect regions in specified bed `file`") expandRegions := flags.Int("expand-regions", 0, "expand specified regions by `N` base pairs on each side`") + flags.IntVar(&cmd.threads, "threads", 16, "number of memory-hungry assembly threads") cmd.filter.Flags(flags) err = flags.Parse(args) if err == flag.ErrHelp { @@ -67,8 +70,8 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s Name: "lightning slice-numpy", Client: arvados.NewClientFromEnv(), ProjectUUID: *projectUUID, - RAM: 250000000000, - VCPUs: 32, + RAM: 650000000000, + VCPUs: 96, Priority: *priority, KeepCache: 2, APIAccess: true, @@ -81,6 +84,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s "-pprof", ":6060", "-input-dir", *inputDir, "-output-dir", "/mnt/output", + "-threads", fmt.Sprintf("%d", cmd.threads), "-regions", *regionsFilename, "-expand-regions", fmt.Sprintf("%d", *expandRegions), } @@ -106,10 +110,18 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s var cgnames []string var refseq map[string][]tileLibRef + var reftiledata = make(map[tileLibRef][]byte, 11000000) in0, err := open(infiles[0]) if err != nil { return 1 } + + matchGenome, err := regexp.Compile(cmd.filter.MatchGenome) + if err != nil { + err = fmt.Errorf("-match-genome: invalid regexp: %q", cmd.filter.MatchGenome) + return 1 + } + taglen := -1 DecodeLibrary(in0, strings.HasSuffix(infiles[0], ".gz"), func(ent *LibraryEntry) error { if len(ent.TagSet) > 0 { @@ -121,7 +133,14 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s } } for _, cg := range ent.CompactGenomes { - cgnames = append(cgnames, cg.Name) + if matchGenome.MatchString(cg.Name) { + cgnames = append(cgnames, cg.Name) + } + } + for _, tv := range ent.TileVariants { + if tv.Ref { + reftiledata[tileLibRef{tv.Tag, tv.Variant}] = tv.Sequence + } } return nil }) @@ -137,6 +156,10 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s err = fmt.Errorf("tagset not found") return 1 } + if len(cgnames) == 0 { + err = fmt.Errorf("no genomes found matching regexp %q", cmd.filter.MatchGenome) + return 1 + } sort.Strings(cgnames) { @@ -162,7 +185,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s } } - log.Info("building list of reference tiles to load") // TODO: more efficient if we had saved all ref tiles in slice0 + log.Info("indexing reference tiles") type reftileinfo struct { variant tileVariantID seqname string // chr1 @@ -172,74 +195,70 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s reftile := map[tagID]*reftileinfo{} for seqname, cseq := range refseq { for _, libref := range cseq { - reftile[libref.Tag] = &reftileinfo{seqname: seqname, variant: libref.Variant} - } - } - log.Info("loading reference tiles from all slices") - throttle := throttle{Max: runtime.GOMAXPROCS(0)} - for _, infile := range infiles { - infile := infile - throttle.Go(func() error { - defer log.Infof("%s: done", infile) - f, err := open(infile) - if err != nil { - return err + reftile[libref.Tag] = &reftileinfo{ + seqname: seqname, + variant: libref.Variant, + tiledata: reftiledata[libref], } - defer f.Close() - return DecodeLibrary(f, strings.HasSuffix(infile, ".gz"), func(ent *LibraryEntry) error { - for _, tv := range ent.TileVariants { - if dst, ok := reftile[tv.Tag]; ok && dst.variant == tv.Variant { - dst.tiledata = tv.Sequence - } - } - return nil - }) - }) + } } - throttle.Wait() + throttleCPU := throttle{Max: runtime.GOMAXPROCS(0)} log.Info("reconstructing reference sequences") for seqname, cseq := range refseq { seqname, cseq := seqname, cseq - throttle.Go(func() error { + throttleCPU.Go(func() error { defer log.Printf("... %s done", seqname) pos := 0 for _, libref := range cseq { rt := reftile[libref.Tag] rt.pos = pos + if len(rt.tiledata) == 0 { + return fmt.Errorf("missing tiledata for tag %d variant %d in %s in ref", libref.Tag, libref.Variant, seqname) + } pos += len(rt.tiledata) - taglen } return nil }) } - throttle.Wait() + throttleCPU.Wait() log.Info("TODO: determining which tiles intersect given regions") + throttleMem := throttle{Max: cmd.threads} // TODO: estimate using mem and data size + throttleNumpyMem := throttle{Max: cmd.threads/2 + 1} log.Info("generating annotations and numpy matrix for each slice") var done int64 for infileIdx, infile := range infiles { infileIdx, infile := infileIdx, infile - throttle.Go(func() error { - seq := map[tagID][][]byte{} + throttleMem.Go(func() error { + seq := make(map[tagID][]TileVariant, 50000) cgs := make(map[string]CompactGenome, len(cgnames)) f, err := open(infile) if err != nil { return err } defer f.Close() - log.Infof("reading %s", infile) + log.Infof("%04d: reading %s", infileIdx, infile) err = DecodeLibrary(f, strings.HasSuffix(infile, ".gz"), func(ent *LibraryEntry) error { for _, tv := range ent.TileVariants { + if tv.Ref { + continue + } variants := seq[tv.Tag] + if len(variants) == 0 { + variants = make([]TileVariant, 100) + } for len(variants) <= int(tv.Variant) { - variants = append(variants, nil) + variants = append(variants, TileVariant{}) } - variants[int(tv.Variant)] = tv.Sequence + variants[int(tv.Variant)] = tv seq[tv.Tag] = variants } for _, cg := range ent.CompactGenomes { - cgs[cg.Name] = cg + if matchGenome.MatchString(cg.Name) { + cgs[cg.Name] = cg + } } return nil }) @@ -250,10 +269,59 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s tagend := cgs[cgnames[0]].EndTag // TODO: filters - // TODO: tidy/renumber + + log.Infof("%04d: renumber/dedup variants for tags %d-%d", infileIdx, tagstart, tagend) + variantRemap := make([][]tileVariantID, tagend-tagstart) + throttleCPU := throttle{Max: runtime.GOMAXPROCS(0)} + for tag, variants := range seq { + tag, variants := tag, variants + throttleCPU.Acquire() + go func() { + defer throttleCPU.Release() + count := make(map[[blake2b.Size256]byte]int, len(variants)) + for _, cg := range cgs { + idx := (tag - tagstart) * 2 + if int(idx) < len(cg.Variants) { + count[variants[cg.Variants[idx]].Blake2b]++ + count[variants[cg.Variants[idx+1]].Blake2b]++ + } + } + // hash[i] will be the hash of + // the variant(s) that should + // be at rank i (0-based). + hash := make([][blake2b.Size256]byte, 0, len(count)) + for b := range count { + hash = append(hash, b) + } + sort.Slice(hash, func(i, j int) bool { + bi, bj := &hash[i], &hash[j] + if ci, cj := count[*bi], count[*bj]; ci != cj { + return ci > cj + } else { + return bytes.Compare((*bi)[:], (*bj)[:]) < 0 + } + }) + // rank[b] will be the 1-based + // new variant number for + // variants whose hash is b. + rank := make(map[[blake2b.Size256]byte]tileVariantID, len(hash)) + for i, h := range hash { + rank[h] = tileVariantID(i + 1) + } + // remap[v] will be the new + // variant number for original + // variant number v. + remap := make([]tileVariantID, len(variants)) + for i, tv := range variants { + remap[i] = rank[tv.Blake2b] + } + variantRemap[tag-tagstart] = remap + }() + } + throttleCPU.Wait() annotationsFilename := fmt.Sprintf("%s/matrix.%04d.annotations.csv", *outputDir, infileIdx) - log.Infof("writing %s", annotationsFilename) + log.Infof("%04d: writing %s", infileIdx, annotationsFilename) annof, err := os.Create(annotationsFilename) if err != nil { return err @@ -270,17 +338,18 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s } outcol := tag - tagID(tagstart) reftilestr := strings.ToUpper(string(rt.tiledata)) - for v, seq := range variants { - if len(seq) == 0 || !bytes.HasSuffix(rt.tiledata, seq[len(seq)-taglen:]) { + remap := variantRemap[tag-tagstart] + for v, tv := range variants { + if len(tv.Sequence) < taglen || !bytes.HasSuffix(rt.tiledata, tv.Sequence[len(tv.Sequence)-taglen:]) { continue } - if lendiff := len(rt.tiledata) - len(seq); lendiff < -1000 || lendiff > 1000 { + if lendiff := len(rt.tiledata) - len(tv.Sequence); lendiff < -1000 || lendiff > 1000 { continue } - diffs, _ := hgvs.Diff(reftilestr, strings.ToUpper(string(seq)), 0) + diffs, _ := hgvs.Diff(reftilestr, strings.ToUpper(string(tv.Sequence)), 0) for _, diff := range diffs { diff.Position += rt.pos - fmt.Fprintf(annow, "%d,%d,%d,%s:g.%s\n", tag, outcol, v, rt.seqname, diff.String()) + fmt.Fprintf(annow, "%d,%d,%d,%s:g.%s,%s,%d,%s,%s\n", tag, outcol, remap[v], rt.seqname, diff.String(), rt.seqname, diff.Position, diff.Ref, diff.New) } } } @@ -293,20 +362,23 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s return err } - log.Infof("%s: preparing numpy", infile) + throttleNumpyMem.Acquire() + log.Infof("%04d: preparing numpy", infileIdx) rows := len(cgnames) cols := 2 * int(tagend-tagstart) out := make([]int16, rows*cols) for row, name := range cgnames { out := out[row*cols:] for col, v := range cgs[name].Variants { - if variants, ok := seq[tagstart+tagID(col/2)]; ok && len(variants) > int(v) && len(variants[v]) > 0 { - out[col] = int16(v) + if variants, ok := seq[tagstart+tagID(col/2)]; ok && len(variants) > int(v) && len(variants[v].Sequence) > 0 { + out[col] = int16(variantRemap[col/2][v]) } else { out[col] = -1 } } } + seq = nil + throttleNumpyMem.Release() fnm := fmt.Sprintf("%s/matrix.%04d.npy", *outputDir, infileIdx) output, err := os.Create(fnm) @@ -319,11 +391,11 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s if err != nil { return err } - log.WithFields(logrus.Fields{ + log.WithFields(log.Fields{ "filename": fnm, "rows": rows, "cols": cols, - }).Info("writing numpy") + }).Infof("%04d: writing numpy", infileIdx) npw.Shape = []int{rows, cols} npw.WriteInt16(out) err = bufw.Flush() @@ -338,7 +410,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s return nil }) } - if err = throttle.Wait(); err != nil { + if err = throttleMem.Wait(); err != nil { return 1 } return 0