numpycomvar.go

   1 package main
   2
   3 import (
   4         "flag"
   5         "fmt"
   6         "io"
   7         _ "net/http/pprof"
   8
   9         "git.arvados.org/arvados.git/sdk/go/arvados"
  10 )
  11
  12 type numpyComVar struct{}
  13
  14 func (cmd *numpyComVar) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
  15         var err error
  16         defer func() {
  17                 if err != nil {
  18                         fmt.Fprintf(stderr, "%s\n", err)
  19                 }
  20         }()
  21         flags := flag.NewFlagSet("", flag.ContinueOnError)
  22         flags.SetOutput(stderr)
  23         projectUUID := flags.String("project", "", "project `UUID` for output data")
  24         inputFilename := flags.String("i", "-", "numpy matrix `file`")
  25         priority := flags.Int("priority", 500, "container request priority")
  26         annotationsFilename := flags.String("annotations", "", "annotations tsv `file`")
  27         maxResults := flags.Int("max-results", 256, "maximum number of tile variants to output")
  28         minFrequency := flags.Float64("min-frequency", 0.4, "minimum allele frequency")
  29         maxFrequency := flags.Float64("max-frequency", 0.6, "maximum allele frequency")
  30         err = flags.Parse(args)
  31         if err == flag.ErrHelp {
  32                 err = nil
  33                 return 0
  34         } else if err != nil {
  35                 return 2
  36         }
  37
  38         runner := arvadosContainerRunner{
  39                 Name:        "lightning numpy-comvar",
  40                 Client:      arvados.NewClientFromEnv(),
  41                 ProjectUUID: *projectUUID,
  42                 RAM:         120000000000,
  43                 VCPUs:       2,
  44                 Priority:    *priority,
  45         }
  46         err = runner.TranslatePaths(inputFilename, annotationsFilename)
  47         if err != nil {
  48                 return 1
  49         }
  50         runner.Prog = "python3"
  51         runner.Args = []string{"-c", `import sys
  52 import scipy
  53 import sys
  54 import csv
  55
  56 numpyFile = sys.argv[1]
  57 annotationsFile = sys.argv[2]
  58 outputFile = sys.argv[3]
  59 maxResults = int(sys.argv[4])
  60 minFrequency = float(sys.argv[5])
  61 maxFrequency = float(sys.argv[6])
  62
  63 out = open(outputFile, 'w')
  64
  65 m = scipy.load(numpyFile)
  66
  67 commonvariants = {}
  68 mincount = m.shape[0] * 2 * minFrequency
  69 maxcount = m.shape[0] * 2 * maxFrequency
  70 for tag in range(m.shape[1] // 2):
  71   example = {}
  72   counter = [0, 0, 0, 0, 0]
  73   for genome in range(m.shape[0]):
  74     for phase in range(2):
  75       variant = m[genome][tag*2+phase]
  76       if variant > 0 and variant < len(counter):
  77         counter[variant] += 1
  78         example[variant] = genome
  79   for variant, count in enumerate(counter):
  80     if count >= mincount and count <= maxcount:
  81       commonvariants[tag,variant] = example[variant]
  82       # sys.stderr.write('tag {} variant {} count {} example {} have {} commonvariants\n'.format(tag, variant, count, example[variant], len(commonvariants)))
  83   if len(commonvariants) >= maxResults:
  84     break
  85
  86 found = {}
  87 with open(annotationsFile, newline='') as tsvfile:
  88   rdr = csv.reader(tsvfile, delimiter='\t', quotechar='"')
  89   for row in rdr:
  90     tag = int(row[0])
  91     variant = int(row[1])
  92     if (tag, variant) in commonvariants:
  93       found[tag, variant] = True
  94       out.write(','.join(row + [str(commonvariants[tag, variant])]) + '\n')
  95     elif len(found) >= len(commonvariants):
  96       sys.stderr.write('done\n')
  97       break
  98
  99 out.close()
 100 `, *inputFilename, *annotationsFilename, "/mnt/output/commonvariants.csv", fmt.Sprintf("%d", *maxResults), fmt.Sprintf("%f", *minFrequency), fmt.Sprintf("%f", *maxFrequency)}
 101         var output string
 102         output, err = runner.Run()
 103         if err != nil {
 104                 return 1
 105         }
 106         fmt.Fprintln(stdout, output+"/commonvariants.csv")
 107         return 0
 108 }