X-Git-Url: https://git.arvados.org/lightning.git/blobdiff_plain/b65f07682c30ec6af796d85b95b07563b65e5146..8778237e68cdbe4fd1d4b43ed5140438fc92f3a6:/slicenumpy.go

diff --git a/slicenumpy.go b/slicenumpy.go
index 8696aa53f3..0fa3f2196f 100644
--- a/slicenumpy.go
+++ b/slicenumpy.go
@@ -12,6 +12,7 @@ import (
 	"fmt"
 	"io"
 	"io/ioutil"
+	"math"
 	"net/http"
 	_ "net/http/pprof"
 	"os"
@@ -31,8 +32,14 @@ import (
 )
 
 type sliceNumpy struct {
-	filter  filter
-	threads int
+	filter                filter
+	threads               int
+	chi2CaseControlColumn string
+	chi2CaseControlFile   string
+	chi2Cases             []bool
+	chi2PValue            float64
+	minCoverage           int
+	cgnames               []string
 }
 
 func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
@@ -56,7 +63,12 @@ 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")
+	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)")
+	flags.Float64Var(&cmd.chi2PValue, "chi2-p-value", 1, "do Χ² test and omit columns with p-value above this threshold")
 	cmd.filter.Flags(flags)
 	err = flags.Parse(args)
 	if err == flag.ErrHelp {
@@ -72,6 +84,11 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 		}()
 	}
 
+	if cmd.chi2PValue != 1 && (cmd.chi2CaseControlFile == "" || cmd.chi2CaseControlColumn == "") {
+		log.Errorf("cannot use provided -chi2-p-value=%f because -chi2-case-control-file= or -chi2-case-control-column= value is empty", cmd.chi2PValue)
+		return 2
+	}
+
 	if !*runlocal {
 		runner := arvadosContainerRunner{
 			Name:        "lightning slice-numpy",
@@ -83,7 +100,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 			KeepCache:   2,
 			APIAccess:   true,
 		}
-		err = runner.TranslatePaths(inputDir, regionsFilename)
+		err = runner.TranslatePaths(inputDir, regionsFilename, &cmd.chi2CaseControlFile)
 		if err != nil {
 			return 1
 		}
@@ -97,6 +114,11 @@ 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,
+			"-chi2-p-value=" + fmt.Sprintf("%f", cmd.chi2PValue),
 		}
 		runner.Args = append(runner.Args, cmd.filter.Args()...)
 		var output string
@@ -108,7 +130,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 		return 0
 	}
 
-	infiles, err := allGobFiles(*inputDir)
+	infiles, err := allFiles(*inputDir, matchGobFile)
 	if err != nil {
 		return 1
 	}
@@ -118,7 +140,6 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 	}
 	sort.Strings(infiles)
 
-	var cgnames []string
 	var refseq map[string][]tileLibRef
 	var reftiledata = make(map[tileLibRef][]byte, 11000000)
 	in0, err := open(infiles[0])
@@ -132,6 +153,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 		return 1
 	}
 
+	cmd.cgnames = nil
 	taglen := -1
 	DecodeLibrary(in0, strings.HasSuffix(infiles[0], ".gz"), func(ent *LibraryEntry) error {
 		if len(ent.TagSet) > 0 {
@@ -144,7 +166,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 		}
 		for _, cg := range ent.CompactGenomes {
 			if matchGenome.MatchString(cg.Name) {
-				cgnames = append(cgnames, cg.Name)
+				cmd.cgnames = append(cmd.cgnames, cg.Name)
 			}
 		}
 		for _, tv := range ent.TileVariants {
@@ -166,14 +188,19 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 		err = fmt.Errorf("tagset not found")
 		return 1
 	}
-	if len(cgnames) == 0 {
+	if len(cmd.cgnames) == 0 {
 		err = fmt.Errorf("no genomes found matching regexp %q", cmd.filter.MatchGenome)
 		return 1
 	}
-	sort.Strings(cgnames)
+	sort.Strings(cmd.cgnames)
+	err = cmd.useCaseControlFiles()
+	if err != nil {
+		return 1
+	}
+	cmd.minCoverage = int(math.Ceil(cmd.filter.MinCoverage * float64(len(cmd.cgnames))))
 
 	{
-		labelsFilename := *outputDir + "/labels.csv"
+		labelsFilename := *outputDir + "/samples.csv"
 		log.Infof("writing labels to %s", labelsFilename)
 		var f *os.File
 		f, err = os.Create(labelsFilename)
@@ -181,8 +208,12 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 			return 1
 		}
 		defer f.Close()
-		for i, name := range cgnames {
-			_, err = fmt.Fprintf(f, "%d,%q\n", i, trimFilenameForLabel(name))
+		for i, name := range cmd.cgnames {
+			cc := 0
+			if cmd.chi2Cases != nil && cmd.chi2Cases[i] {
+				cc = 1
+			}
+			_, err = fmt.Fprintf(f, "%d,%q,%d\n", i, trimFilenameForLabel(name), cc)
 			if err != nil {
 				err = fmt.Errorf("write %s: %w", labelsFilename, err)
 				return 1
@@ -285,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}
@@ -294,7 +331,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 		infileIdx, infile := infileIdx, infile
 		throttleMem.Go(func() error {
 			seq := make(map[tagID][]TileVariant, 50000)
-			cgs := make(map[string]CompactGenome, len(cgnames))
+			cgs := make(map[string]CompactGenome, len(cmd.cgnames))
 			f, err := open(infile)
 			if err != nil {
 				return err
@@ -339,8 +376,8 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 			if err != nil {
 				return err
 			}
-			tagstart := cgs[cgnames[0]].StartTag
-			tagend := cgs[cgnames[0]].EndTag
+			tagstart := cgs[cmd.cgnames[0]].StartTag
+			tagend := cgs[cmd.cgnames[0]].EndTag
 
 			// TODO: filters
 
@@ -405,6 +442,9 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 			}
 			throttleCPU.Wait()
 
+			var onehotChunk [][]int8
+			var onehotXref []onehotXref
+
 			annotationsFilename := fmt.Sprintf("%s/matrix.%04d.annotations.csv", *outputDir, infileIdx)
 			log.Infof("%04d: writing %s", infileIdx, annotationsFilename)
 			annof, err := os.Create(annotationsFilename)
@@ -426,9 +466,6 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 					// mention it in annotations?)
 					continue
 				}
-				fmt.Fprintf(annow, "%d,%d,%d,=,%s,%d,,,\n", tag, outcol, rt.variant, rt.seqname, rt.pos)
-				variants := seq[tag]
-				reftilestr := strings.ToUpper(string(rt.tiledata))
 				remap := variantRemap[tag-tagstart]
 				maxv := tileVariantID(0)
 				for _, v := range remap {
@@ -436,6 +473,15 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 						maxv = v
 					}
 				}
+				if *onehotChunked || *onehotSingle {
+					onehot, xrefs := cmd.tv2homhet(cgs, maxv, remap, tag, tagstart)
+					onehotChunk = append(onehotChunk, onehot...)
+					onehotXref = append(onehotXref, xrefs...)
+				}
+				fmt.Fprintf(annow, "%d,%d,%d,=,%s,%d,,,\n", tag, outcol, rt.variant, rt.seqname, rt.pos)
+				variants := seq[tag]
+				reftilestr := strings.ToUpper(string(rt.tiledata))
+
 				done := make([]bool, maxv+1)
 				variantDiffs := make([][]hgvs.Variant, maxv+1)
 				for v, tv := range variants {
@@ -454,8 +500,10 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 						continue
 					}
 					diffs, _ := hgvs.Diff(reftilestr, strings.ToUpper(string(tv.Sequence)), 0)
+					for i := range diffs {
+						diffs[i].Position += rt.pos
+					}
 					for _, diff := range diffs {
-						diff.Position += rt.pos
 						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)
 					}
 					if *hgvsChunked {
@@ -477,12 +525,12 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 								continue
 							}
 							hgvsCol[diff] = [2][]int8{
-								make([]int8, len(cgnames)),
-								make([]int8, len(cgnames)),
+								make([]int8, len(cmd.cgnames)),
+								make([]int8, len(cmd.cgnames)),
 							}
 						}
 					}
-					for row, name := range cgnames {
+					for row, name := range cmd.cgnames {
 						variants := cgs[name].Variants[(tag-tagstart)*2:]
 						for ph := 0; ph < 2; ph++ {
 							v := variants[ph]
@@ -505,7 +553,15 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 							}
 						}
 					}
-					encodeHGVSTodo[rt.seqname] <- hgvsCol
+					for diff, colpair := range hgvsCol {
+						allele2homhet(colpair)
+						if !cmd.filterHGVScolpair(colpair) {
+							delete(hgvsCol, diff)
+						}
+					}
+					if len(hgvsCol) > 0 {
+						encodeHGVSTodo[rt.seqname] <- hgvsCol
+					}
 				}
 				outcol++
 			}
@@ -518,44 +574,69 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 				return err
 			}
 
-			log.Infof("%04d: preparing numpy", infileIdx)
-			throttleNumpyMem.Acquire()
-			rows := len(cgnames)
-			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 {
-					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[outcol] = -1
-					}
-					outcol++
+			if *onehotChunked {
+				// 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, len(cmd.cgnames), len(onehotChunk), len(cmd.cgnames)*len(onehotChunk))
+				throttleNumpyMem.Acquire()
+				out := onehotcols2int8(onehotChunk)
+				fnm := fmt.Sprintf("%s/onehot.%04d.npy", *outputDir, infileIdx)
+				err = writeNumpyInt8(fnm, out, rows, cols)
+				if err != nil {
+					return err
 				}
-			}
-			seq = nil
-			cgs = nil
-			debug.FreeOSMemory()
-			throttleNumpyMem.Release()
-
-			if *mergeOutput || *hgvsSingle {
-				log.Infof("%04d: matrix fragment %d rows x %d cols", infileIdx, rows, cols)
-				toMerge[infileIdx] = out
-			}
-			if !*mergeOutput {
-				fnm := fmt.Sprintf("%s/matrix.%04d.npy", *outputDir, infileIdx)
-				err = writeNumpyInt16(fnm, out, rows, cols)
+				fnm = fmt.Sprintf("%s/onehot-columns.%04d.npy", *outputDir, infileIdx)
+				err = writeNumpyInt32(fnm, onehotXref2int32(onehotXref), 4, len(onehotXref))
 				if err != nil {
 					return err
 				}
 				debug.FreeOSMemory()
+				throttleNumpyMem.Release()
+			}
+			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)
+				cols := 2 * outcol
+				out := make([]int16, rows*cols)
+				for row, name := range cmd.cgnames {
+					out := out[row*cols:]
+					outcol := 0
+					for col, v := range cgs[name].Variants {
+						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[outcol] = -1
+						}
+						outcol++
+					}
+				}
+				seq = nil
+				cgs = nil
+				debug.FreeOSMemory()
+				throttleNumpyMem.Release()
+				if *mergeOutput || *hgvsSingle {
+					log.Infof("%04d: matrix fragment %d rows x %d cols", infileIdx, rows, cols)
+					toMerge[infileIdx] = out
+				}
+				if !*mergeOutput && !*onehotChunked && !*onehotSingle {
+					fnm := fmt.Sprintf("%s/matrix.%04d.npy", *outputDir, infileIdx)
+					err = writeNumpyInt16(fnm, out, rows, cols)
+					if err != nil {
+						return err
+					}
+				}
 			}
+			debug.FreeOSMemory()
 			log.Infof("%s: done (%d/%d)", infile, int(atomic.AddInt64(&done, 1)), len(infiles))
 			return nil
 		})
@@ -603,13 +684,13 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 					return vi.New < vj.New
 				}
 			})
-			rows := len(cgnames)
+			rows := len(cmd.cgnames)
 			cols := len(variants) * 2
 			log.Infof("%s: building hgvs matrix (rows=%d, cols=%d, mem=%d)", seqname, rows, cols, rows*cols)
 			out := make([]int8, rows*cols)
 			for varIdx, variant := range variants {
 				hgvsCols := hgvsCols[variant]
-				for row := range cgnames {
+				for row := range cmd.cgnames {
 					for ph := 0; ph < 2; ph++ {
 						out[row*cols+varIdx+ph] = hgvsCols[ph][row]
 					}
@@ -646,7 +727,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 			annow = bufio.NewWriterSize(annof, 1<<20)
 		}
 
-		rows := len(cgnames)
+		rows := len(cmd.cgnames)
 		cols := 0
 		for _, chunk := range toMerge {
 			cols += len(chunk) / rows
@@ -715,7 +796,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
 					// change to 1 ("hgvs variant
 					// present") below, either on
 					// this line or a future line.
-					hgvsColPair = [2][]int16{make([]int16, len(cgnames)), make([]int16, len(cgnames))}
+					hgvsColPair = [2][]int16{make([]int16, len(cmd.cgnames)), make([]int16, len(cmd.cgnames))}
 					rt, ok := reftile[tagID(tag)]
 					if !ok {
 						err = fmt.Errorf("bug: seeing annotations for tag %d, but it has no reftile entry", tag)
@@ -804,9 +885,157 @@ 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
 }
 
+// Read case/control files, remove non-case/control entries from
+// cmd.cgnames, and build cmd.chi2Cases.
+func (cmd *sliceNumpy) useCaseControlFiles() error {
+	if cmd.chi2CaseControlFile == "" {
+		return nil
+	}
+	infiles, err := allFiles(cmd.chi2CaseControlFile, nil)
+	if err != nil {
+		return err
+	}
+	// index in cmd.cgnames => case(true) / control(false)
+	cc := map[int]bool{}
+	for _, infile := range infiles {
+		f, err := open(infile)
+		if err != nil {
+			return err
+		}
+		buf, err := io.ReadAll(f)
+		f.Close()
+		if err != nil {
+			return err
+		}
+		ccCol := -1
+		for _, tsv := range bytes.Split(buf, []byte{'\n'}) {
+			if len(tsv) == 0 {
+				continue
+			}
+			split := strings.Split(string(tsv), "\t")
+			if ccCol < 0 {
+				// header row
+				for col, name := range split {
+					if name == cmd.chi2CaseControlColumn {
+						ccCol = col
+						break
+					}
+				}
+				if ccCol < 0 {
+					return fmt.Errorf("%s: no column named %q in header row %q", infile, cmd.chi2CaseControlColumn, tsv)
+				}
+				continue
+			}
+			if len(split) <= ccCol {
+				continue
+			}
+			pattern := split[0]
+			found := -1
+			for i, name := range cmd.cgnames {
+				if strings.Contains(name, pattern) {
+					if found >= 0 {
+						log.Warnf("pattern %q in %s matches multiple genome IDs (%qs, %q)", pattern, infile, cmd.cgnames[found], name)
+					}
+					found = i
+				}
+			}
+			if found < 0 {
+				log.Warnf("pattern %q in %s does not match any genome IDs", pattern, infile)
+				continue
+			}
+			if split[ccCol] == "0" {
+				cc[found] = false
+			}
+			if split[ccCol] == "1" {
+				cc[found] = true
+			}
+		}
+	}
+	allnames := cmd.cgnames
+	cmd.cgnames = nil
+	cmd.chi2Cases = nil
+	ncases := 0
+	for i, name := range allnames {
+		if cc, ok := cc[i]; ok {
+			cmd.cgnames = append(cmd.cgnames, name)
+			cmd.chi2Cases = append(cmd.chi2Cases, cc)
+			if cc {
+				ncases++
+			}
+		}
+	}
+	log.Printf("%d cases, %d controls, %d neither (dropped)", ncases, len(cmd.cgnames)-ncases, len(allnames)-len(cmd.cgnames))
+	return nil
+}
+
+func (cmd *sliceNumpy) filterHGVScolpair(colpair [2][]int8) bool {
+	if cmd.chi2PValue >= 1 {
+		return true
+	}
+	col0 := make([]bool, 0, len(cmd.chi2Cases))
+	col1 := make([]bool, 0, len(cmd.chi2Cases))
+	cases := make([]bool, 0, len(cmd.chi2Cases))
+	for i, c := range cmd.chi2Cases {
+		if colpair[0][i] < 0 {
+			continue
+		}
+		col0 = append(col0, colpair[0][i] != 0)
+		col1 = append(col1, colpair[1][i] != 0)
+		cases = append(cases, c)
+	}
+	return len(cases) >= cmd.minCoverage &&
+		(pvalue(cases, col0) <= cmd.chi2PValue || pvalue(cases, col1) <= cmd.chi2PValue)
+}
+
+func writeNumpyInt32(fnm string, out []int32, 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,
+		"bytes":    rows * cols * 4,
+	}).Infof("writing numpy: %s", fnm)
+	npw.Shape = []int{rows, cols}
+	npw.WriteInt32(out)
+	err = bufw.Flush()
+	if err != nil {
+		return err
+	}
+	return output.Close()
+}
+
 func writeNumpyInt16(fnm string, out []int16, rows, cols int) error {
 	output, err := os.Create(fnm)
 	if err != nil {
@@ -822,6 +1051,7 @@ func writeNumpyInt16(fnm string, out []int16, rows, cols int) error {
 		"filename": fnm,
 		"rows":     rows,
 		"cols":     cols,
+		"bytes":    rows * cols * 2,
 	}).Infof("writing numpy: %s", fnm)
 	npw.Shape = []int{rows, cols}
 	npw.WriteInt16(out)
@@ -847,6 +1077,7 @@ func writeNumpyInt8(fnm string, out []int8, rows, cols int) error {
 		"filename": fnm,
 		"rows":     rows,
 		"cols":     cols,
+		"bytes":    rows * cols,
 	}).Infof("writing numpy: %s", fnm)
 	npw.Shape = []int{rows, cols}
 	npw.WriteInt8(out)
@@ -856,3 +1087,134 @@ func writeNumpyInt8(fnm string, out []int8, rows, cols int) error {
 	}
 	return output.Close()
 }
+
+func allele2homhet(colpair [2][]int8) {
+	a, b := colpair[0], colpair[1]
+	for i, av := range a {
+		bv := b[i]
+		if av < 0 || bv < 0 {
+			// no-call
+			a[i], b[i] = -1, -1
+		} else if av > 0 && bv > 0 {
+			// hom
+			a[i], b[i] = 1, 0
+		} else if av > 0 || bv > 0 {
+			// het
+			a[i], b[i] = 0, 1
+		} else {
+			// ref (or a different variant in same position)
+			// (this is a no-op) a[i], b[i] = 0, 0
+		}
+	}
+}
+
+type onehotXref struct {
+	tag     tagID
+	variant tileVariantID
+	het     bool
+	pvalue  float64
+}
+
+// Build onehot matrix (m[variant*2+isHet][genome] == 0 or 1) for all
+// variants of a single tile/tag#.
+//
+// Return nil if no tile variant passes Χ² filter.
+func (cmd *sliceNumpy) tv2homhet(cgs map[string]CompactGenome, maxv tileVariantID, remap []tileVariantID, tag, chunkstarttag tagID) ([][]int8, []onehotXref) {
+	if maxv < 2 {
+		// everyone has the most common variant
+		return nil, nil
+	}
+	tagoffset := tag - chunkstarttag
+	coverage := 0
+	for _, cg := range cgs {
+		if cg.Variants[tagoffset*2] > 0 && cg.Variants[tagoffset*2+1] > 0 {
+			coverage++
+		}
+	}
+	if coverage < cmd.minCoverage {
+		return nil, nil
+	}
+	obs := make([][]bool, (maxv+1)*2) // 2 slices (hom + het) for each variant#
+	for i := range obs {
+		obs[i] = make([]bool, len(cmd.cgnames))
+	}
+	for cgid, name := range cmd.cgnames {
+		cgvars := cgs[name].Variants
+		for v := tileVariantID(2); v <= maxv; v++ {
+			if remap[cgvars[tagoffset*2]] == v && remap[cgvars[tagoffset*2+1]] == v {
+				obs[v*2][cgid] = true
+			} else if remap[cgvars[tagoffset*2]] == v || remap[cgvars[tagoffset*2+1]] == v {
+				obs[v*2+1][cgid] = true
+			}
+		}
+	}
+	var onehot [][]int8
+	var xref []onehotXref
+	for homcol := 4; homcol < len(obs); homcol += 2 {
+		p := [2]float64{
+			pvalue(cmd.chi2Cases, obs[homcol]),
+			pvalue(cmd.chi2Cases, obs[homcol+1]),
+		}
+		if cmd.chi2PValue < 1 && !(p[0] < cmd.chi2PValue || p[1] < cmd.chi2PValue) {
+			continue
+		}
+		for het := 0; het < 2; het++ {
+			onehot = append(onehot, bool2int8(obs[homcol+het]))
+			xref = append(xref, onehotXref{
+				tag:     tag,
+				variant: tileVariantID(homcol / 2),
+				het:     het == 1,
+				pvalue:  p[het],
+			})
+		}
+	}
+	return onehot, xref
+}
+
+func bool2int8(in []bool) []int8 {
+	out := make([]int8, len(in))
+	for i, v := range in {
+		if v {
+			out[i] = 1
+		}
+	}
+	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
+}