From: Tom Clegg Date: Mon, 17 Jan 2022 16:04:43 +0000 (-0500) Subject: Add -single-onehot output option. X-Git-Url: https://git.arvados.org/lightning.git/commitdiff_plain/8778237e68cdbe4fd1d4b43ed5140438fc92f3a6 Add -single-onehot output option. refs #18581 Arvados-DCO-1.1-Signed-off-by: Tom Clegg --- diff --git a/slice_test.go b/slice_test.go index 10202a4f0a..4278f13e20 100644 --- a/slice_test.go +++ b/slice_test.go @@ -244,7 +244,7 @@ pipeline1dup/input2 0 `) } - c.Log("=== slice-numpy + onehot ===") + c.Log("=== slice-numpy + onehotChunked ===") { err = ioutil.WriteFile(tmpdir+"/casecontrol.tsv", []byte(`SampleID CC pipeline1/input1 1 @@ -287,4 +287,48 @@ pipeline1dup/input2 0 }) } } + + c.Log("=== slice-numpy + onehotSingle ===") + { + err = ioutil.WriteFile(tmpdir+"/casecontrol.tsv", []byte(`SampleID CC +pipeline1/input1 1 +pipeline1/input2 0 +pipeline1dup/input1 1 +pipeline1dup/input2 0 +`), 0600) + c.Assert(err, check.IsNil) + npydir := c.MkDir() + exited := (&sliceNumpy{}).RunCommand("slice-numpy", []string{ + "-local=true", + "-single-onehot=true", + "-chi2-case-control-file=" + tmpdir + "/casecontrol.tsv", + "-chi2-case-control-column=CC", + "-chi2-p-value=0.05", + "-min-coverage=0.75", + "-input-dir=" + slicedir, + "-output-dir=" + npydir, + }, nil, os.Stderr, os.Stderr) + c.Check(exited, check.Equals, 0) + out, _ := exec.Command("find", npydir, "-ls").CombinedOutput() + c.Logf("%s", out) + + f, err := os.Open(npydir + "/onehot.npy") + c.Assert(err, check.IsNil) + defer f.Close() + npy, err := gonpy.NewReader(f) + c.Assert(err, check.IsNil) + c.Check(npy.Shape, check.DeepEquals, []int{4, 16}) + onehot, err := npy.GetInt8() + if c.Check(err, check.IsNil) { + for r := 0; r < npy.Shape[0]; r++ { + c.Logf("%v", onehot[r*npy.Shape[1]:(r+1)*npy.Shape[1]]) + } + c.Check(onehot, check.DeepEquals, []int8{ + 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, // input1 + 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, // input2 + 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, // dup/input1 + 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, // dup/input2 + }) + } + } } diff --git a/slicenumpy.go b/slicenumpy.go index 25224b2617..0fa3f2196f 100644 --- a/slicenumpy.go +++ b/slicenumpy.go @@ -63,7 +63,8 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s mergeOutput := flags.Bool("merge-output", false, "merge output into one matrix.npy and one matrix.annotations.csv") hgvsSingle := flags.Bool("single-hgvs-matrix", false, "also generate hgvs-based matrix") hgvsChunked := flags.Bool("chunked-hgvs-matrix", false, "also generate hgvs-based matrix per chromosome") - onehotChunked := flags.Bool("chunked-onehot", false, "generate one-hot tile-based matrix") + onehotSingle := flags.Bool("single-onehot", false, "generate one-hot tile-based matrix") + onehotChunked := flags.Bool("chunked-onehot", false, "generate one-hot tile-based matrix per input chunk") flags.IntVar(&cmd.threads, "threads", 16, "number of memory-hungry assembly threads") flags.StringVar(&cmd.chi2CaseControlFile, "chi2-case-control-file", "", "tsv file or directory indicating cases and controls for Χ² test (if directory, all .tsv files will be read)") flags.StringVar(&cmd.chi2CaseControlColumn, "chi2-case-control-column", "", "name of case/control column in case-control files for Χ² test (value must be 0 for control, 1 for case)") @@ -113,6 +114,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s "-merge-output=" + fmt.Sprintf("%v", *mergeOutput), "-single-hgvs-matrix=" + fmt.Sprintf("%v", *hgvsSingle), "-chunked-hgvs-matrix=" + fmt.Sprintf("%v", *hgvsChunked), + "-single-onehot=" + fmt.Sprintf("%v", *onehotSingle), "-chunked-onehot=" + fmt.Sprintf("%v", *onehotChunked), "-chi2-case-control-file=" + cmd.chi2CaseControlFile, "-chi2-case-control-column=" + cmd.chi2CaseControlColumn, @@ -314,6 +316,12 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s if *mergeOutput || *hgvsSingle { toMerge = make([][]int16, len(infiles)) } + var onehotChunks [][][]int8 + var onehotXrefs [][]onehotXref + if *onehotSingle { + onehotChunks = make([][][]int8, len(infiles)) + onehotXrefs = make([][]onehotXref, len(infiles)) + } throttleMem := throttle{Max: cmd.threads} // TODO: estimate using mem and data size throttleNumpyMem := throttle{Max: cmd.threads/2 + 1} @@ -435,7 +443,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s throttleCPU.Wait() var onehotChunk [][]int8 - var onehotXrefs []onehotXref + var onehotXref []onehotXref annotationsFilename := fmt.Sprintf("%s/matrix.%04d.annotations.csv", *outputDir, infileIdx) log.Infof("%04d: writing %s", infileIdx, annotationsFilename) @@ -465,10 +473,10 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s maxv = v } } - if *onehotChunked { + if *onehotChunked || *onehotSingle { onehot, xrefs := cmd.tv2homhet(cgs, maxv, remap, tag, tagstart) onehotChunk = append(onehotChunk, onehot...) - onehotXrefs = append(onehotXrefs, xrefs...) + onehotXref = append(onehotXref, xrefs...) } fmt.Fprintf(annow, "%d,%d,%d,=,%s,%d,,,\n", tag, outcol, rt.variant, rt.seqname, rt.pos) variants := seq[tag] @@ -570,43 +578,27 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s // transpose onehotChunk[col][row] to numpy[row*ncols+col] rows := len(cmd.cgnames) cols := len(onehotChunk) - log.Infof("%04d: preparing onehot numpy (rows=%d, cols=%d, mem=%d)", infileIdx, rows, cols, rows*cols) + log.Infof("%04d: preparing onehot numpy (rows=%d, cols=%d, mem=%d)", infileIdx, len(cmd.cgnames), len(onehotChunk), len(cmd.cgnames)*len(onehotChunk)) throttleNumpyMem.Acquire() - out := make([]int8, rows*cols) - for row := range cmd.cgnames { - out := out[row*cols:] - for colnum, values := range onehotChunk { - out[colnum] = values[row] - } - } - seq = nil - cgs = nil - debug.FreeOSMemory() - throttleNumpyMem.Release() - + out := onehotcols2int8(onehotChunk) fnm := fmt.Sprintf("%s/onehot.%04d.npy", *outputDir, infileIdx) err = writeNumpyInt8(fnm, out, rows, cols) if err != nil { return err } - fnm = fmt.Sprintf("%s/onehot-columns.%04d.npy", *outputDir, infileIdx) - xcols := len(onehotXrefs) - xdata := make([]int32, 4*xcols) - for i, xref := range onehotXrefs { - xdata[i] = int32(xref.tag) - xdata[xcols+i] = int32(xref.variant) - if xref.het { - xdata[xcols*2+i] = 1 - } - xdata[xcols*3+i] = int32(xref.pvalue * 1000000) - } - err = writeNumpyInt32(fnm, xdata, 4, xcols) + err = writeNumpyInt32(fnm, onehotXref2int32(onehotXref), 4, len(onehotXref)) if err != nil { return err } + debug.FreeOSMemory() + throttleNumpyMem.Release() } - if !*onehotChunked || *mergeOutput || *hgvsSingle { + if *onehotSingle { + onehotChunks[infileIdx] = onehotChunk + onehotXrefs[infileIdx] = onehotXref + } + if !(*onehotSingle || *onehotChunked) || *mergeOutput || *hgvsSingle { log.Infof("%04d: preparing numpy", infileIdx) throttleNumpyMem.Acquire() rows := len(cmd.cgnames) @@ -636,7 +628,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s log.Infof("%04d: matrix fragment %d rows x %d cols", infileIdx, rows, cols) toMerge[infileIdx] = out } - if !*mergeOutput && !*onehotChunked { + if !*mergeOutput && !*onehotChunked && !*onehotSingle { fnm := fmt.Sprintf("%s/matrix.%04d.npy", *outputDir, infileIdx) err = writeNumpyInt16(fnm, out, rows, cols) if err != nil { @@ -893,6 +885,26 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s } } } + if *onehotSingle { + var onehot [][]int8 + var xrefs []onehotXref + for i := range onehotChunks { + onehot = append(onehot, onehotChunks[i]...) + onehotChunks[i] = nil + xrefs = append(xrefs, onehotXrefs[i]...) + onehotXrefs[i] = nil + } + fnm := fmt.Sprintf("%s/onehot.npy", *outputDir) + err = writeNumpyInt8(fnm, onehotcols2int8(onehot), len(cmd.cgnames), len(onehot)) + if err != nil { + return 1 + } + fnm = fmt.Sprintf("%s/onehot-columns.npy", *outputDir) + err = writeNumpyInt32(fnm, onehotXref2int32(xrefs), 4, len(xrefs)) + if err != nil { + return 1 + } + } return 0 } @@ -1168,3 +1180,41 @@ func bool2int8(in []bool) []int8 { } return out } + +// convert a []onehotXref with length N to a numpy-style []int32 +// matrix with N columns, one row per field of onehotXref struct. +// +// Hom/het row contains hom=0, het=1. +// +// P-value row contains 1000000x actual p-value. +func onehotXref2int32(xrefs []onehotXref) []int32 { + xcols := len(xrefs) + xdata := make([]int32, 4*xcols) + for i, xref := range xrefs { + xdata[i] = int32(xref.tag) + xdata[xcols+i] = int32(xref.variant) + if xref.het { + xdata[xcols*2+i] = 1 + } + xdata[xcols*3+i] = int32(xref.pvalue * 1000000) + } + return xdata +} + +// transpose onehot data from in[col][row] to numpy-style +// out[row*cols+col]. +func onehotcols2int8(in [][]int8) []int8 { + if len(in) == 0 { + return nil + } + cols := len(in) + rows := len(in[0]) + out := make([]int8, rows*cols) + for row := 0; row < rows; row++ { + outrow := out[row*cols:] + for col, incol := range in { + outrow[col] = incol[row] + } + } + return out +}