"net/http"
_ "net/http/pprof"
"os"
+ "regexp"
"runtime"
"sort"
"strings"
"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 {
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 {
Name: "lightning slice-numpy",
Client: arvados.NewClientFromEnv(),
ProjectUUID: *projectUUID,
- RAM: 250000000000,
- VCPUs: 32,
+ RAM: 650000000000,
+ VCPUs: 96,
Priority: *priority,
KeepCache: 2,
APIAccess: true,
"-pprof", ":6060",
"-input-dir", *inputDir,
"-output-dir", "/mnt/output",
+ "-threads", fmt.Sprintf("%d", cmd.threads),
"-regions", *regionsFilename,
"-expand-regions", fmt.Sprintf("%d", *expandRegions),
}
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 {
}
}
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
})
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)
{
}
}
- 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
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")
- throttle1 := throttle{Max: runtime.GOMAXPROCS(0)}
- for _, infile := range infiles {
- infile := infile
- throttle1.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
- })
- })
+ }
}
- throttle1.Wait()
+ throttleCPU := throttle{Max: runtime.GOMAXPROCS(0)}
log.Info("reconstructing reference sequences")
for seqname, cseq := range refseq {
seqname, cseq := seqname, cseq
- throttle1.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
})
}
- throttle1.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
- throttle1.Go(func() error {
+ throttleMem.Go(func() error {
seq := make(map[tagID][]TileVariant, 50000)
cgs := make(map[string]CompactGenome, len(cgnames))
f, err := open(infile)
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)
seq[tv.Tag] = variants
}
for _, cg := range ent.CompactGenomes {
- cgs[cg.Name] = cg
+ if matchGenome.MatchString(cg.Name) {
+ cgs[cg.Name] = cg
+ }
}
return nil
})
// TODO: filters
- log.Infof("renumber/dedup variants for tags %d-%d", tagstart, tagend)
+ log.Infof("%04d: renumber/dedup variants for tags %d-%d", infileIdx, tagstart, tagend)
variantRemap := make([][]tileVariantID, tagend-tagstart)
- throttle2 := throttle{Max: runtime.GOMAXPROCS(0)}
+ throttleCPU := throttle{Max: runtime.GOMAXPROCS(0)}
for tag, variants := range seq {
tag, variants := tag, variants
- throttle2.Acquire()
+ throttleCPU.Acquire()
go func() {
- defer throttle2.Release()
+ defer throttleCPU.Release()
count := make(map[[blake2b.Size256]byte]int, len(variants))
for _, cg := range cgs {
idx := (tag - tagstart) * 2
variantRemap[tag-tagstart] = remap
}()
}
- throttle2.Wait()
+ 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
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\n", tag, outcol, remap[v], rt.seqname, diff.String(), rt.seqname, diff.Position, diff.Ref, diff.New)
}
}
}
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)
}
}
}
+ seq = nil
+ throttleNumpyMem.Release()
fnm := fmt.Sprintf("%s/matrix.%04d.npy", *outputDir, infileIdx)
output, err := os.Create(fnm)
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()
return nil
})
}
- if err = throttle1.Wait(); err != nil {
+ if err = throttleMem.Wait(); err != nil {
return 1
}
return 0
projects / lightning.git / blobdiff