19526: Output tile locations and pvalues at specified threshold.

Arvados-DCO-1.1-Signed-off-by: Tom Clegg <tom@curii.com>

diff --git a/manhattan.go b/manhattan.go
index 4ae44796be3ffba00626889e15951dc48b84dcfa..a6a58fb5cf06dbf8e7964c4366e8e68329726728 100644 (file)
--- a/manhattan.go
+++ b/manhattan.go
@@ -32,6 +32,8 @@ func (cmd *manhattanPlot) RunCommand(prog string, args []string, stdin io.Reader
        projectUUID := flags.String("project", "", "project `UUID` for output data")
        inputDirectory := flags.String("i", "-", "input `directory` (output of slice-numpy -single-onehot)")
        outputFilename := flags.String("o", "", "output `filename` (e.g., './plot.png')")
+       csvOutputFilename := flags.String("csv-output", "", "csv output `filename` (e.g., './tile-locations-pvalues.csv')")
+       csvOutputThreshold := flags.Float64("csv-output-threshold", 0, "logpvalue threshold for csv output (0 for none)")
        priority := flags.Int("priority", 500, "container request priority")
        runlocal := flags.Bool("local", false, "run on local host (default: run in an arvados container)")
        err = flags.Parse(args)
@@ -65,10 +67,13 @@ func (cmd *manhattanPlot) RunCommand(prog string, args []string, stdin io.Reader
                        return 1
                }
                *outputFilename = "/mnt/output/plot.png"
+               *csvOutputFilename = "/mnt/output/tile-locations-pvalues.csv"
        }
        args = []string{
                *inputDirectory,
                *outputFilename,
+               fmt.Sprintf("%g", *csvOutputThreshold),
+               *csvOutputFilename,
        }
        if *runlocal {
                if *outputFilename == "" {

diff --git a/manhattan.py b/manhattan.py
index 95689d972af021d3b662d4262259b2e7fd9612ce..08a4da400a544acad2e6ae73ad002b6362f87482 100644 (file)
--- a/manhattan.py
+++ b/manhattan.py
@@ -14,7 +14,11 @@ import qmplot
 (_,
  input_path,
  output_path,
+ csv_threshold_str,
+ csv_output_path,
  ) = sys.argv
+csv_threshold = float(csv_threshold_str)
+
 columns = numpy.load(os.path.join(input_path, 'onehot-columns.npy'))
 
 # pvalue maps tag# => [pvalue1, pvalue2, ...] (one het p-value and one hom p-value for each tile variant)
@@ -35,6 +39,13 @@ for dirent in os.scandir(input_path):
                 if annotation[3] == "=":
                     tilepos[int(annotation[0])] = (annotation[4], int(annotation[5]))
 
+if csv_threshold > 0 and csv_output_path != "":
+    with open(csv_output_path, 'wt') as f:
+        for tag, chrpos in sorted(tilepos.items(), key=lambda item: (item[1][0][-1] > '9', item[1][0].lstrip('chr').zfill(2), item[1][1])):
+            for p in pvalue.get(tag, []):
+                if p < pow(10, -csv_threshold):
+                    print(f'{tag},{chrpos[0]},{chrpos[1]},{p}', file=f)
+
 series = {"#CHROM": [], "POS": [], "P": []}
 for tag, chrpos in sorted(tilepos.items(), key=lambda item: (item[1][0][-1] > '9', item[1][0].lstrip('chr').zfill(2), item[1][1])):
     for p in pvalue.get(tag, []):