Read pheno JSON file and write out pheno file for gemma

3 files changed, 42 insertions, 28 deletions

diff --git a/gn/runner/gemma.scm b/gn/runner/gemma.scm
index 131284a..4c1c374 100644
--- a/gn/runner/gemma.scm
+++ b/gn/runner/gemma.scm
@@ -10,16 +10,37 @@
   #:use-module (rnrs base)
 
   #:export (
+            write-pheno-file
             run-gemma
             ))
 
+(define (write-pheno-file fn traits)
+  (define bxd-inds (geno-inds-bxd "BXD.json"))
+  (assert (= 235 (length bxd-inds)))
+  (display bxd-inds)
+  (call-with-output-file fn
+    (lambda (port)
+      (for-each
+       (lambda (ind)
+         (begin
+           (let* [(value (assoc-ref traits ind))
+                  (outvalue (if value
+                                value
+                                "NA"))]
+             (format #t "~s ~s" ind outvalue)
+             (newline)
+             (display outvalue port)
+             (newline port))))
+       bxd-inds)
+      (close port)
+      ))
+)
+
 (define (run-gemma population data-id probeset-id probesetfreeze-id name trait-name traits)
   "Run gemma-wrapper to compute GRM and GWA. On failure the run will stop(!)"
-  (define bxd-inds (geno-inds-bxd "BXD.json"))
   (define ids (string-append (int-to-string data-id) ","
                              (int-to-string probeset-id) ","
                              (int-to-string probesetfreeze-id)))
-  (assert (= 235 (length bxd-inds)))
   (if name
       (display (string-append "WE HAVE OUR " population " DATASET " name " and trait " trait-name " for precompute!\n")))
   (display ids)
@@ -31,24 +52,6 @@
          (pheno-fn (string-append tmpdir "/pheno.txt"))
          (k-json-fn (string-append tmpdir "/K.json"))
          (gwa-json-fn (string-append tmpdir "/GWA.json"))]
-    (call-with-output-file pheno-fn
-      (lambda (port)
-        (for-each (lambda (ind)
-                    (begin
-                      (let* [(value (assoc-ref traits ind))
-                             (outvalue (if value
-                                           value
-                                           "NA"))]
-                        (display outvalue)
-                        (newline)
-                        (display outvalue port)
-                        (newline port))))
-                  bxd-inds)
-        (close port)
-       ))
-
-    ;; set up with ./.guix-shell -- guile -L . -s ./scripts/precompute/precompute-hits.scm
-
     ;; ---- to start GEMMA precompute inside container
     ;; env TMPDIR=. LD_LIBRARY_PATH=$GUIX_ENVIRONMENT/lib/ guile -L . -s ./scripts/precompute/precompute-hits.scm
     ;; --- First we compute K - control output goes to K.json
diff --git a/scripts/precompute/precompute-hits.scm b/scripts/precompute/precompute-hits.scm
index a18b33d..203dfbb 100755
--- a/scripts/precompute/precompute-hits.scm
+++ b/scripts/precompute/precompute-hits.scm
@@ -1,8 +1,8 @@
 #! Precompute hits
 
-This is the obsolete script to compute GEMMA directly from the DB. We have split out
-this step into fetching traits and running gemma without DB. This allows scaling up
-on HPC etc.
+This is the *obsolete* script to compute GEMMA directly from the DB. We have split out
+this step into separate steps for (1) fetching traits and (2) running gemma without DB.
+This allows scaling up on HPC etc.
 
 Run from base dir with
 
diff --git a/scripts/precompute/run-gemma.scm b/scripts/precompute/run-gemma.scm
index b14927a..36edf8b 100755
--- a/scripts/precompute/run-gemma.scm
+++ b/scripts/precompute/run-gemma.scm
@@ -5,9 +5,9 @@ GEMMA on those.
 
 Run from base dir with
 
-guix shell guile guile-dbi -- guile -L ../.. -s run-gemma.scm --help
+ ~/opt/guix-pull/bin/guix shell guile guile-dbi guile-json -- guile -L . -e main -s scripts/precompute/run-gemma.scm 115476.json
 
-and with some extra paths
+and with some extra paths (for gemma)
 
 . .guix-shell ruby --expose=/home/wrk/iwrk/opensource/code/genetics/gemma-wrapper/=/gemma-wrapper --expose=/home/wrk/iwrk/opensource/code/genetics/gemma/=/gemma -- env TMPDIR=tmp guile -L . -s ./scripts/precompute/run-gemma.scm
 
@@ -33,10 +33,21 @@ and with some extra paths
                          (help    (single-char #\h) (value #f))))
          (options (getopt-long args option-spec))
          (help-wanted (option-ref options 'help #f))]
+    (display "RUNNING")
     (if help-wanted
         (format #t "list-traits-to-compute writes JSON traits files from the GN DB
 Usage: list-traits-to-compute [options...]
   -h, --help          Display this help
-")
-  (display "gemma-run")
-)))
+"))
+    (let [(trait-name "115475")]
+      (call-with-input-file "115475.json"
+        (lambda (port)
+          (let* [(json (json->scm port))
+                 (dataset (assoc-ref json "traits"))
+                 (dataset-name (car (car dataset)))
+                 (traits (assoc-ref dataset dataset-name))
+                 ]
+            (display dataset)
+            (write-pheno-file "pheno.txt" traits)
+                                        ; (run-gemma "BXD" data-id probeset-id probesetfreeze-id name trait-name traits)))
+          ))))))