X-Git-Url: https://git.arvados.org/lightning.git/blobdiff_plain/6e78b5927f608c1e090ef7f067b52287fcefebc6..d5bfb6ea362a5b906a537c315aed6c3a7ceaf0d9:/slicenumpy.go diff --git a/slicenumpy.go b/slicenumpy.go index 87cb4f6d8b..34cd777458 100644 --- a/slicenumpy.go +++ b/slicenumpy.go @@ -8,6 +8,7 @@ import ( "bufio" "bytes" "encoding/gob" + "encoding/json" "errors" "flag" "fmt" @@ -28,40 +29,56 @@ import ( "git.arvados.org/arvados.git/sdk/go/arvados" "github.com/arvados/lightning/hgvs" + "github.com/james-bowman/nlp" "github.com/kshedden/gonpy" "github.com/sirupsen/logrus" log "github.com/sirupsen/logrus" "golang.org/x/crypto/blake2b" + "gonum.org/v1/gonum/mat" ) const annotationMaxTileSpan = 100 type sliceNumpy struct { - filter filter - threads int - chi2CaseControlColumn string - chi2CaseControlFile string - chi2Cases []bool - chi2PValue float64 - minCoverage int - cgnames []string - includeVariant1 bool - debugTag tagID + filter filter + threads int + chi2Cases []bool + chi2PValue float64 + pvalueMinFrequency float64 + maxFrequency float64 + pcaComponents int + minCoverage int + minCoverageAll bool + includeVariant1 bool + debugTag tagID + + cgnames []string + samples []sampleInfo + trainingSet []int // samples index => training set index, or -1 if not in training set + trainingSetSize int + pvalue func(onehot []bool) float64 + pvalueCallCount int64 } func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int { - var err error - defer func() { - if err != nil { - fmt.Fprintf(stderr, "%s\n", err) - } - }() + err := cmd.run(prog, args, stdin, stdout, stderr) + if err != nil { + fmt.Fprintf(stderr, "%s\n", err) + return 1 + } + return 0 +} + +func (cmd *sliceNumpy) run(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) error { flags := flag.NewFlagSet("", flag.ContinueOnError) flags.SetOutput(stderr) pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`") runlocal := flags.Bool("local", false, "run on local host (default: run in an arvados container)") + arvadosRAM := flags.Int("arvados-ram", 750000000000, "amount of memory to request for arvados container (`bytes`)") + arvadosVCPUs := flags.Int("arvados-vcpus", 96, "number of VCPUs to request for arvados container") projectUUID := flags.String("project", "", "project `UUID` for output data") priority := flags.Int("priority", 500, "container request priority") + preemptible := flags.Bool("preemptible", true, "request preemptible instance") inputDir := flags.String("input-dir", "./in", "input `directory`") outputDir := flags.String("output-dir", "./out", "output `directory`") ref := flags.String("ref", "", "reference name (if blank, choose last one that appears in input)") @@ -72,20 +89,26 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s 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") - onlyPCA := flags.Bool("pca", false, "generate pca matrix") + samplesFilename := flags.String("samples", "", "`samples.csv` file with training/validation and case/control groups (see 'lightning choose-samples')") + caseControlOnly := flags.Bool("case-control-only", false, "drop samples that are not in case/control groups") + onlyPCA := flags.Bool("pca", false, "run principal component analysis, write components to pca.npy and samples.csv") + flags.IntVar(&cmd.pcaComponents, "pca-components", 4, "number of PCA components to compute / use in logistic regression") + maxPCATiles := flags.Int("max-pca-tiles", 0, "maximum tiles to use as PCA input (filter, then drop every 2nd colum pair until below max)") debugTag := flags.Int("debug-tag", -1, "log debugging details about specified tag") - 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") + flags.BoolVar(&cmd.minCoverageAll, "min-coverage-all", false, "apply -min-coverage filter based on all samples, not just training set") + flags.IntVar(&cmd.threads, "threads", 16, "number of memory-hungry assembly threads, and number of VCPUs to request for arvados container") + flags.Float64Var(&cmd.chi2PValue, "chi2-p-value", 1, "do Χ² test (or logistic regression if -samples file has PCA components) and omit columns with p-value above this threshold") + flags.Float64Var(&cmd.pvalueMinFrequency, "pvalue-min-frequency", 0.01, "skip p-value calculation on tile variants below this frequency in the training set") + flags.Float64Var(&cmd.maxFrequency, "max-frequency", 1, "do not output variants above this frequency in the training set") flags.BoolVar(&cmd.includeVariant1, "include-variant-1", false, "include most common variant when building one-hot matrix") cmd.filter.Flags(flags) - err = flags.Parse(args) + err := flags.Parse(args) if err == flag.ErrHelp { - err = nil - return 0 + return nil } else if err != nil { - return 2 + return err + } else if flags.NArg() > 0 { + return fmt.Errorf("errant command line arguments after parsed flags: %v", flags.Args()) } if *pprof != "" { @@ -94,9 +117,8 @@ 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 cmd.chi2PValue != 1 && *samplesFilename == "" { + return fmt.Errorf("cannot use provided -chi2-p-value=%f because -samples= value is empty", cmd.chi2PValue) } cmd.debugTag = tagID(*debugTag) @@ -106,15 +128,16 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s Name: "lightning slice-numpy", Client: arvados.NewClientFromEnv(), ProjectUUID: *projectUUID, - RAM: 750000000000, - VCPUs: 96, + RAM: int64(*arvadosRAM), + VCPUs: *arvadosVCPUs, Priority: *priority, KeepCache: 2, APIAccess: true, + Preemptible: *preemptible, } - err = runner.TranslatePaths(inputDir, regionsFilename, &cmd.chi2CaseControlFile) + err = runner.TranslatePaths(inputDir, regionsFilename, samplesFilename) if err != nil { - return 1 + return err } runner.Args = []string{"slice-numpy", "-local=true", "-pprof=:6060", @@ -128,10 +151,15 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s "-chunked-hgvs-matrix=" + fmt.Sprintf("%v", *hgvsChunked), "-single-onehot=" + fmt.Sprintf("%v", *onehotSingle), "-chunked-onehot=" + fmt.Sprintf("%v", *onehotChunked), + "-samples=" + *samplesFilename, + "-case-control-only=" + fmt.Sprintf("%v", *caseControlOnly), + "-min-coverage-all=" + fmt.Sprintf("%v", cmd.minCoverageAll), "-pca=" + fmt.Sprintf("%v", *onlyPCA), - "-chi2-case-control-file=" + cmd.chi2CaseControlFile, - "-chi2-case-control-column=" + cmd.chi2CaseControlColumn, + "-pca-components=" + fmt.Sprintf("%d", cmd.pcaComponents), + "-max-pca-tiles=" + fmt.Sprintf("%d", *maxPCATiles), "-chi2-p-value=" + fmt.Sprintf("%f", cmd.chi2PValue), + "-pvalue-min-frequency=" + fmt.Sprintf("%f", cmd.pvalueMinFrequency), + "-max-frequency=" + fmt.Sprintf("%f", cmd.maxFrequency), "-include-variant-1=" + fmt.Sprintf("%v", cmd.includeVariant1), "-debug-tag=" + fmt.Sprintf("%d", cmd.debugTag), } @@ -139,19 +167,19 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s var output string output, err = runner.Run() if err != nil { - return 1 + return err } fmt.Fprintln(stdout, output) - return 0 + return nil } infiles, err := allFiles(*inputDir, matchGobFile) if err != nil { - return 1 + return err } if len(infiles) == 0 { err = fmt.Errorf("no input files found in %s", *inputDir) - return 1 + return err } sort.Strings(infiles) @@ -159,13 +187,22 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s var reftiledata = make(map[tileLibRef][]byte, 11000000) in0, err := open(infiles[0]) if err != nil { - return 1 + return err } matchGenome, err := regexp.Compile(cmd.filter.MatchGenome) if err != nil { err = fmt.Errorf("-match-genome: invalid regexp: %q", cmd.filter.MatchGenome) - return 1 + return err + } + + if *samplesFilename != "" { + cmd.samples, err = loadSampleInfo(*samplesFilename) + if err != nil { + return err + } + } else if *caseControlOnly { + return fmt.Errorf("-case-control-only does not make sense without -samples") } cmd.cgnames = nil @@ -192,75 +229,114 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s return nil }) if err != nil { - return 1 + return err } in0.Close() if refseq == nil { err = fmt.Errorf("%s: reference sequence not found", infiles[0]) - return 1 + return err } if len(tagset) == 0 { err = fmt.Errorf("tagset not found") - return 1 + return err } taglib := &tagLibrary{} err = taglib.setTags(tagset) if err != nil { - return 1 + return err } taglen := taglib.TagLen() - - if len(cmd.cgnames) == 0 { - err = fmt.Errorf("no genomes found matching regexp %q", cmd.filter.MatchGenome) - return 1 - } sort.Strings(cmd.cgnames) - err = cmd.useCaseControlFiles() - if err != nil { - return 1 - } + if len(cmd.cgnames) == 0 { - err = fmt.Errorf("fatal: 0 cases, 0 controls, nothing to do") - return 1 - } - if cmd.filter.MinCoverage == 1 { - // In the generic formula below, floating point - // arithmetic can effectively push the coverage - // threshold above 1.0, which is impossible/useless. - // 1.0 needs to mean exactly 100% coverage. - cmd.minCoverage = len(cmd.cgnames) - } else { - cmd.minCoverage = int(math.Ceil(cmd.filter.MinCoverage * float64(len(cmd.cgnames)))) + return fmt.Errorf("fatal: 0 matching samples in library, nothing to do") } - - { - labelsFilename := *outputDir + "/samples.csv" - log.Infof("writing labels to %s", labelsFilename) - var f *os.File - f, err = os.Create(labelsFilename) - if err != nil { - return 1 + cmd.trainingSet = make([]int, len(cmd.cgnames)) + if *samplesFilename == "" { + cmd.trainingSetSize = len(cmd.cgnames) + for i, name := range cmd.cgnames { + cmd.samples = append(cmd.samples, sampleInfo{ + id: trimFilenameForLabel(name), + isTraining: true, + }) + cmd.trainingSet[i] = i } - defer f.Close() + } else if len(cmd.cgnames) != len(cmd.samples) { + return fmt.Errorf("mismatched sample list: %d samples in library, %d in %s", len(cmd.cgnames), len(cmd.samples), *samplesFilename) + } else { for i, name := range cmd.cgnames { - cc := 0 - if cmd.chi2Cases != nil && cmd.chi2Cases[i] { - cc = 1 + if s := trimFilenameForLabel(name); s != cmd.samples[i].id { + return fmt.Errorf("mismatched sample list: sample %d is %q in library, %q in %s", i, s, cmd.samples[i].id, *samplesFilename) } - _, 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 + } + if *caseControlOnly { + for i := 0; i < len(cmd.samples); i++ { + if !cmd.samples[i].isTraining && !cmd.samples[i].isValidation { + if i+1 < len(cmd.samples) { + copy(cmd.samples[i:], cmd.samples[i+1:]) + copy(cmd.cgnames[i:], cmd.cgnames[i+1:]) + } + cmd.samples = cmd.samples[:len(cmd.samples)-1] + cmd.cgnames = cmd.cgnames[:len(cmd.cgnames)-1] + i-- + } } } - err = f.Close() + cmd.chi2Cases = nil + cmd.trainingSetSize = 0 + for i := range cmd.cgnames { + if cmd.samples[i].isTraining { + cmd.trainingSet[i] = cmd.trainingSetSize + cmd.trainingSetSize++ + cmd.chi2Cases = append(cmd.chi2Cases, cmd.samples[i].isCase) + } else { + cmd.trainingSet[i] = -1 + } + } + if cmd.pvalue == nil { + cmd.pvalue = func(onehot []bool) float64 { + return pvalue(onehot, cmd.chi2Cases) + } + } + } + + if cmd.minCoverageAll { + cmd.minCoverage = len(cmd.cgnames) + } else { + cmd.minCoverage = cmd.trainingSetSize + } + if cmd.filter.MinCoverage < 1 { + cmd.minCoverage = int(math.Ceil(cmd.filter.MinCoverage * float64(cmd.minCoverage))) + } + + if len(cmd.samples[0].pcaComponents) > 0 { + cmd.pvalue = glmPvalueFunc(cmd.samples, cmd.pcaComponents) + // Unfortunately, statsmodel/glm lib logs stuff to + // os.Stdout when it panics on an unsolvable + // problem. We recover() from the panic in glm.go, but + // we also need to commandeer os.Stdout to avoid + // producing large quantities of logs. + stdoutWas := os.Stdout + defer func() { os.Stdout = stdoutWas }() + os.Stdout, err = os.Open(os.DevNull) if err != nil { - err = fmt.Errorf("close %s: %w", labelsFilename, err) - return 1 + return err } } + // cgnamemap[name]==true for samples that we are including in + // output + cgnamemap := map[string]bool{} + for _, name := range cmd.cgnames { + cgnamemap[name] = true + } + + err = writeSampleInfo(cmd.samples, *outputDir) + if err != nil { + return err + } + log.Info("indexing reference tiles") type reftileinfo struct { variant tileVariantID @@ -282,10 +358,10 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s tiledata := reftiledata[libref] if len(tiledata) == 0 { err = fmt.Errorf("missing tiledata for tag %d variant %d in %s in ref", libref.Tag, libref.Variant, seqname) - return 1 + return err } foundthistag := false - taglib.FindAll(tiledata[:len(tiledata)-1], func(tagid tagID, offset, _ int) { + taglib.FindAll(bufio.NewReader(bytes.NewReader(tiledata[:len(tiledata)-1])), nil, func(tagid tagID, offset, _ int) { if !foundthistag && tagid == libref.Tag { foundthistag = true return @@ -328,7 +404,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s log.Printf("loading regions from %s", *regionsFilename) mask, err = makeMask(*regionsFilename, *expandRegions) if err != nil { - return 1 + return err } log.Printf("before applying mask, len(reftile) == %d", len(reftile)) log.Printf("deleting reftile entries for regions outside %d intervals", mask.Len()) @@ -349,7 +425,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s var f *os.File f, err = os.Create(*outputDir + "/tmp." + seqname + ".gob") if err != nil { - return 1 + return err } defer os.Remove(f.Name()) bufw := bufio.NewWriterSize(f, 1<<24) @@ -442,7 +518,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s if cmd.filter.MaxTag >= 0 && cg.StartTag > tagID(cmd.filter.MaxTag) { return errSkip } - if !matchGenome.MatchString(cg.Name) { + if !cgnamemap[cg.Name] { continue } // pad to full slice size @@ -458,7 +534,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s if err == errSkip { return nil } else if err != nil { - return err + return fmt.Errorf("%04d: DecodeLibrary(%s): %w", infileIdx, infile, err) } tagstart := cgs[cmd.cgnames[0]].StartTag tagend := cgs[cmd.cgnames[0]].EndTag @@ -480,7 +556,11 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s count[blake2b.Sum256(rt.tiledata)] = 0 } - for cgname, cg := range cgs { + for cgidx, cgname := range cmd.cgnames { + if !cmd.minCoverageAll && !cmd.samples[cgidx].isTraining { + continue + } + cg := cgs[cgname] idx := int(tag-tagstart) * 2 for allele := 0; allele < 2; allele++ { v := cg.Variants[idx+allele] @@ -598,6 +678,12 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s break } remap := variantRemap[tag-tagstart] + if remap == nil { + // was not assigned above, + // because minCoverage + outcol++ + continue + } maxv := tileVariantID(0) for _, v := range remap { if maxv < v { @@ -786,7 +872,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s if cmd.filter.MaxTag >= 0 && tag > tagID(cmd.filter.MaxTag) { break } - if rt := reftile[tag]; rt == nil || rt.excluded { + if rt := reftile[tag]; mask != nil && (rt == nil || rt.excluded) { continue } if v == 0 { @@ -824,7 +910,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s }) } if err = throttleMem.Wait(); err != nil { - return 1 + return err } if *hgvsChunked { @@ -834,14 +920,14 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s } err = encodeHGVS.Wait() if err != nil { - return 1 + return err } for seqname := range refseq { log.Infof("%s: reading hgvsCols from temp file", seqname) f := tmpHGVSCols[seqname] _, err = f.Seek(0, io.SeekStart) if err != nil { - return 1 + return err } var hgvsCols hgvsColSet dec := gob.NewDecoder(bufio.NewReaderSize(f, 1<<24)) @@ -849,7 +935,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s err = dec.Decode(&hgvsCols) } if err != io.EOF { - return 1 + return err } log.Infof("%s: sorting %d hgvs variants", seqname, len(hgvsCols)) variants := make([]hgvs.Variant, 0, len(hgvsCols)) @@ -880,7 +966,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s } err = writeNumpyInt8(fmt.Sprintf("%s/hgvs.%s.npy", *outputDir, seqname), out, rows, cols) if err != nil { - return 1 + return err } out = nil @@ -892,7 +978,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s } err = ioutil.WriteFile(fnm, hgvsLabels.Bytes(), 0666) if err != nil { - return 1 + return err } } } @@ -904,7 +990,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s annoFilename := fmt.Sprintf("%s/matrix.annotations.csv", *outputDir) annof, err = os.Create(annoFilename) if err != nil { - return 1 + return err } annow = bufio.NewWriterSize(annof, 1<<20) } @@ -934,12 +1020,12 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s log.Infof("reading %s", annotationsFilename) buf, err := os.ReadFile(annotationsFilename) if err != nil { - return 1 + return err } if *mergeOutput { err = os.Remove(annotationsFilename) if err != nil { - return 1 + return err } } for _, line := range bytes.Split(buf, []byte{'\n'}) { @@ -982,7 +1068,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s rt, ok := reftile[tagID(tag)] if !ok { err = fmt.Errorf("bug: seeing annotations for tag %d, but it has no reftile entry", tag) - return 1 + return err } for ph := 0; ph < 2; ph++ { for row := 0; row < rows; row++ { @@ -1022,15 +1108,15 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s if *mergeOutput { err = annow.Flush() if err != nil { - return 1 + return err } err = annof.Close() if err != nil { - return 1 + return err } err = writeNumpyInt16(fmt.Sprintf("%s/matrix.npy", *outputDir), out, rows, cols) if err != nil { - return 1 + return err } } out = nil @@ -1056,18 +1142,18 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s } err = writeNumpyInt16(fmt.Sprintf("%s/hgvs.npy", *outputDir), out, rows, cols) if err != nil { - return 1 + 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 1 + return err } } } - if *onehotSingle { + if *onehotSingle || *onlyPCA { nzCount := 0 for _, part := range onehotIndirect { nzCount += len(part[0]) @@ -1092,19 +1178,113 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s onehotXrefs[i] = nil debug.FreeOSMemory() } - fnm := fmt.Sprintf("%s/onehot.npy", *outputDir) - err = writeNumpyUint32(fnm, onehot, 2, nzCount) - if err != nil { - return 1 - } - fnm = fmt.Sprintf("%s/onehot-columns.npy", *outputDir) - err = writeNumpyInt32(fnm, onehotXref2int32(xrefs), 5, len(xrefs)) - if err != nil { - return 1 + 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 { + 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[c] = 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" @@ -1112,105 +1292,137 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s var f *os.File f, err = os.Create(tagoffsetFilename) if err != nil { - return 1 + return err } defer f.Close() for idx, offset := range chunkStartTag { _, err = fmt.Fprintf(f, "%q,%d\n", fmt.Sprintf("matrix.%04d.npy", idx), offset) if err != nil { err = fmt.Errorf("write %s: %w", tagoffsetFilename, err) - return 1 + return err } } err = f.Close() if err != nil { err = fmt.Errorf("close %s: %w", tagoffsetFilename, err) - return 1 + return err } } - return 0 + + return nil } -// 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 +type sampleInfo struct { + id string + isCase bool + isControl bool + isTraining bool + isValidation bool + pcaComponents []float64 +} + +// Read samples.csv file with case/control and training/validation +// flags. +func loadSampleInfo(samplesFilename string) ([]sampleInfo, error) { + var si []sampleInfo + f, err := open(samplesFilename) + if err != nil { + return nil, err } - infiles, err := allFiles(cmd.chi2CaseControlFile, nil) + buf, err := io.ReadAll(f) + f.Close() if err != nil { - return err + return nil, 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 + lineNum := 0 + for _, csv := range bytes.Split(buf, []byte{'\n'}) { + lineNum++ + if len(csv) == 0 { + continue } - buf, err := io.ReadAll(f) - f.Close() - if err != nil { - return err + split := strings.Split(string(csv), ",") + if len(split) < 4 { + return nil, fmt.Errorf("%d fields < 4 in %s line %d: %q", len(split), samplesFilename, lineNum, csv) } - 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 + if split[0] == "Index" && split[1] == "SampleID" && split[2] == "CaseControl" && split[3] == "TrainingValidation" { + continue + } + idx, err := strconv.Atoi(split[0]) + if err != nil { + if lineNum == 1 { + return nil, fmt.Errorf("header does not look right: %q", csv) } - 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 + return nil, fmt.Errorf("%s line %d: index: %s", samplesFilename, lineNum, err) + } + if idx != len(si) { + return nil, fmt.Errorf("%s line %d: index %d out of order", samplesFilename, lineNum, idx) + } + var pcaComponents []float64 + if len(split) > 4 { + for _, s := range split[4:] { + f, err := strconv.ParseFloat(s, 64) + if err != nil { + return nil, fmt.Errorf("%s line %d: cannot parse float %q: %s", samplesFilename, lineNum, s, err) } - } - 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 + pcaComponents = append(pcaComponents, f) } } + si = append(si, sampleInfo{ + id: split[1], + isCase: split[2] == "1", + isControl: split[2] == "0", + isTraining: split[3] == "1", + isValidation: split[3] == "0" && len(split[2]) > 0, // fix errant 0s in input + pcaComponents: pcaComponents, + }) } - 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++ - } + return si, nil +} + +func writeSampleInfo(samples []sampleInfo, outputDir string) error { + fnm := outputDir + "/samples.csv" + log.Infof("writing sample metadata to %s", fnm) + f, err := os.Create(fnm) + if err != nil { + return err + } + defer f.Close() + pcaLabels := "" + if len(samples) > 0 { + for i := range samples[0].pcaComponents { + pcaLabels += fmt.Sprintf(",PCA%d", i) + } + } + _, err = fmt.Fprintf(f, "Index,SampleID,CaseControl,TrainingValidation%s\n", pcaLabels) + if err != nil { + return err + } + for i, si := range samples { + var cc, tv string + if si.isCase { + cc = "1" + } else if si.isControl { + cc = "0" + } + if si.isTraining { + tv = "1" + } else if si.isValidation { + tv = "0" + } + var pcavals string + for _, pcaval := range si.pcaComponents { + pcavals += fmt.Sprintf(",%f", pcaval) + } + _, err = fmt.Fprintf(f, "%d,%s,%s,%s%s\n", i, si.id, cc, tv, pcavals) + if err != nil { + return fmt.Errorf("write %s: %w", fnm, err) } } - log.Printf("%d cases, %d controls, %d neither (dropped)", ncases, len(cmd.cgnames)-ncases, len(allnames)-len(cmd.cgnames)) + err = f.Close() + if err != nil { + return fmt.Errorf("close %s: %w", fnm, err) + } + log.Print("done") return nil } @@ -1362,6 +1574,7 @@ type onehotXref struct { variant tileVariantID hom bool pvalue float64 + maf float64 } const onehotXrefSize = unsafe.Sizeof(onehotXref{}) @@ -1390,7 +1603,11 @@ func (cmd *sliceNumpy) tv2homhet(cgs map[string]CompactGenome, maxv tileVariantI } tagoffset := tag - chunkstarttag coverage := 0 - for _, cg := range cgs { + for cgidx, cgname := range cmd.cgnames { + if !cmd.minCoverageAll && !cmd.samples[cgidx].isTraining { + continue + } + cg := cgs[cgname] alleles := 0 for _, v := range cg.Variants[tagoffset*2 : tagoffset*2+2] { if v > 0 && int(v) < len(seq[tag]) && len(seq[tag][v].Sequence) > 0 { @@ -1404,23 +1621,36 @@ func (cmd *sliceNumpy) tv2homhet(cgs map[string]CompactGenome, maxv tileVariantI if coverage < cmd.minCoverage { return nil, nil } + // "observed" array for p-value calculation (training set + // only) obs := make([][]bool, (maxv+1)*2) // 2 slices (hom + het) for each variant# + // one-hot output (all samples) + outcols := make([][]int8, (maxv+1)*2) for i := range obs { - obs[i] = make([]bool, len(cmd.cgnames)) + obs[i] = make([]bool, cmd.trainingSetSize) + outcols[i] = make([]int8, len(cmd.cgnames)) } for cgid, name := range cmd.cgnames { + tsid := cmd.trainingSet[cgid] cgvars := cgs[name].Variants[tagoffset*2:] tv0, tv1 := remap[cgvars[0]], remap[cgvars[1]] for v := tileVariantID(1); v <= maxv; v++ { if tv0 == v && tv1 == v { - obs[v*2][cgid] = true + if tsid >= 0 { + obs[v*2][tsid] = true + } + outcols[v*2][cgid] = 1 } else if tv0 == v || tv1 == v { - obs[v*2+1][cgid] = true + if tsid >= 0 { + obs[v*2+1][tsid] = true + } + outcols[v*2+1][cgid] = 1 } } } var onehot [][]int8 var xref []onehotXref + var maf float64 for col := 2; col < len(obs); col++ { // col 0,1 correspond to tile variant 0, i.e., // no-call; col 2,3 correspond to the most common @@ -1428,29 +1658,53 @@ func (cmd *sliceNumpy) tv2homhet(cgs map[string]CompactGenome, maxv tileVariantI if col < 4 && !cmd.includeVariant1 { continue } - p := pvalue(obs[col], cmd.chi2Cases) + if col&1 == 0 { + maf = homhet2maf(obs[col : col+2]) + if maf < cmd.pvalueMinFrequency { + // Skip both columns (hom and het) if + // allele frequency is below threshold + col++ + continue + } + if maf > cmd.maxFrequency { + // Skip both columns if allele + // frequency is above threshold + col++ + continue + } + } + atomic.AddInt64(&cmd.pvalueCallCount, 1) + p := cmd.pvalue(obs[col]) if cmd.chi2PValue < 1 && !(p < cmd.chi2PValue) { continue } - onehot = append(onehot, bool2int8(obs[col])) + onehot = append(onehot, outcols[col]) xref = append(xref, onehotXref{ tag: tag, variant: tileVariantID(col >> 1), hom: col&1 == 0, pvalue: p, + maf: maf, }) } return onehot, xref } -func bool2int8(in []bool) []int8 { - out := make([]int8, len(in)) - for i, v := range in { - if v { - out[i] = 1 +func homhet2maf(onehot [][]bool) float64 { + if len(onehot[0]) == 0 { + return 0 + } + n := 0 + for i := range onehot[0] { + if onehot[0][i] { + // hom + n += 2 + } else if onehot[1][i] { + // het + n += 1 } } - return out + return float64(n) / float64(len(onehot[0])*2) } // convert a []onehotXref with length N to a numpy-style []int32 @@ -1461,7 +1715,7 @@ func bool2int8(in []bool) []int8 { // P-value row contains 1000000x actual p-value. func onehotXref2int32(xrefs []onehotXref) []int32 { xcols := len(xrefs) - xdata := make([]int32, 5*xcols) + xdata := make([]int32, 6*xcols) for i, xref := range xrefs { xdata[i] = int32(xref.tag) xdata[xcols+i] = int32(xref.variant) @@ -1470,6 +1724,7 @@ func onehotXref2int32(xrefs []onehotXref) []int32 { } xdata[xcols*3+i] = int32(xref.pvalue * 1000000) xdata[xcols*4+i] = int32(-math.Log10(xref.pvalue) * 1000000) + xdata[xcols*5+i] = int32(xref.maf * 1000000) } return xdata }