Add numpy-common-variants sanity check.

[lightning.git] / numpycomvar.go

diff --git a/numpycomvar.go b/numpycomvar.go
new file mode 100644 (file)
index 0000000..7f99168
--- /dev/null
+++ b/numpycomvar.go
@@ -0,0 +1,108 @@
+package main
+
+import (
+       "flag"
+       "fmt"
+       "io"
+       _ "net/http/pprof"
+
+       "git.arvados.org/arvados.git/sdk/go/arvados"
+)
+
+type numpyComVar struct{}
+
+func (cmd *numpyComVar) 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)
+               }
+       }()
+       flags := flag.NewFlagSet("", flag.ContinueOnError)
+       flags.SetOutput(stderr)
+       projectUUID := flags.String("project", "", "project `UUID` for output data")
+       inputFilename := flags.String("i", "-", "numpy matrix `file`")
+       priority := flags.Int("priority", 500, "container request priority")
+       annotationsFilename := flags.String("annotations", "", "annotations tsv `file`")
+       maxResults := flags.Int("max-results", 256, "maximum number of tile variants to output")
+       minFrequency := flags.Float64("min-frequency", 0.4, "minimum allele frequency")
+       maxFrequency := flags.Float64("max-frequency", 0.6, "maximum allele frequency")
+       err = flags.Parse(args)
+       if err == flag.ErrHelp {
+               err = nil
+               return 0
+       } else if err != nil {
+               return 2
+       }
+
+       runner := arvadosContainerRunner{
+               Name:        "lightning numpy-comvar",
+               Client:      arvados.NewClientFromEnv(),
+               ProjectUUID: *projectUUID,
+               RAM:         120000000000,
+               VCPUs:       2,
+               Priority:    *priority,
+       }
+       err = runner.TranslatePaths(inputFilename, annotationsFilename)
+       if err != nil {
+               return 1
+       }
+       runner.Prog = "python3"
+       runner.Args = []string{"-c", `import sys
+import scipy
+import sys
+import csv
+
+numpyFile = sys.argv[1]
+annotationsFile = sys.argv[2]
+outputFile = sys.argv[3]
+maxResults = int(sys.argv[4])
+minFrequency = float(sys.argv[5])
+maxFrequency = float(sys.argv[6])
+
+out = open(outputFile, 'w')
+
+m = scipy.load(numpyFile)
+
+commonvariants = {}
+mincount = m.shape[0] * 2 * minFrequency
+maxcount = m.shape[0] * 2 * maxFrequency
+for tag in range(m.shape[1] // 2):
+  example = {}
+  counter = [0, 0, 0, 0, 0]
+  for genome in range(m.shape[0]):
+    for phase in range(2):
+      variant = m[genome][tag*2+phase]
+      if variant > 0 and variant < len(counter):
+        counter[variant] += 1
+        example[variant] = genome
+  for variant, count in enumerate(counter):
+    if count >= mincount and count <= maxcount:
+      commonvariants[tag,variant] = example[variant]
+      # sys.stderr.write('tag {} variant {} count {} example {} have {} commonvariants\n'.format(tag, variant, count, example[variant], len(commonvariants)))
+  if len(commonvariants) >= maxResults:
+    break
+
+found = {}
+with open(annotationsFile, newline='') as tsvfile:
+  rdr = csv.reader(tsvfile, delimiter='\t', quotechar='"')
+  for row in rdr:
+    tag = int(row[0])
+    variant = int(row[1])
+    if (tag, variant) in commonvariants:
+      found[tag, variant] = True
+      out.write(','.join(row + [str(commonvariants[tag, variant])]) + '\n')
+    elif len(found) >= len(commonvariants):
+      sys.stderr.write('done\n')
+      break
+
+out.close()
+`, *inputFilename, *annotationsFilename, "/mnt/output/commonvariants.csv", fmt.Sprintf("%d", *maxResults), fmt.Sprintf("%f", *minFrequency), fmt.Sprintf("%f", *maxFrequency)}
+       var output string
+       output, err = runner.Run()
+       if err != nil {
+               return 1
+       }
+       fmt.Fprintln(stdout, output+"/commonvariants.csv")
+       return 0
+}