X-Git-Url: https://git.arvados.org/lightning.git/blobdiff_plain/4c5a84272388853c36b290f687a76df53c6c4b3d..9e753a193772550aec797c8cc55d2c75421b9b18:/slicenumpy.go diff --git a/slicenumpy.go b/slicenumpy.go index c48a5d54a9..fe1a77497e 100644 --- a/slicenumpy.go +++ b/slicenumpy.go @@ -13,15 +13,16 @@ import ( "net/http" _ "net/http/pprof" "os" + "regexp" "runtime" "sort" + "strconv" "strings" "sync/atomic" "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" ) @@ -49,6 +50,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`") + mergeOutput := flags.Bool("merge-output", false, "merge output into one matrix.npy and one matrix.annotations.csv") flags.IntVar(&cmd.threads, "threads", 16, "number of memory-hungry assembly threads") cmd.filter.Flags(flags) err = flags.Parse(args) @@ -70,8 +72,8 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s Name: "lightning slice-numpy", Client: arvados.NewClientFromEnv(), ProjectUUID: *projectUUID, - RAM: 240000000000, - VCPUs: 64, + RAM: 650000000000, + VCPUs: 96, Priority: *priority, KeepCache: 2, APIAccess: true, @@ -81,12 +83,13 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s return 1 } runner.Args = []string{"slice-numpy", "-local=true", - "-pprof", ":6060", - "-input-dir", *inputDir, - "-output-dir", "/mnt/output", - "-threads", fmt.Sprintf("%d", cmd.threads), - "-regions", *regionsFilename, - "-expand-regions", fmt.Sprintf("%d", *expandRegions), + "-pprof=:6060", + "-input-dir=" + *inputDir, + "-output-dir=/mnt/output", + "-threads=" + fmt.Sprintf("%d", cmd.threads), + "-regions=" + *regionsFilename, + "-expand-regions=" + fmt.Sprintf("%d", *expandRegions), + "-merge-output=" + fmt.Sprintf("%v", *mergeOutput), } runner.Args = append(runner.Args, cmd.filter.Args()...) var output string @@ -116,6 +119,12 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s 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 { @@ -127,7 +136,9 @@ 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 { @@ -148,6 +159,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) { @@ -211,7 +226,24 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s } throttleCPU.Wait() - log.Info("TODO: determining which tiles intersect given regions") + var mask *mask + if *regionsFilename != "" { + mask, err = makeMask(*regionsFilename, *expandRegions) + if err != nil { + return 1 + } + // Delete reftile entries for masked-out regions. + for tag, rt := range reftile { + if !mask.Check(strings.TrimPrefix(rt.seqname, "chr"), rt.pos, rt.pos+len(rt.tiledata)) { + delete(reftile, tag) + } + } + } + + var toMerge [][]int16 + if *mergeOutput { + toMerge = make([][]int16, len(infiles)) + } throttleMem := throttle{Max: cmd.threads} // TODO: estimate using mem and data size throttleNumpyMem := throttle{Max: cmd.threads/2 + 1} @@ -233,6 +265,12 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s if tv.Ref { continue } + if mask != nil && reftile[tv.Tag] == nil { + // Don't waste + // time/memory on + // masked-out tiles. + continue + } variants := seq[tv.Tag] if len(variants) == 0 { variants = make([]TileVariant, 100) @@ -244,7 +282,9 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s seq[tv.Tag] = variants } for _, cg := range ent.CompactGenomes { - cgs[cg.Name] = cg + if matchGenome.MatchString(cg.Name) { + cgs[cg.Name] = cg + } } return nil }) @@ -266,10 +306,14 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s 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]++ + idx := int(tag-tagstart) * 2 + if idx < len(cg.Variants) { + for allele := 0; allele < 2; allele++ { + v := cg.Variants[idx+allele] + if v > 0 && len(variants[v].Sequence) > 0 { + count[variants[v].Blake2b]++ + } + } } } // hash[i] will be the hash of @@ -313,19 +357,35 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s return err } annow := bufio.NewWriterSize(annof, 1<<20) - for tag, variants := range seq { + outcol := 0 + for tag := tagstart; tag < tagend; tag++ { rt, ok := reftile[tag] if !ok { - // Reference does not use any - // variant of this tile. - // TODO: log this? mention it - // in annotations? + if mask == nil { + outcol++ + } + // Excluded by specified + // regions, or reference does + // not use any variant of this + // tile. (TODO: log this? + // mention it in annotations?) + continue + } + variants, ok := seq[tag] + if !ok { + outcol++ continue } - outcol := tag - tagID(tagstart) reftilestr := strings.ToUpper(string(rt.tiledata)) remap := variantRemap[tag-tagstart] + done := make([]bool, len(variants)) for v, tv := range variants { + v := remap[v] + if done[v] { + continue + } else { + done[v] = true + } if len(tv.Sequence) < taglen || !bytes.HasSuffix(rt.tiledata, tv.Sequence[len(tv.Sequence)-taglen:]) { continue } @@ -335,9 +395,10 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s 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, remap[v], rt.seqname, diff.String()) + fmt.Fprintf(annow, "%d,%d,%d,%s:g.%s,%s,%d,%s,%s,%s\n", tag, outcol, v, rt.seqname, diff.String(), rt.seqname, diff.Position, diff.Ref, diff.New, diff.Left) } } + outcol++ } err = annow.Flush() if err != nil { @@ -348,49 +409,39 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s return err } - throttleNumpyMem.Acquire() log.Infof("%04d: preparing numpy", infileIdx) + throttleNumpyMem.Acquire() rows := len(cgnames) - cols := 2 * int(tagend-tagstart) + cols := 2 * outcol out := make([]int16, rows*cols) for row, name := range cgnames { out := out[row*cols:] + outcol := 0 for col, v := range cgs[name].Variants { - 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]) + tag := tagstart + tagID(col/2) + if mask != nil && reftile[tag] == nil { + continue + } + if variants, ok := seq[tag]; ok && len(variants) > int(v) && len(variants[v].Sequence) > 0 { + out[outcol] = int16(variantRemap[tag-tagstart][v]) } else { - out[col] = -1 + out[outcol] = -1 } + outcol++ } } seq = nil throttleNumpyMem.Release() - fnm := fmt.Sprintf("%s/matrix.%04d.npy", *outputDir, infileIdx) - output, err := os.Create(fnm) - if err != nil { - return err - } - defer output.Close() - bufw := bufio.NewWriterSize(output, 1<<26) - npw, err := gonpy.NewWriter(nopCloser{bufw}) - if err != nil { - return err - } - log.WithFields(logrus.Fields{ - "filename": fnm, - "rows": rows, - "cols": cols, - }).Infof("%04d: writing numpy", infileIdx) - npw.Shape = []int{rows, cols} - npw.WriteInt16(out) - err = bufw.Flush() - if err != nil { - return err - } - err = output.Close() - if err != nil { - return err + if *mergeOutput { + log.Infof("%04d: matrix fragment %d rows x %d cols", infileIdx, rows, cols) + toMerge[infileIdx] = out + } else { + fnm := fmt.Sprintf("%s/matrix.%04d.npy", *outputDir, infileIdx) + err = writeNumpyInt16(fnm, out, rows, cols) + if err != nil { + return err + } } log.Infof("%s: done (%d/%d)", infile, int(atomic.AddInt64(&done, 1)), len(infiles)) return nil @@ -399,5 +450,88 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s if err = throttleMem.Wait(); err != nil { return 1 } + if *mergeOutput { + log.Info("merging output matrix and annotations") + + annoFilename := fmt.Sprintf("%s/matrix.annotations.csv", *outputDir) + annof, err := os.Create(annoFilename) + if err != nil { + return 1 + } + annow := bufio.NewWriterSize(annof, 1<<20) + + rows := len(cgnames) + cols := 0 + for _, chunk := range toMerge { + cols += len(chunk) / rows + } + out := make([]int16, rows*cols) + startcol := 0 + for outIdx, chunk := range toMerge { + chunkcols := len(chunk) / rows + for row := 0; row < rows; row++ { + copy(out[row*cols+startcol:], chunk[row*chunkcols:(row+1)*chunkcols]) + } + toMerge[outIdx] = nil + + annotationsFilename := fmt.Sprintf("%s/matrix.%04d.annotations.csv", *outputDir, outIdx) + log.Infof("reading %s", annotationsFilename) + buf, err := os.ReadFile(annotationsFilename) + if err != nil { + return 1 + } + err = os.Remove(annotationsFilename) + if err != nil { + return 1 + } + for _, line := range bytes.Split(buf, []byte{'\n'}) { + if len(line) == 0 { + continue + } + fields := bytes.SplitN(line, []byte{','}, 3) + incol, _ := strconv.Atoi(string(fields[1])) + fmt.Fprintf(annow, "%s,%d,%s\n", fields[0], incol+startcol/2, fields[2]) + } + + startcol += chunkcols + } + err = annow.Flush() + if err != nil { + return 1 + } + err = annof.Close() + if err != nil { + return 1 + } + err = writeNumpyInt16(fmt.Sprintf("%s/matrix.npy", *outputDir), out, rows, cols) + if err != nil { + return 1 + } + } return 0 } + +func writeNumpyInt16(fnm string, out []int16, rows, cols int) error { + output, err := os.Create(fnm) + if err != nil { + return err + } + defer output.Close() + bufw := bufio.NewWriterSize(output, 1<<26) + npw, err := gonpy.NewWriter(nopCloser{bufw}) + if err != nil { + return err + } + log.WithFields(log.Fields{ + "filename": fnm, + "rows": rows, + "cols": cols, + }).Infof("writing numpy: %s", fnm) + npw.Shape = []int{rows, cols} + npw.WriteInt16(out) + err = bufw.Flush() + if err != nil { + return err + } + return output.Close() +}