+ 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)
+ }
+ err = output.Close()
+ if err != nil {
+ return err
+ }
+ log.Print("done")
+
+ log.Print("copying pca components to sampleInfo")
+ for i := range cmd.samples {
+ cmd.samples[i].pcaComponents = make([]float64, outcols)
+ for c := 0; c < outcols; c++ {
+ cmd.samples[i].pcaComponents[i] = pca.At(i, c)
+ }
+ }
+ log.Print("done")
+
+ err = writeSampleInfo(cmd.samples, *outputDir)
+ if err != nil {
+ return err
+ }
+ }
+ }
+ if !*mergeOutput && !*onehotChunked && !*onehotSingle && !*onlyPCA {
+ tagoffsetFilename := *outputDir + "/chunk-tag-offset.csv"
+ log.Infof("writing tag offsets to %s", tagoffsetFilename)
+ var f *os.File
+ f, err = os.Create(tagoffsetFilename)