Default to lmm 9 and output lmdb values

genetics-statistics · Jan 7, 2024 · e6d3a8b · e6d3a8b
1 parent fd7c5aa
commit e6d3a8b
Show file tree

Hide file tree

Showing 3 changed files with 11 additions and 5 deletions.
diff --git a/bin/gemma-values.py b/bin/gemma-values.py
@@ -16,6 +16,7 @@
 import argparse
 import json
 import lmdb
+import math
 from struct import *
 
 parser = argparse.ArgumentParser(description="Fetch GEMMA lmdb values.")
@@ -51,4 +52,9 @@
                         chr = "Y"
                     else:
                         chr = str(chr2)
-                    print(chr,pos)
+                    rec = txn.get(key)
+                    af,beta,se,l_mle,p_lrt = unpack('=fffff',rec)
+                    effect = -beta/2.0
+                    minusLogP = -math.log(p_lrt,10)
+                    print(",".join([chr,str(pos),str(round(af,4)),str(round(effect,4)),str(round(se,4)),str(l_mle),str(round(minusLogP,2))]))
+                    sys.exit(2)
diff --git a/bin/gemma-wrapper b/bin/gemma-wrapper
@@ -32,7 +32,7 @@ GEMMA wrapper example:
         -p test/data/input/BXD_pheno.txt \
         -c test/data/input/BXD_covariates2.txt \
         -a test/data/input/BXD_snps.txt \
-        -lmm 2 -maf 0.1 \
+        -lmm 9 -maf 0.1 \
         -debug > GWA.json
 
 Gemma gets used from the path. You can override by setting

diff --git a/bin/gemma2lmdb.py b/bin/gemma2lmdb.py
@@ -43,7 +43,7 @@
             with open(fn) as f:
                 with env.begin(write=True) as txn:
                     for line in f.readlines():
-                        chr,rs,pos,miss,a1,a0,af,logl_H1,l_mle,p_lrt = line.rstrip('\n').split('\t')
+                        chr,rs,pos,miss,a1,a0,af,beta,se,l_mle,p_lrt = line.rstrip('\n').split('\t')
                         if chr=='chr':
                             continue
                         if (chr =='X'):
@@ -62,8 +62,8 @@
                         assert test_pos == int(pos)
                         test_chr = unpack('c',chr_c)
                         # assert test_chr == int(chr), f"{test_chr} vs {int(chr)} - {chr}"
-                        val = pack('=ffff', float(af), float(logl_H1), float(l_mle), float(p_lrt))
-                        assert len(val)==16, f"Packed size is expected to be 16, but is {len(val)}"
+                        val = pack('=fffff', float(af), float(beta), float(se), float(l_mle), float(p_lrt))
+                        assert len(val)==20, f"Packed size is expected to be 20, but is {len(val)}"
                         res = txn.put(key, bytes(val), dupdata=False, overwrite=False)
                         if res == 0:
                             print(f"WARNING: failed to update lmdb record with key {key} -- probably a duplicate {chr}:{pos} ({test_chr_c}:{test_pos})")