Replace missing values with MAF

3 files changed, 42 insertions, 18 deletions

diff --git a/wqflask/wqflask/my_pylmm/pyLMM/genotype.py b/wqflask/wqflask/my_pylmm/pyLMM/genotype.py
index 19b0957b..e2457f6b 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/genotype.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/genotype.py
@@ -17,20 +17,35 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 import numpy as np
+from collections import Counter
 
-def normalizeGenotype(g):
+def replace_missing_with_MAF(snp_g):
+    """
+    Replace the missing genotype with the minor allele frequency (MAF)
+    in the snp row
+    """
+    g1 = np.copy(snp_g)
+    cnt = Counter(g1)
+    print cnt
+    min_val = min(cnt.itervalues())
+    print "min_val=",min_val
+    l = [k for k, v in cnt.iteritems() if v == min_val and not np.isnan(k)]
+    print "l=",l[0]
+    return [l[0] if np.isnan(snp) else snp for snp in g1] 
+    
+def normalize(ind_g):
     """
     Run for every SNP list (for one individual) and return
     normalized SNP genotype values with missing data filled in
     """
-    g1 = np.copy(g)          # avoid side effects
-    x = True - np.isnan(g)   # Matrix of True/False
-    m = g[x].mean()          # Global mean value
-    s = np.sqrt(g[x].var())  # Global stddev
-    g1[np.isnan(g)] = m      # Plug-in mean values for missing data
+    g1 = np.copy(ind_g)          # avoid side effects
+    x = True - np.isnan(ind_g)   # Matrix of True/False
+    m = ind_g[x].mean()          # Global mean value
+    s = np.sqrt(ind_g[x].var())  # Global stddev
+    g1[np.isnan(ind_g)] = m      # Plug-in mean values for missing data
     if s == 0:
-        g1 = g1 - m          # Subtract the mean
+        g1 = g1 - m              # Subtract the mean
     else:
-        g1 = (g1 - m) / s    # Normalize the deviation
+        g1 = (g1 - m) / s        # Normalize the deviation
     return g1
 
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py b/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py
index c22da761..bb620052 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py
@@ -19,7 +19,7 @@
 import sys
 import numpy as np
 
-def removeMissingPhenotypes(y,g,verbose=False):
+def remove_missing(y,g,verbose=False):
     """
     Remove missing data from matrices, make sure the genotype data has
     individuals as rows
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
index 35f6e9a9..ffe25fcf 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
@@ -23,8 +23,8 @@ import tsvreader
 import numpy as np
 from lmm import gn2_load_redis, calculate_kinship
 from kinship import kinship, kinship_full
-from genotype import normalizeGenotype
-from phenotype import removeMissingPhenotypes
+import genotype
+import phenotype
 
 usage = """
 python runlmm.py [options] command
@@ -51,6 +51,9 @@ parser.add_option("--pheno",dest="pheno",
                   help="Phenotype file format 1.0")
 parser.add_option("--geno",dest="geno",
                   help="Genotype file format 1.0")
+parser.add_option("--maf-normalization",
+                  action="store_true", dest="maf_normalization", default=False,
+                  help="Apply MAF genotype normalization")
 parser.add_option("--skip-genotype-normalization",
                   action="store_true", dest="skip_genotype_normalization", default=False,
                   help="Skip genotype normalization")
@@ -97,10 +100,10 @@ if cmd == 'redis':
     # Emulating the redis setup of GN2
     gn = []
     for ind_g in g:
-        gn.append( normalizeGenotype(ind_g) )
+        gn.append( genotype.normalize(ind_g) )
     gnt = np.array(gn).T
     if y:
-        Y,G = removeMissingPhenotypes(y,gnt,options.verbose)
+        Y,G = phenotype.remove_missing(y,gnt,options.verbose)
         print "G",G.shape,G
     else:
         G = gnt
@@ -111,14 +114,20 @@ if cmd == 'redis':
     print round(ps[-1],4)
     assert(options.testing and round(ps[-1],4)==0.3461)
 elif cmd == 'kinship':
-    gn = []
+    G = g
+    print G.shape, "\n", G
+    if options.maf_normalization:
+        g1 = np.apply_along_axis( genotype.replace_missing_with_MAF, axis=0, arr=g )
+        print "MAF: ",g1
+    sys.exit()
     for ind_g in g:
         if len(gn)>=8000: break
         if options.skip_genotype_normalization:
-          gn.append(ind_g)
+            gn.append(ind_g)
         else:
-            gn.append( normalizeGenotype(ind_g) )
-    G = np.array(gn)
+            gn.append( genotype.normalize(ind_g) )
+        G = np.array(gn)
+            
     print G.shape, "\n", G
     K = kinship_full(G,options)
     print "first Kinship method",K.shape,"\n",K
@@ -128,7 +137,7 @@ elif cmd == 'kinship':
     print "third Kinship method",K3.shape,"\n",K3
     sys.exit(1)
     gnt = np.array(gn).T
-    Y,g = removeMissingPhenotypes(y,gnt,options.verbose)
+    Y,g = remove_missing_phenotypes(y,gnt,options.verbose)
     G = g
     print G.shape,G
     K = calculate_kinship(np.copy(G),temp_data=None,is_testing=options.testing)