Run gemma

4 files changed, 31 insertions, 14 deletions

diff --git a/.gitignore b/.gitignore
index 047b170..fb840e4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,2 @@
 BXD.*
+pheno.txt
diff --git a/.guix-shell b/.guix-shell
index 988e529..f3c4ebc 100755
--- a/.guix-shell
+++ b/.guix-shell
@@ -4,4 +4,4 @@
 
 echo "Create a shell to run tools. In the container"
 
-guix shell -C -D -F --network coreutils guile guile-dbi guile-dbd-mysql guile-fibers guile-json guile-gnutls guile-readline guile-redis openssl nss-certs gemma $*
+guix shell -C -D -F --network coreutils guile guile-dbi guile-dbd-mysql guile-fibers guile-json guile-gnutls guile-readline guile-redis openssl nss-certs gemma parallel $*
diff --git a/gn/runner/gemma.scm b/gn/runner/gemma.scm
index 67dae54..ff70f57 100644
--- a/gn/runner/gemma.scm
+++ b/gn/runner/gemma.scm
@@ -23,14 +23,22 @@
   ;; ---- write phenotype file
   (call-with-output-file "pheno.txt"
     (lambda (port)
-      (write 12 port)
-      (newline port)
-      (write "HELLO" 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)))
+
   ;; set up with ./.guix-shell -- guile -L . -s ./scripts/precompute/precompute-hits.scm
 
   ;; ---- to start GEMMA precompute inside container
-  ;; env LD_LIBRARY_PATH=$GUIX_ENVIRONMENT/lib/ guile -L . -s ./scripts/precompute/precompute-hits.scm
+  ;; env TMPDIR=. LD_LIBRARY_PATH=$GUIX_ENVIRONMENT/lib/ guile -L . -s ./scripts/precompute/precompute-hits.scm
   ;; --- First we compute K
-  (system (string-append "env GEMMA_COMMAND=/gemma/bin/gemma /gemma-wrapper/bin/gemma-wrapper --debug -- -gk -g BXD.8_geno.txt.gz -p pheno.txt -a BXD.8_snps.txt"  ))
+  (system (string-append "env GEMMA_COMMAND=gemma /gemma-wrapper/bin/gemma-wrapper --verbose --loco --json --debug --parallel -- -gk -g BXD.8_geno.txt.gz -p pheno.txt -a BXD.8_snps.txt"  ))
   )
diff --git a/scripts/precompute/precompute-hits.scm b/scripts/precompute/precompute-hits.scm
index b9dd746..8762a6a 100755
--- a/scripts/precompute/precompute-hits.scm
+++ b/scripts/precompute/precompute-hits.scm
@@ -11,7 +11,7 @@
              (gn data strains)
              (gn util convert)
              (gn runner gemma)
-             (rnrs base)
+             ; (rnrs base)
              (ice-9 match)
              (srfi srfi-1)
              )
@@ -64,26 +64,34 @@
               (probeset-id (assoc-ref hit "ProbeSetId"))
               (trait (get-trait-name db probeset-id))
               (trait-name (assoc-ref trait "Name"))
+              (name (dataset-name db probesetfreeze-id))
               ]
          (display hit)
          (display data-id)
          (newline)
-         (define name (dataset-name db probesetfreeze-id))
          ;; ---- Get strains and phenotypes for this dataset
          (dbi-query db (string-append "SELECT StrainId,value from ProbeSetData WHERE Id=" data-id-str))
-         (define traits (get-rows-apply db
+         (define id_traits (get-rows-apply db
                                         (lambda (r) `(,(assoc-ref r "StrainId") . ,(assoc-ref r "value")))
                                         '()))
          ;; ---- Now we need to make sure that all strains belong to BXD
          (define non-bxd (fold
                           (lambda (strain lst)
-                            (let ([id (car strain)])
-                              (if (assoc id bxd-strains)
+                            (let* [(id (car strain))
+                                   (name (assoc id bxd-strains))]
+                              (if name
                                   lst
-                                  (append lst `(,id)))))
+                                  (append lst `(,name)))))
 
                           '()
-                          traits))
+                          id_traits))
+         (define traits (map
+                         (lambda (t)
+                           (match t
+                             ((id . value) (cons (assoc-ref bxd-strains id) value)
+                             )))
+                         id_traits))
+         (display traits)
          ;; (if (= 0 (length non-bxd))
          (if (eq? non-bxd '())
              (begin