+ }
+ }
+
+ if *mergeOutput || *hgvsSingle {
+ var annow *bufio.Writer
+ var annof *os.File
+ if *mergeOutput {
+ annoFilename := fmt.Sprintf("%s/matrix.annotations.csv", *outputDir)
+ annof, err = os.Create(annoFilename)
+ if err != nil {
+ return err
+ }
+ annow = bufio.NewWriterSize(annof, 1<<20)
+ }
+
+ rows := len(cmd.cgnames)
+ cols := 0
+ for _, chunk := range toMerge {
+ cols += len(chunk) / rows
+ }
+ log.Infof("merging output matrix (rows=%d, cols=%d, mem=%d) and annotations", rows, cols, rows*cols*2)
+ var out []int16
+ if *mergeOutput {
+ out = make([]int16, rows*cols)
+ }
+ hgvsCols := map[string][2][]int16{} // hgvs -> [[g0,g1,g2,...], [g0,g1,g2,...]] (slice of genomes for each phase)
+ startcol := 0
+ for outIdx, chunk := range toMerge {
+ chunkcols := len(chunk) / rows
+ if *mergeOutput {
+ 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 err
+ }
+ if *mergeOutput {
+ err = os.Remove(annotationsFilename)
+ if err != nil {
+ return err
+ }
+ }
+ for _, line := range bytes.Split(buf, []byte{'\n'}) {
+ if len(line) == 0 {
+ continue
+ }
+ fields := bytes.SplitN(line, []byte{','}, 9)
+ tag, _ := strconv.Atoi(string(fields[0]))
+ incol, _ := strconv.Atoi(string(fields[1]))
+ tileVariant, _ := strconv.Atoi(string(fields[2]))
+ hgvsID := string(fields[3])
+ seqname := string(fields[4])
+ pos, _ := strconv.Atoi(string(fields[5]))
+ refseq := fields[6]
+ if hgvsID == "" {
+ // Null entry for un-diffable
+ // tile variant
+ continue
+ }
+ if hgvsID == "=" {
+ // Null entry for ref tile
+ continue
+ }
+ if mask != nil && !mask.Check(strings.TrimPrefix(seqname, "chr"), pos, pos+len(refseq)) {
+ // The tile intersects one of
+ // the selected regions, but
+ // this particular HGVS
+ // variant does not.
+ continue
+ }
+ hgvsColPair := hgvsCols[hgvsID]
+ if hgvsColPair[0] == nil {
+ // values in new columns start
+ // out as -1 ("no data yet")
+ // or 0 ("=ref") here, may
+ // change to 1 ("hgvs variant
+ // present") below, either on
+ // this line or a future line.
+ 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)
+ return err
+ }
+ for ph := 0; ph < 2; ph++ {
+ for row := 0; row < rows; row++ {
+ v := chunk[row*chunkcols+incol*2+ph]
+ if tileVariantID(v) == rt.variant {
+ hgvsColPair[ph][row] = 0
+ } else {
+ hgvsColPair[ph][row] = -1
+ }
+ }
+ }
+ hgvsCols[hgvsID] = hgvsColPair
+ if annow != nil {
+ hgvsref := hgvs.Variant{
+ Position: pos,
+ Ref: string(refseq),
+ New: string(refseq),
+ }
+ fmt.Fprintf(annow, "%d,%d,%d,%s:g.%s,%s,%d,%s,%s,%s\n", tag, incol+startcol/2, rt.variant, seqname, hgvsref.String(), seqname, pos, refseq, refseq, fields[8])
+ }
+ }
+ if annow != nil {
+ fmt.Fprintf(annow, "%d,%d,%d,%s,%s,%d,%s,%s,%s\n", tag, incol+startcol/2, tileVariant, hgvsID, seqname, pos, refseq, fields[7], fields[8])
+ }
+ for ph := 0; ph < 2; ph++ {
+ for row := 0; row < rows; row++ {
+ v := chunk[row*chunkcols+incol*2+ph]
+ if int(v) == tileVariant {
+ hgvsColPair[ph][row] = 1
+ }
+ }
+ }
+ }
+
+ startcol += chunkcols
+ }
+ if *mergeOutput {
+ err = annow.Flush()
+ if err != nil {
+ return err
+ }
+ err = annof.Close()
+ if err != nil {
+ return err
+ }
+ err = writeNumpyInt16(fmt.Sprintf("%s/matrix.npy", *outputDir), out, rows, cols)
+ if err != nil {
+ return err
+ }
+ }
+ out = nil
+
+ if *hgvsSingle {
+ cols = len(hgvsCols) * 2
+ log.Printf("building hgvs-based matrix: %d rows x %d cols", rows, cols)
+ out = make([]int16, rows*cols)
+ hgvsIDs := make([]string, 0, cols/2)
+ for hgvsID := range hgvsCols {
+ hgvsIDs = append(hgvsIDs, hgvsID)
+ }
+ sort.Strings(hgvsIDs)
+ var hgvsLabels bytes.Buffer
+ for idx, hgvsID := range hgvsIDs {
+ fmt.Fprintf(&hgvsLabels, "%d,%s\n", idx, hgvsID)
+ for ph := 0; ph < 2; ph++ {
+ hgvscol := hgvsCols[hgvsID][ph]
+ for row, val := range hgvscol {
+ out[row*cols+idx*2+ph] = val
+ }
+ }
+ }
+ err = writeNumpyInt16(fmt.Sprintf("%s/hgvs.npy", *outputDir), out, rows, cols)
+ if err != nil {
+ return err
+ }
+
+ fnm := fmt.Sprintf("%s/hgvs.annotations.csv", *outputDir)
+ log.Printf("writing hgvs labels: %s", fnm)
+ err = ioutil.WriteFile(fnm, hgvsLabels.Bytes(), 0777)
+ if err != nil {
+ return err
+ }
+ }
+ }
+ if *onehotSingle || *onlyPCA {
+ nzCount := 0
+ for _, part := range onehotIndirect {
+ nzCount += len(part[0])
+ }
+ onehot := make([]uint32, nzCount*2) // [r,r,r,...,c,c,c,...]
+ var xrefs []onehotXref
+ chunkOffset := uint32(0)
+ outcol := 0
+ for i, part := range onehotIndirect {
+ for i := range part[1] {
+ part[1][i] += chunkOffset
+ }
+ copy(onehot[outcol:], part[0])
+ copy(onehot[outcol+nzCount:], part[1])
+ xrefs = append(xrefs, onehotXrefs[i]...)
+
+ outcol += len(part[0])
+ chunkOffset += onehotChunkSize[i]
+
+ part[0] = nil
+ part[1] = nil
+ onehotXrefs[i] = nil
+ debug.FreeOSMemory()
+ }
+ if *onehotSingle {
+ fnm := fmt.Sprintf("%s/onehot.npy", *outputDir)
+ err = writeNumpyUint32(fnm, onehot, 2, nzCount)
+ if err != nil {
+ return err
+ }
+ fnm = fmt.Sprintf("%s/onehot-columns.npy", *outputDir)
+ err = writeNumpyInt32(fnm, onehotXref2int32(xrefs), 5, len(xrefs))
+ if err != nil {
+ return err
+ }
+ fnm = fmt.Sprintf("%s/stats.json", *outputDir)
+ j, err := json.Marshal(map[string]interface{}{
+ "pvalueCallCount": cmd.pvalueCallCount,
+ })
+ if err != nil {
+ return err
+ }
+ err = os.WriteFile(fnm, j, 0777)
+ if err != nil {
+ return err
+ }
+ }
+ if *onlyPCA {
+ cols := 0
+ for _, c := range onehot[nzCount:] {
+ if int(c) >= cols {
+ cols = int(c) + 1
+ }
+ }
+ if cols == 0 {
+ return fmt.Errorf("cannot do PCA: one-hot matrix is empty")
+ }
+ log.Printf("have %d one-hot cols", cols)
+ stride := 1
+ for *maxPCATiles > 0 && cols > *maxPCATiles*2 {
+ cols = (cols + 1) / 2
+ stride = stride * 2
+ }
+ if cols%2 == 1 {
+ // we work with pairs of columns
+ cols++
+ }
+ log.Printf("creating full matrix (%d rows) and training matrix (%d rows) with %d cols, stride %d", len(cmd.cgnames), cmd.trainingSetSize, cols, stride)
+ mtxFull := mat.NewDense(len(cmd.cgnames), cols, nil)
+ mtxTrain := mat.NewDense(cmd.trainingSetSize, cols, nil)
+ for i, c := range onehot[nzCount:] {
+ if int(c/2)%stride == 0 {
+ outcol := int(c/2)/stride*2 + int(c)%2
+ mtxFull.Set(int(onehot[i]), outcol, 1)
+ if trainRow := cmd.trainingSet[int(onehot[i])]; trainRow >= 0 {
+ mtxTrain.Set(trainRow, outcol, 1)
+ }
+ }
+ }
+ log.Print("fitting")
+ transformer := nlp.NewPCA(cmd.pcaComponents)
+ transformer.Fit(mtxTrain.T())
+ log.Printf("transforming")
+ pca, err := transformer.Transform(mtxFull.T())
+ if err != nil {
+ return err
+ }
+ pca = pca.T()
+ outrows, outcols := pca.Dims()
+ log.Printf("copying result to numpy output array: %d rows, %d cols", outrows, outcols)
+ out := make([]float64, outrows*outcols)
+ for i := 0; i < outrows; i++ {
+ for j := 0; j < outcols; j++ {
+ out[i*outcols+j] = pca.At(i, j)
+ }
+ }
+ fnm := fmt.Sprintf("%s/pca.npy", *outputDir)
+ log.Printf("writing numpy: %s", fnm)
+ output, err := os.OpenFile(fnm, os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0777)
+ if err != nil {
+ return err
+ }
+ npw, err := gonpy.NewWriter(nopCloser{output})
+ if err != nil {
+ return fmt.Errorf("gonpy.NewWriter: %w", err)
+ }
+ npw.Shape = []int{outrows, outcols}
+ err = npw.WriteFloat64(out)
+ if err != nil {
+ return fmt.Errorf("WriteFloat64: %w", err)
+ }