aboutsummaryrefslogtreecommitdiff
path: root/wqflask/wqflask
diff options
context:
space:
mode:
authorzsloan2015-05-12 16:50:18 +0000
committerzsloan2015-05-12 16:50:18 +0000
commita544231453be26339d9b1c45888f1e7fffa1ca2e (patch)
tree1eeb20db35015c43cae291b8150bf6712fb8a90a /wqflask/wqflask
parent12b4663b495fed487d52c76656b64ca08273e9a9 (diff)
parent35a81d6a96df1d2fbb315d65af4f7f864b1e5b76 (diff)
downloadgenenetwork2-a544231453be26339d9b1c45888f1e7fffa1ca2e.tar.gz
Resolved conflict before commit
Diffstat (limited to 'wqflask/wqflask')
-rw-r--r--wqflask/wqflask/heatmap/heatmap.py635
-rwxr-xr-xwqflask/wqflask/marker_regression/marker_regression.py6
-rw-r--r--wqflask/wqflask/marker_regression/marker_regression_old.py576
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/__init__.py1
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/benchmark.py44
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/chunks.py96
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/convertlmm.py184
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/genotype.py51
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/gn2.py110
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/gwas.py165
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/input.py267
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/kinship.py168
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/lmm.py995
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/lmm2.py433
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/optmatrix.py55
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/phenotype.py65
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/plink.py107
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/runlmm.py229
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/standalone.py110
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/temp_data.py25
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py122
-rwxr-xr-xwqflask/wqflask/templates/collections/view.html3
22 files changed, 324 insertions, 4123 deletions
diff --git a/wqflask/wqflask/heatmap/heatmap.py b/wqflask/wqflask/heatmap/heatmap.py
index 9b6b1b69..035736fd 100644
--- a/wqflask/wqflask/heatmap/heatmap.py
+++ b/wqflask/wqflask/heatmap/heatmap.py
@@ -1,317 +1,318 @@
-from __future__ import absolute_import, print_function, division
-
-import sys
-sys.path.append(".")
-
-import gc
-import string
-import cPickle
-import os
-import datetime
-import time
-import pp
-import math
-import collections
-import resource
-
-import scipy
-import numpy as np
-from scipy import linalg
-
-from pprint import pformat as pf
-
-from htmlgen import HTMLgen2 as HT
-import reaper
-
-from base.trait import GeneralTrait
-from base import data_set
-from base import species
-from base import webqtlConfig
-from utility import webqtlUtil
-from wqflask.my_pylmm.data import prep_data
-from wqflask.my_pylmm.pyLMM import lmm
-from wqflask.my_pylmm.pyLMM import input
-from utility import helper_functions
-from utility import Plot, Bunch
-from utility import temp_data
-
-from MySQLdb import escape_string as escape
-
-import cPickle as pickle
-import simplejson as json
-
-from pprint import pformat as pf
-
-from redis import Redis
-Redis = Redis()
-
-from flask import Flask, g
-
-class Heatmap(object):
-
- def __init__(self, start_vars, temp_uuid):
-
- trait_db_list = [trait.strip() for trait in start_vars['trait_list'].split(',')]
-
- helper_functions.get_trait_db_obs(self, trait_db_list)
-
- self.temp_uuid = temp_uuid
- self.num_permutations = 5000
- self.dataset = self.trait_list[0][1]
-
- self.json_data = {} #The dictionary that will be used to create the json object that contains all the data needed to create the figure
-
- self.all_sample_list = []
- self.traits = []
-
- chrnames = []
- self.species = species.TheSpecies(dataset=self.trait_list[0][1])
- for key in self.species.chromosomes.chromosomes.keys():
- chrnames.append([self.species.chromosomes.chromosomes[key].name, self.species.chromosomes.chromosomes[key].mb_length])
-
- for trait_db in self.trait_list:
-
- this_trait = trait_db[0]
- self.traits.append(this_trait.name)
- this_sample_data = this_trait.data
-
- for sample in this_sample_data:
- if sample not in self.all_sample_list:
- self.all_sample_list.append(sample)
-
- self.sample_data = []
- for trait_db in self.trait_list:
- this_trait = trait_db[0]
- this_sample_data = this_trait.data
-
- #self.sample_data[this_trait.name] = []
- this_trait_vals = []
- for sample in self.all_sample_list:
- if sample in this_sample_data:
- this_trait_vals.append(this_sample_data[sample].value)
- #self.sample_data[this_trait.name].append(this_sample_data[sample].value)
- else:
- this_trait_vals.append('')
- #self.sample_data[this_trait.name].append('')
- self.sample_data.append(this_trait_vals)
-
- self.gen_reaper_results()
- #self.gen_pylmm_results()
-
- #chrnames = []
- lodnames = []
- chr_pos = []
- pos = []
- markernames = []
-
- for trait in self.trait_results.keys():
- lodnames.append(trait)
-
- for marker in self.dataset.group.markers.markers:
- #if marker['chr'] not in chrnames:
- # chr_ob = [marker['chr'], "filler"]
- # chrnames.append(chr_ob)
- chr_pos.append(marker['chr'])
- pos.append(marker['Mb'])
- markernames.append(marker['name'])
-
- self.json_data['chrnames'] = chrnames
- self.json_data['lodnames'] = lodnames
- self.json_data['chr'] = chr_pos
- self.json_data['pos'] = pos
- self.json_data['markernames'] = markernames
-
- for trait in self.trait_results:
- self.json_data[trait] = self.trait_results[trait]
-
- self.js_data = dict(
- json_data = self.json_data
- )
-
- print("self.js_data:", self.js_data)
-
-
- def gen_reaper_results(self):
- self.trait_results = {}
- for trait_db in self.trait_list:
- self.dataset.group.get_markers()
- this_trait = trait_db[0]
- #this_db = trait_db[1]
- genotype = self.dataset.group.read_genotype_file()
- samples, values, variances = this_trait.export_informative()
-
- trimmed_samples = []
- trimmed_values = []
- for i in range(0, len(samples)):
- if samples[i] in self.dataset.group.samplelist:
- trimmed_samples.append(samples[i])
- trimmed_values.append(values[i])
-
- self.lrs_array = genotype.permutation(strains = trimmed_samples,
- trait = trimmed_values,
- nperm= self.num_permutations)
-
- #self.suggestive = self.lrs_array[int(self.num_permutations*0.37-1)]
- #self.significant = self.lrs_array[int(self.num_permutations*0.95-1)]
-
- reaper_results = genotype.regression(strains = trimmed_samples,
- trait = trimmed_values)
-
-
- lrs_values = [float(qtl.lrs) for qtl in reaper_results]
- print("lrs_values:", lrs_values)
- #self.dataset.group.markers.add_pvalues(p_values)
-
- self.trait_results[this_trait.name] = []
- for qtl in reaper_results:
- if qtl.additive > 0:
- self.trait_results[this_trait.name].append(-float(qtl.lrs))
- else:
- self.trait_results[this_trait.name].append(float(qtl.lrs))
- #for lrs in lrs_values:
- # if
- # self.trait_results[this_trait.name].append(lrs)
-
-
- #this_db_samples = self.dataset.group.samplelist
- #this_sample_data = this_trait.data
- ##print("this_sample_data", this_sample_data)
- #this_trait_vals = []
- #for index, sample in enumerate(this_db_samples):
- # if sample in this_sample_data:
- # sample_value = this_sample_data[sample].value
- # this_trait_vals.append(sample_value)
- # else:
- # this_trait_vals.append("x")
-
- #pheno_vector = np.array([val == "x" and np.nan or float(val) for val in this_trait_vals])
-
- #key = "pylmm:input:" + str(self.temp_uuid)
- #print("key is:", pf(key))
-
- #genotype_data = [marker['genotypes'] for marker in self.dataset.group.markers.markers]
-
- #no_val_samples = self.identify_empty_samples(this_trait_vals)
- #trimmed_genotype_data = self.trim_genotypes(genotype_data, no_val_samples)
-
- #genotype_matrix = np.array(trimmed_genotype_data).T
-
- #print("genotype_matrix:", str(genotype_matrix.tolist()))
- #print("pheno_vector:", str(pheno_vector.tolist()))
-
- #params = dict(pheno_vector = pheno_vector.tolist(),
- # genotype_matrix = genotype_matrix.tolist(),
- # restricted_max_likelihood = True,
- # refit = False,
- # temp_uuid = str(self.temp_uuid),
- #
- # # meta data
- # timestamp = datetime.datetime.now().isoformat(),
- # )
- #
- #json_params = json.dumps(params)
- ##print("json_params:", json_params)
- #Redis.set(key, json_params)
- #Redis.expire(key, 60*60)
- #print("before printing command")
- #
- #command = 'python /home/zas1024/gene/wqflask/wqflask/my_pylmm/pyLMM/lmm.py --key {} --species {}'.format(key,
- # "other")
- #print("command is:", command)
- #print("after printing command")
- #
- #os.system(command)
- #
- #json_results = Redis.blpop("pylmm:results:" + str(self.temp_uuid), 45*60)
-
- def gen_pylmm_results(self):
- self.trait_results = {}
- for trait_db in self.trait_list:
- this_trait = trait_db[0]
- #this_db = trait_db[1]
- self.dataset.group.get_markers()
-
- this_db_samples = self.dataset.group.samplelist
- this_sample_data = this_trait.data
- #print("this_sample_data", this_sample_data)
- this_trait_vals = []
- for index, sample in enumerate(this_db_samples):
- if sample in this_sample_data:
- sample_value = this_sample_data[sample].value
- this_trait_vals.append(sample_value)
- else:
- this_trait_vals.append("x")
-
- pheno_vector = np.array([val == "x" and np.nan or float(val) for val in this_trait_vals])
-
- key = "pylmm:input:" + str(self.temp_uuid)
- #print("key is:", pf(key))
-
- genotype_data = [marker['genotypes'] for marker in self.dataset.group.markers.markers]
-
- no_val_samples = self.identify_empty_samples(this_trait_vals)
- trimmed_genotype_data = self.trim_genotypes(genotype_data, no_val_samples)
-
- genotype_matrix = np.array(trimmed_genotype_data).T
-
- #print("genotype_matrix:", str(genotype_matrix.tolist()))
- #print("pheno_vector:", str(pheno_vector.tolist()))
-
- params = dict(pheno_vector = pheno_vector.tolist(),
- genotype_matrix = genotype_matrix.tolist(),
- restricted_max_likelihood = True,
- refit = False,
- temp_uuid = str(self.temp_uuid),
-
- # meta data
- timestamp = datetime.datetime.now().isoformat(),
- )
-
- json_params = json.dumps(params)
- #print("json_params:", json_params)
- Redis.set(key, json_params)
- Redis.expire(key, 60*60)
- print("before printing command")
-
- command = 'python /home/zas1024/gene/wqflask/wqflask/my_pylmm/pyLMM/lmm.py --key {} --species {}'.format(key,
- "other")
- print("command is:", command)
- print("after printing command")
-
- os.system(command)
-
- json_results = Redis.blpop("pylmm:results:" + str(self.temp_uuid), 45*60)
- results = json.loads(json_results[1])
- p_values = [float(result) for result in results['p_values']]
- #print("p_values:", p_values)
- self.dataset.group.markers.add_pvalues(p_values)
-
- self.trait_results[this_trait.name] = []
- for marker in self.dataset.group.markers.markers:
- self.trait_results[this_trait.name].append(marker['lod_score'])
-
-
- def identify_empty_samples(self, values):
- no_val_samples = []
- for sample_count, val in enumerate(values):
- if val == "x":
- no_val_samples.append(sample_count)
- return no_val_samples
-
- def trim_genotypes(self, genotype_data, no_value_samples):
- trimmed_genotype_data = []
- for marker in genotype_data:
- new_genotypes = []
- for item_count, genotype in enumerate(marker):
- if item_count in no_value_samples:
- continue
- try:
- genotype = float(genotype)
- except ValueError:
- genotype = np.nan
- pass
- new_genotypes.append(genotype)
- trimmed_genotype_data.append(new_genotypes)
- return trimmed_genotype_data
-
- \ No newline at end of file
+from __future__ import absolute_import, print_function, division
+
+import sys
+sys.path.append(".")
+
+import gc
+import string
+import cPickle
+import os
+import datetime
+import time
+import pp
+import math
+import collections
+import resource
+
+import scipy
+import numpy as np
+from scipy import linalg
+
+from pprint import pformat as pf
+
+from htmlgen import HTMLgen2 as HT
+import reaper
+
+from base.trait import GeneralTrait
+from base import data_set
+from base import species
+from base import webqtlConfig
+from utility import webqtlUtil
+from wqflask.my_pylmm.data import prep_data
+# from wqflask.my_pylmm.pyLMM import lmm
+# from wqflask.my_pylmm.pyLMM import input
+from utility import helper_functions
+from utility import Plot, Bunch
+from utility import temp_data
+
+from MySQLdb import escape_string as escape
+
+import cPickle as pickle
+import simplejson as json
+
+from pprint import pformat as pf
+
+from redis import Redis
+Redis = Redis()
+
+from flask import Flask, g
+
+class Heatmap(object):
+
+ def __init__(self, start_vars, temp_uuid):
+
+ trait_db_list = [trait.strip() for trait in start_vars['trait_list'].split(',')]
+
+ helper_functions.get_trait_db_obs(self, trait_db_list)
+
+ self.temp_uuid = temp_uuid
+ self.num_permutations = 5000
+ self.dataset = self.trait_list[0][1]
+
+ self.json_data = {} #The dictionary that will be used to create the json object that contains all the data needed to create the figure
+
+ self.all_sample_list = []
+ self.traits = []
+
+ chrnames = []
+ self.species = species.TheSpecies(dataset=self.trait_list[0][1])
+ for key in self.species.chromosomes.chromosomes.keys():
+ chrnames.append([self.species.chromosomes.chromosomes[key].name, self.species.chromosomes.chromosomes[key].mb_length])
+
+ for trait_db in self.trait_list:
+
+ this_trait = trait_db[0]
+ self.traits.append(this_trait.name)
+ this_sample_data = this_trait.data
+
+ for sample in this_sample_data:
+ if sample not in self.all_sample_list:
+ self.all_sample_list.append(sample)
+
+ self.sample_data = []
+ for trait_db in self.trait_list:
+ this_trait = trait_db[0]
+ this_sample_data = this_trait.data
+
+ #self.sample_data[this_trait.name] = []
+ this_trait_vals = []
+ for sample in self.all_sample_list:
+ if sample in this_sample_data:
+ this_trait_vals.append(this_sample_data[sample].value)
+ #self.sample_data[this_trait.name].append(this_sample_data[sample].value)
+ else:
+ this_trait_vals.append('')
+ #self.sample_data[this_trait.name].append('')
+ self.sample_data.append(this_trait_vals)
+
+ self.gen_reaper_results()
+ #self.gen_pylmm_results()
+
+ #chrnames = []
+ lodnames = []
+ chr_pos = []
+ pos = []
+ markernames = []
+
+ for trait in self.trait_results.keys():
+ lodnames.append(trait)
+
+ for marker in self.dataset.group.markers.markers:
+ #if marker['chr'] not in chrnames:
+ # chr_ob = [marker['chr'], "filler"]
+ # chrnames.append(chr_ob)
+ chr_pos.append(marker['chr'])
+ pos.append(marker['Mb'])
+ markernames.append(marker['name'])
+
+ self.json_data['chrnames'] = chrnames
+ self.json_data['lodnames'] = lodnames
+ self.json_data['chr'] = chr_pos
+ self.json_data['pos'] = pos
+ self.json_data['markernames'] = markernames
+
+ for trait in self.trait_results:
+ self.json_data[trait] = self.trait_results[trait]
+
+ self.js_data = dict(
+ json_data = self.json_data
+ )
+
+ print("self.js_data:", self.js_data)
+
+
+ def gen_reaper_results(self):
+ self.trait_results = {}
+ for trait_db in self.trait_list:
+ self.dataset.group.get_markers()
+ this_trait = trait_db[0]
+ #this_db = trait_db[1]
+ genotype = self.dataset.group.read_genotype_file()
+ samples, values, variances = this_trait.export_informative()
+
+ trimmed_samples = []
+ trimmed_values = []
+ for i in range(0, len(samples)):
+ if samples[i] in self.dataset.group.samplelist:
+ trimmed_samples.append(samples[i])
+ trimmed_values.append(values[i])
+
+ self.lrs_array = genotype.permutation(strains = trimmed_samples,
+ trait = trimmed_values,
+ nperm= self.num_permutations)
+
+ #self.suggestive = self.lrs_array[int(self.num_permutations*0.37-1)]
+ #self.significant = self.lrs_array[int(self.num_permutations*0.95-1)]
+
+ reaper_results = genotype.regression(strains = trimmed_samples,
+ trait = trimmed_values)
+
+
+ lrs_values = [float(qtl.lrs) for qtl in reaper_results]
+ print("lrs_values:", lrs_values)
+ #self.dataset.group.markers.add_pvalues(p_values)
+
+ self.trait_results[this_trait.name] = []
+ for qtl in reaper_results:
+ if qtl.additive > 0:
+ self.trait_results[this_trait.name].append(-float(qtl.lrs))
+ else:
+ self.trait_results[this_trait.name].append(float(qtl.lrs))
+ #for lrs in lrs_values:
+ # if
+ # self.trait_results[this_trait.name].append(lrs)
+
+
+ #this_db_samples = self.dataset.group.samplelist
+ #this_sample_data = this_trait.data
+ ##print("this_sample_data", this_sample_data)
+ #this_trait_vals = []
+ #for index, sample in enumerate(this_db_samples):
+ # if sample in this_sample_data:
+ # sample_value = this_sample_data[sample].value
+ # this_trait_vals.append(sample_value)
+ # else:
+ # this_trait_vals.append("x")
+
+ #pheno_vector = np.array([val == "x" and np.nan or float(val) for val in this_trait_vals])
+
+ #key = "pylmm:input:" + str(self.temp_uuid)
+ #print("key is:", pf(key))
+
+ #genotype_data = [marker['genotypes'] for marker in self.dataset.group.markers.markers]
+
+ #no_val_samples = self.identify_empty_samples(this_trait_vals)
+ #trimmed_genotype_data = self.trim_genotypes(genotype_data, no_val_samples)
+
+ #genotype_matrix = np.array(trimmed_genotype_data).T
+
+ #print("genotype_matrix:", str(genotype_matrix.tolist()))
+ #print("pheno_vector:", str(pheno_vector.tolist()))
+
+ #params = dict(pheno_vector = pheno_vector.tolist(),
+ # genotype_matrix = genotype_matrix.tolist(),
+ # restricted_max_likelihood = True,
+ # refit = False,
+ # temp_uuid = str(self.temp_uuid),
+ #
+ # # meta data
+ # timestamp = datetime.datetime.now().isoformat(),
+ # )
+ #
+ #json_params = json.dumps(params)
+ ##print("json_params:", json_params)
+ #Redis.set(key, json_params)
+ #Redis.expire(key, 60*60)
+ #print("before printing command")
+ #
+ #command = 'python /home/zas1024/gene/wqflask/wqflask/my_pylmm/pyLMM/lmm.py --key {} --species {}'.format(key,
+ # "other")
+ #print("command is:", command)
+ #print("after printing command")
+ #
+ #os.system(command)
+ #
+ #json_results = Redis.blpop("pylmm:results:" + str(self.temp_uuid), 45*60)
+
+ def gen_pylmm_results(self):
+ # This function is NOT used. If it is, we should use a shared function with marker_regression.py
+ self.trait_results = {}
+ for trait_db in self.trait_list:
+ this_trait = trait_db[0]
+ #this_db = trait_db[1]
+ self.dataset.group.get_markers()
+
+ this_db_samples = self.dataset.group.samplelist
+ this_sample_data = this_trait.data
+ #print("this_sample_data", this_sample_data)
+ this_trait_vals = []
+ for index, sample in enumerate(this_db_samples):
+ if sample in this_sample_data:
+ sample_value = this_sample_data[sample].value
+ this_trait_vals.append(sample_value)
+ else:
+ this_trait_vals.append("x")
+
+ pheno_vector = np.array([val == "x" and np.nan or float(val) for val in this_trait_vals])
+
+ key = "pylmm:input:" + str(self.temp_uuid)
+ #print("key is:", pf(key))
+
+ genotype_data = [marker['genotypes'] for marker in self.dataset.group.markers.markers]
+
+ no_val_samples = self.identify_empty_samples(this_trait_vals)
+ trimmed_genotype_data = self.trim_genotypes(genotype_data, no_val_samples)
+
+ genotype_matrix = np.array(trimmed_genotype_data).T
+
+ #print("genotype_matrix:", str(genotype_matrix.tolist()))
+ #print("pheno_vector:", str(pheno_vector.tolist()))
+
+ params = dict(pheno_vector = pheno_vector.tolist(),
+ genotype_matrix = genotype_matrix.tolist(),
+ restricted_max_likelihood = True,
+ refit = False,
+ temp_uuid = str(self.temp_uuid),
+
+ # meta data
+ timestamp = datetime.datetime.now().isoformat(),
+ )
+
+ json_params = json.dumps(params)
+ #print("json_params:", json_params)
+ Redis.set(key, json_params)
+ Redis.expire(key, 60*60)
+ print("before printing command")
+
+ command = 'python /home/zas1024/gene/wqflask/wqflask/my_pylmm/pyLMM/lmm.py --key {} --species {}'.format(key,
+ "other")
+ print("command is:", command)
+ print("after printing command")
+
+ os.system(command)
+
+ json_results = Redis.blpop("pylmm:results:" + str(self.temp_uuid), 45*60)
+ results = json.loads(json_results[1])
+ p_values = [float(result) for result in results['p_values']]
+ #print("p_values:", p_values)
+ self.dataset.group.markers.add_pvalues(p_values)
+
+ self.trait_results[this_trait.name] = []
+ for marker in self.dataset.group.markers.markers:
+ self.trait_results[this_trait.name].append(marker['lod_score'])
+
+
+ def identify_empty_samples(self, values):
+ no_val_samples = []
+ for sample_count, val in enumerate(values):
+ if val == "x":
+ no_val_samples.append(sample_count)
+ return no_val_samples
+
+ def trim_genotypes(self, genotype_data, no_value_samples):
+ trimmed_genotype_data = []
+ for marker in genotype_data:
+ new_genotypes = []
+ for item_count, genotype in enumerate(marker):
+ if item_count in no_value_samples:
+ continue
+ try:
+ genotype = float(genotype)
+ except ValueError:
+ genotype = np.nan
+ pass
+ new_genotypes.append(genotype)
+ trimmed_genotype_data.append(new_genotypes)
+ return trimmed_genotype_data
+
+
diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py
index b3f6d0bb..999a32b8 100755
--- a/wqflask/wqflask/marker_regression/marker_regression.py
+++ b/wqflask/wqflask/marker_regression/marker_regression.py
@@ -37,8 +37,8 @@ from utility import webqtlUtil
from wqflask.marker_regression import gemma_mapping
#from wqflask.marker_regression import rqtl_mapping
from wqflask.my_pylmm.data import prep_data
-from wqflask.my_pylmm.pyLMM import lmm
-from wqflask.my_pylmm.pyLMM import input
+# from wqflask.my_pylmm.pyLMM import lmm
+# from wqflask.my_pylmm.pyLMM import input
from utility import helper_functions
from utility import Plot, Bunch
from utility import temp_data
@@ -51,7 +51,7 @@ if os.environ.get('PYLMM_PATH') is None:
if PYLMM_PATH is None:
PYLMM_PATH=os.environ['HOME']+'/gene/wqflask/wqflask/my_pylmm/pyLMM'
if not os.path.isfile(PYLMM_PATH+'/lmm.py'):
- raise 'PYLMM_PATH unknown or faulty'
+ raise Exception('PYLMM_PATH '+PYLMM_PATH+' unknown or faulty')
PYLMM_COMMAND= 'python '+PYLMM_PATH+'/lmm.py'
class MarkerRegression(object):
diff --git a/wqflask/wqflask/marker_regression/marker_regression_old.py b/wqflask/wqflask/marker_regression/marker_regression_old.py
deleted file mode 100644
index 36331250..00000000
--- a/wqflask/wqflask/marker_regression/marker_regression_old.py
+++ /dev/null
@@ -1,576 +0,0 @@
-from __future__ import absolute_import, print_function, division
-
-from base.trait import GeneralTrait
-from base import data_set #import create_dataset
-
-from pprint import pformat as pf
-
-import string
-import sys
-import datetime
-import os
-import collections
-import uuid
-
-import numpy as np
-from scipy import linalg
-
-import cPickle as pickle
-
-import simplejson as json
-
-from redis import Redis
-Redis = Redis()
-
-from flask import Flask, g
-
-from base.trait import GeneralTrait
-from base import data_set
-from base import species
-from base import webqtlConfig
-from utility import webqtlUtil
-from wqflask.my_pylmm.data import prep_data
-from wqflask.my_pylmm.pyLMM import lmm
-from wqflask.my_pylmm.pyLMM import input
-from utility import helper_functions
-from utility import Plot, Bunch
-from utility import temp_data
-
-from utility.benchmark import Bench
-
-
-class MarkerRegression(object):
-
- def __init__(self, start_vars, temp_uuid):
-
- helper_functions.get_species_dataset_trait(self, start_vars)
-
- #tempdata = temp_data.TempData(temp_uuid)
-
- self.samples = [] # Want only ones with values
- self.vals = []
-
- for sample in self.dataset.group.samplelist:
- value = start_vars['value:' + sample]
- self.samples.append(str(sample))
- self.vals.append(value)
-
- self.mapping_method = start_vars['method']
- self.maf = start_vars['maf'] # Minor allele frequency
- print("self.maf:", self.maf)
-
- self.dataset.group.get_markers()
- if self.mapping_method == "gemma":
- qtl_results = self.run_gemma()
- elif self.mapping_method == "plink":
- qtl_results = self.run_plink()
- #print("qtl_results:", pf(qtl_results))
- elif self.mapping_method == "pylmm":
- print("RUNNING PYLMM")
- #self.qtl_results = self.gen_data(tempdata)
- qtl_results = self.gen_data(str(temp_uuid))
- else:
- print("RUNNING NOTHING")
-
- self.lod_cutoff = 2
- self.filtered_markers = []
- for marker in qtl_results:
- if marker['chr'] > 0:
- self.filtered_markers.append(marker)
-
- #Get chromosome lengths for drawing the manhattan plot
- chromosome_mb_lengths = {}
- for key in self.species.chromosomes.chromosomes.keys():
- chromosome_mb_lengths[key] = self.species.chromosomes.chromosomes[key].mb_length
-
- self.js_data = dict(
- this_trait = self.this_trait.name,
- data_set = self.dataset.name,
- maf = self.maf,
- chromosomes = chromosome_mb_lengths,
- qtl_results = self.filtered_markers,
- )
-
- def run_gemma(self):
- """Generates p-values for each marker using GEMMA"""
-
- #filename = webqtlUtil.genRandStr("{}_{}_".format(self.dataset.group.name, self.this_trait.name))
- self.gen_pheno_txt_file()
-
- os.chdir("/home/zas1024/gene/web/gemma")
-
- gemma_command = './gemma -bfile %s -k output_%s.cXX.txt -lmm 1 -o %s_output' % (
- self.dataset.group.name,
- self.dataset.group.name,
- self.dataset.group.name)
- print("gemma_command:" + gemma_command)
-
- os.system(gemma_command)
-
- included_markers, p_values = self.parse_gemma_output()
-
- self.dataset.group.get_specified_markers(markers = included_markers)
-
- #for marker in self.dataset.group.markers.markers:
- # if marker['name'] not in included_markers:
- # print("marker:", marker)
- # self.dataset.group.markers.markers.remove(marker)
- # #del self.dataset.group.markers.markers[marker]
-
- self.dataset.group.markers.add_pvalues(p_values)
-
- return self.dataset.group.markers.markers
-
-
- def parse_gemma_output(self):
- included_markers = []
- p_values = []
- with open("/home/zas1024/gene/web/gemma/output/{}_output.assoc.txt".format(self.dataset.group.name)) as output_file:
- for line in output_file:
- if line.startswith("chr"):
- continue
- else:
- included_markers.append(line.split("\t")[1])
- p_values.append(float(line.split("\t")[10]))
- #p_values[line.split("\t")[1]] = float(line.split("\t")[10])
- print("p_values: ", p_values)
- return included_markers, p_values
-
- def gen_pheno_txt_file(self):
- """Generates phenotype file for GEMMA"""
-
- #with open("/home/zas1024/gene/web/gemma/tmp_pheno/{}.txt".format(filename), "w") as outfile:
- # for sample, i in enumerate(self.samples):
- # print("sample:" + str(i))
- # print("self.vals[i]:" + str(self.vals[sample]))
- # outfile.write(str(i) + "\t" + str(self.vals[sample]) + "\n")
-
- with open("/home/zas1024/gene/web/gemma/{}.fam".format(self.dataset.group.name), "w") as outfile:
- for i, sample in enumerate(self.samples):
- outfile.write(str(sample) + " " + str(sample) + " 0 0 0 " + str(self.vals[i]) + "\n")
-
- #def gen_plink_for_gemma(self, filename):
- #
- # make_bed = "/home/zas1024/plink/plink --file /home/zas1024/plink/%s --noweb --no-fid --no-parents --no-sex --no-pheno --pheno %s%s.txt --out %s%s --make-bed" % (webqtlConfig.HTMLPATH,
- # webqtlConfig.HTMLPATH,
- # self.dataset.group.name,
- # webqtlConfig.TMPDIR,
- # filename,
- # webqtlConfig.TMPDIR,
- # filename)
- #
- #
-
- def run_plink(self):
-
- os.chdir("/home/zas1024/plink")
-
- plink_output_filename = webqtlUtil.genRandStr("%s_%s_"%(self.dataset.group.name, self.this_trait.name))
-
- self.gen_pheno_txt_file_plink(pheno_filename = plink_output_filename)
-
- plink_command = './plink --noweb --ped %s.ped --no-fid --no-parents --no-sex --no-pheno --map %s.map --pheno %s/%s.txt --pheno-name %s --maf %s --missing-phenotype -9999 --out %s%s --assoc ' % (self.dataset.group.name, self.dataset.group.name, webqtlConfig.TMPDIR, plink_output_filename, self.this_trait.name, self.maf, webqtlConfig.TMPDIR, plink_output_filename)
-
- os.system(plink_command)
-
- count, p_values = self.parse_plink_output(plink_output_filename)
- #gemma_command = './gemma -bfile %s -k output_%s.cXX.txt -lmm 1 -o %s_output' % (
- # self.dataset.group.name,
- # self.dataset.group.name,
- # self.dataset.group.name)
- #print("gemma_command:" + gemma_command)
- #
- #os.system(gemma_command)
- #
- #included_markers, p_values = self.parse_gemma_output()
- #
- #self.dataset.group.get_specified_markers(markers = included_markers)
-
- #for marker in self.dataset.group.markers.markers:
- # if marker['name'] not in included_markers:
- # print("marker:", marker)
- # self.dataset.group.markers.markers.remove(marker)
- # #del self.dataset.group.markers.markers[marker]
-
- print("p_values:", pf(p_values))
-
- self.dataset.group.markers.add_pvalues(p_values)
-
- return self.dataset.group.markers.markers
-
-
- def gen_pheno_txt_file_plink(self, pheno_filename = ''):
- ped_sample_list = self.get_samples_from_ped_file()
- output_file = open("%s%s.txt" % (webqtlConfig.TMPDIR, pheno_filename), "wb")
- header = 'FID\tIID\t%s\n' % self.this_trait.name
- output_file.write(header)
-
- new_value_list = []
-
- #if valueDict does not include some strain, value will be set to -9999 as missing value
- for i, sample in enumerate(ped_sample_list):
- try:
- value = self.vals[i]
- value = str(value).replace('value=','')
- value = value.strip()
- except:
- value = -9999
-
- new_value_list.append(value)
-
-
- new_line = ''
- for i, sample in enumerate(ped_sample_list):
- j = i+1
- value = new_value_list[i]
- new_line += '%s\t%s\t%s\n'%(sample, sample, value)
-
- if j%1000 == 0:
- output_file.write(newLine)
- new_line = ''
-
- if new_line:
- output_file.write(new_line)
-
- output_file.close()
-
- # get strain name from ped file in order
- def get_samples_from_ped_file(self):
-
- os.chdir("/home/zas1024/plink")
-
- ped_file= open("{}.ped".format(self.dataset.group.name),"r")
- line = ped_file.readline()
- sample_list=[]
-
- while line:
- lineList = string.split(string.strip(line), '\t')
- lineList = map(string.strip, lineList)
-
- sample_name = lineList[0]
- sample_list.append(sample_name)
-
- line = ped_file.readline()
-
- return sample_list
-
- ################################################################
- # Generate Chr list, Chr OrderId and Retrieve Length Information
- ################################################################
- #def getChrNameOrderIdLength(self,RISet=''):
- # try:
- # query = """
- # Select
- # Chr_Length.Name,Chr_Length.OrderId,Length from Chr_Length, InbredSet
- # where
- # Chr_Length.SpeciesId = InbredSet.SpeciesId AND
- # InbredSet.Name = '%s'
- # Order by OrderId
- # """ % (self.dataset.group.name)
- # results =g.db.execute(query).fetchall()
- # ChrList=[]
- # ChrLengthMbList=[]
- # ChrNameOrderIdDict={}
- # ChrOrderIdNameDict={}
- #
- # for item in results:
- # ChrList.append(item[0])
- # ChrNameOrderIdDict[item[0]]=item[1] # key is chr name, value is orderId
- # ChrOrderIdNameDict[item[1]]=item[0] # key is orderId, value is chr name
- # ChrLengthMbList.append(item[2])
- #
- # except:
- # ChrList=[]
- # ChrNameOrderIdDict={}
- # ChrLengthMbList=[]
- #
- # return ChrList,ChrNameOrderIdDict,ChrOrderIdNameDict,ChrLengthMbList
-
-
- def parse_plink_output(self, output_filename):
- plink_results={}
-
- threshold_p_value = 0.01
-
- result_fp = open("%s%s.qassoc"% (webqtlConfig.TMPDIR, output_filename), "rb")
-
- header_line = result_fp.readline()# read header line
- line = result_fp.readline()
-
- value_list = [] # initialize value list, this list will include snp, bp and pvalue info
- p_value_dict = {}
- count = 0
-
- while line:
- #convert line from str to list
- line_list = self.build_line_list(line=line)
-
- # only keep the records whose chromosome name is in db
- if self.species.chromosomes.chromosomes.has_key(int(line_list[0])) and line_list[-1] and line_list[-1].strip()!='NA':
-
- chr_name = self.species.chromosomes.chromosomes[int(line_list[0])]
- snp = line_list[1]
- BP = line_list[2]
- p_value = float(line_list[-1])
- if threshold_p_value >= 0 and threshold_p_value <= 1:
- if p_value < threshold_p_value:
- p_value_dict[snp] = p_value
-
- if plink_results.has_key(chr_name):
- value_list = plink_results[chr_name]
-
- # pvalue range is [0,1]
- if threshold_p_value >=0 and threshold_p_value <= 1:
- if p_value < threshold_p_value:
- value_list.append((snp, BP, p_value))
- count += 1
-
- plink_results[chr_name] = value_list
- value_list = []
- else:
- if threshold_p_value >= 0 and threshold_p_value <= 1:
- if p_value < threshold_p_value:
- value_list.append((snp, BP, p_value))
- count += 1
-
- if value_list:
- plink_results[chr_name] = value_list
-
- value_list=[]
-
- line = result_fp.readline()
- else:
- line = result_fp.readline()
-
- #if p_value_list:
- # min_p_value = min(p_value_list)
- #else:
- # min_p_value = 0
-
- return count, p_value_dict
-
- ######################################################
- # input: line: str,one line read from file
- # function: convert line from str to list;
- # output: lineList list
- #######################################################
- def build_line_list(self, line=None):
-
- line_list = string.split(string.strip(line),' ')# irregular number of whitespaces between columns
- line_list = [item for item in line_list if item <>'']
- line_list = map(string.strip, line_list)
-
- return line_list
-
- #def gen_data(self, tempdata):
- def gen_data(self, temp_uuid):
- """Generates p-values for each marker"""
-
- pheno_vector = np.array([val == "x" and np.nan or float(val) for val in self.vals])
-
- #lmm_uuid = str(uuid.uuid4())
-
- key = "pylmm:input:" + temp_uuid
- print("key is:", pf(key))
- #with Bench("Loading cache"):
- # result = Redis.get(key)
-
- if self.dataset.group.species == "human":
- p_values, t_stats = self.gen_human_results(pheno_vector, key, temp_uuid)
- #p_values = self.trim_results(p_values)
-
- else:
- print("NOW CWD IS:", os.getcwd())
- genotype_data = [marker['genotypes'] for marker in self.dataset.group.markers.markers]
-
- no_val_samples = self.identify_empty_samples()
- trimmed_genotype_data = self.trim_genotypes(genotype_data, no_val_samples)
-
- genotype_matrix = np.array(trimmed_genotype_data).T
-
- #print("pheno_vector: ", pf(pheno_vector))
- #print("genotype_matrix: ", pf(genotype_matrix))
- #print("genotype_matrix.shape: ", pf(genotype_matrix.shape))
-
- #params = {"pheno_vector": pheno_vector,
- # "genotype_matrix": genotype_matrix,
- # "restricted_max_likelihood": True,
- # "refit": False,
- # "temp_data": tempdata}
-
- params = dict(pheno_vector = pheno_vector.tolist(),
- genotype_matrix = genotype_matrix.tolist(),
- restricted_max_likelihood = True,
- refit = False,
- temp_uuid = temp_uuid,
-
- # meta data
- timestamp = datetime.datetime.now().isoformat(),
- )
-
- json_params = json.dumps(params)
- #print("json_params:", json_params)
- Redis.set(key, json_params)
- Redis.expire(key, 60*60)
- print("before printing command")
-
- command = 'python /home/zas1024/gene/wqflask/wqflask/my_pylmm/pyLMM/lmm.py --key {} --species {}'.format(key,
- "other")
- print("command is:", command)
- print("after printing command")
-
- os.system(command)
-
- #t_stats, p_values = lmm.run(key)
- #lmm.run(key)
-
- json_results = Redis.blpop("pylmm:results:" + temp_uuid, 45*60)
- results = json.loads(json_results[1])
- p_values = [float(result) for result in results['p_values']]
- print("p_values:", p_values)
- #p_values = self.trim_results(p_values)
- t_stats = results['t_stats']
-
- #t_stats, p_values = lmm.run(
- # pheno_vector,
- # genotype_matrix,
- # restricted_max_likelihood=True,
- # refit=False,
- # temp_data=tempdata
- #)
- #print("p_values:", p_values)
-
- self.dataset.group.markers.add_pvalues(p_values)
-
- #self.get_lod_score_cutoff()
-
- return self.dataset.group.markers.markers
-
- def trim_results(self, p_values):
- print("len_p_values:", len(p_values))
- if len(p_values) > 500:
- p_values.sort(reverse=True)
- trimmed_values = p_values[:500]
-
- return trimmed_values
-
- #def gen_human_results(self, pheno_vector, tempdata):
- def gen_human_results(self, pheno_vector, key, temp_uuid):
- file_base = os.path.join(webqtlConfig.PYLMM_PATH, self.dataset.group.name)
-
- plink_input = input.plink(file_base, type='b')
- input_file_name = os.path.join(webqtlConfig.SNP_PATH, self.dataset.group.name + ".snps.gz")
-
- pheno_vector = pheno_vector.reshape((len(pheno_vector), 1))
- covariate_matrix = np.ones((pheno_vector.shape[0],1))
- kinship_matrix = np.fromfile(open(file_base + '.kin','r'),sep=" ")
- kinship_matrix.resize((len(plink_input.indivs),len(plink_input.indivs)))
-
- print("Before creating params")
-
- params = dict(pheno_vector = pheno_vector.tolist(),
- covariate_matrix = covariate_matrix.tolist(),
- input_file_name = input_file_name,
- kinship_matrix = kinship_matrix.tolist(),
- refit = False,
- temp_uuid = temp_uuid,
-
- # meta data
- timestamp = datetime.datetime.now().isoformat(),
- )
-
- print("After creating params")
-
- json_params = json.dumps(params)
- Redis.set(key, json_params)
- Redis.expire(key, 60*60)
-
- print("Before creating the command")
-
- command = 'python /home/zas1024/gene/wqflask/wqflask/my_pylmm/pyLMM/lmm.py --key {} --species {}'.format(key,
- "human")
-
- print("command is:", command)
-
- os.system(command)
-
- json_results = Redis.blpop("pylmm:results:" + temp_uuid, 45*60)
- results = json.loads(json_results[1])
- t_stats = results['t_stats']
- p_values = results['p_values']
-
-
- #p_values, t_stats = lmm.run_human(key)
-
- #p_values, t_stats = lmm.run_human(
- # pheno_vector,
- # covariate_matrix,
- # input_file_name,
- # kinship_matrix,
- # loading_progress=tempdata
- # )
-
- return p_values, t_stats
-
- def get_lod_score_cutoff(self):
- print("INSIDE GET LOD CUTOFF")
- high_qtl_count = 0
- for marker in self.dataset.group.markers.markers:
- if marker['lod_score'] > 1:
- high_qtl_count += 1
-
- if high_qtl_count > 1000:
- return 1
- else:
- return 0
-
- def identify_empty_samples(self):
- no_val_samples = []
- for sample_count, val in enumerate(self.vals):
- if val == "x":
- no_val_samples.append(sample_count)
- return no_val_samples
-
- def trim_genotypes(self, genotype_data, no_value_samples):
- trimmed_genotype_data = []
- for marker in genotype_data:
- new_genotypes = []
- for item_count, genotype in enumerate(marker):
- if item_count in no_value_samples:
- continue
- try:
- genotype = float(genotype)
- except ValueError:
- genotype = np.nan
- pass
- new_genotypes.append(genotype)
- trimmed_genotype_data.append(new_genotypes)
- return trimmed_genotype_data
-
-def create_snp_iterator_file(group):
- plink_file_base = os.path.join(webqtlConfig.PYLMM_PATH, group)
- plink_input = input.plink(plink_file_base, type='b')
-
- data = dict(plink_input = list(plink_input),
- numSNPs = plink_input.numSNPs)
-
- #input_dict = {}
- #
- #input_dict['plink_input'] = list(plink_input)
- #input_dict['numSNPs'] = plink_input.numSNPs
- #
-
- snp_file_base = os.path.join(webqtlConfig.SNP_PATH, group + ".snps.gz")
-
- with gzip.open(snp_file_base, "wb") as fh:
- pickle.dump(data, fh, pickle.HIGHEST_PROTOCOL)
-
-#if __name__ == '__main__':
-# import cPickle as pickle
-# import gzip
-# create_snp_iterator_file("HLC")
-
-if __name__ == '__main__':
- import cPickle as pickle
- import gzip
- create_snp_iterator_file("HLC")
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/__init__.py b/wqflask/wqflask/my_pylmm/pyLMM/__init__.py
deleted file mode 100644
index f33c4e74..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-PYLMM_VERSION="0.51-gn2"
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/benchmark.py b/wqflask/wqflask/my_pylmm/pyLMM/benchmark.py
deleted file mode 100644
index 6c6c9f88..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/benchmark.py
+++ /dev/null
@@ -1,44 +0,0 @@
-from __future__ import print_function, division, absolute_import
-
-import collections
-import inspect
-import time
-
-class Bench(object):
- entries = collections.OrderedDict()
-
- def __init__(self, name=None):
- self.name = name
-
- def __enter__(self):
- if self.name:
- print("Starting benchmark: %s" % (self.name))
- else:
- print("Starting benchmark at: %s [%i]" % (inspect.stack()[1][3], inspect.stack()[1][2]))
- self.start_time = time.time()
-
- def __exit__(self, type, value, traceback):
- if self.name:
- name = self.name
- else:
- name = "That"
-
- time_taken = time.time() - self.start_time
- print(" %s took: %f seconds" % (name, (time_taken)))
-
- if self.name:
- Bench.entries[self.name] = Bench.entries.get(self.name, 0) + time_taken
-
-
- @classmethod
- def report(cls):
- total_time = sum((time_taken for time_taken in cls.entries.itervalues()))
- print("\nTiming report\n")
- for name, time_taken in cls.entries.iteritems():
- percent = int(round((time_taken/total_time) * 100))
- print("[{}%] {}: {}".format(percent, name, time_taken))
- print()
-
- def reset(cls):
- """Reset the entries"""
- cls.entries = collections.OrderedDict()
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/chunks.py b/wqflask/wqflask/my_pylmm/pyLMM/chunks.py
deleted file mode 100644
index 9565fb96..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/chunks.py
+++ /dev/null
@@ -1,96 +0,0 @@
-from __future__ import absolute_import, print_function, division
-
-import math
-import time
-
-
-def divide_into_chunks(the_list, number_chunks):
- """Divides a list into approximately number_chunks smaller lists
-
- >>> divide_into_chunks([1, 2, 7, 3, 22, 8, 5, 22, 333], 3)
- [[1, 2, 7], [3, 22, 8], [5, 22, 333]]
- >>> divide_into_chunks([1, 2, 7, 3, 22, 8, 5, 22, 333], 4)
- [[1, 2, 7], [3, 22, 8], [5, 22, 333]]
- >>> divide_into_chunks([1, 2, 7, 3, 22, 8, 5, 22, 333], 5)
- [[1, 2], [7, 3], [22, 8], [5, 22], [333]]
- >>>
-
- """
- length = len(the_list)
-
- if length == 0:
- return [[]]
-
- if length <= number_chunks:
- number_chunks = length
-
- chunksize = int(math.ceil(length / number_chunks))
-
- chunks = []
- for counter in range(0, length, chunksize):
- chunks.append(the_list[counter:counter+chunksize])
-
- return chunks
-
-def _confirm_chunk(original, result):
- all_chunked = []
- for chunk in result:
- all_chunked.extend(chunk)
- print("length of all chunked:", len(all_chunked))
- assert original == all_chunked, "You didn't chunk right"
-
-
-def _chunk_test(divide_func):
- import random
- random.seed(7)
-
- number_exact = 0
- total_amount_off = 0
-
- for test in range(1, 1001):
- print("\n\ntest:", test)
- number_chunks = random.randint(1, 20)
- number_elements = random.randint(0, 100)
- the_list = list(range(1, number_elements))
- result = divide_func(the_list, number_chunks)
-
- print("Dividing list of length {} into approximately {} chunks - got {} chunks".format(
- len(the_list), number_chunks, len(result)))
- print("result:", result)
-
- _confirm_chunk(the_list, result)
-
- amount_off = abs(number_chunks - len(result))
- if amount_off == 0:
- number_exact += 1
- else:
- total_amount_off += amount_off
-
-
- print("\n{} exact out of {} [Total amount off: {}]".format(number_exact,
- test,
- total_amount_off))
- assert number_exact == 558
- assert total_amount_off == 1580
- return number_exact, total_amount_off
-
-
-def _main():
- info = dict()
- #funcs = (("sam", sam_divide_into_chunks), ("zach", zach_divide_into_chunks))
- funcs = (("only one", divide_into_chunks),)
- for name, func in funcs:
- start = time.time()
- number_exact, total_amount_off = _chunk_test(func)
- took = time.time() - start
- info[name] = dict(number_exact=number_exact,
- total_amount_off=total_amount_off,
- took=took)
-
- print("info is:", info)
-
-if __name__ == '__main__':
- _main()
- print("\nConfirming doctests...")
- import doctest
- doctest.testmod()
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/convertlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/convertlmm.py
deleted file mode 100644
index 4312fed0..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/convertlmm.py
+++ /dev/null
@@ -1,184 +0,0 @@
-# This is a converter for common LMM formats, so as to keep file
-# reader complexity outside the main routines.
-
-# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-from __future__ import print_function
-from optparse import OptionParser
-import sys
-import os
-import numpy as np
-# from lmm import LMM, run_other
-# import input
-import plink
-
-usage = """
-python convertlmm.py [--plink] [--prefix out_basename] [--kinship kfile] [--pheno pname] [--geno gname]
-
- Convert files for runlmm.py processing. Writes to stdout by default.
-
- try --help for more information
-
-Examples:
-
- python ./my_pylmm/pyLMM/convertlmm.py --plink --pheno data/test_snps.132k.clean.noX.fake.phenos > test.pheno
-
- python ./my_pylmm/pyLMM/convertlmm.py --plink --pheno data/test_snps.132k.clean.noX.fake.phenos --geno data/test_snps.132k.clean.noX > test.geno
-"""
-
-# if len(args) == 0:
-# print usage
-# sys.exit(1)
-
-option_parser = OptionParser(usage=usage)
-option_parser.add_option("--kinship", dest="kinship",
- help="Parse a kinship file. This is an nxn plain text file and can be computed with the pylmmKinship program")
-option_parser.add_option("--pheno", dest="pheno",
- help="Parse a phenotype file (use with --plink only)")
-option_parser.add_option("--geno", dest="geno",
- help="Parse a genotype file (use with --plink only)")
-option_parser.add_option("--plink", dest="plink", action="store_true", default=False,
- help="Parse PLINK style")
-# option_parser.add_option("--kinship",action="store_false", dest="kinship", default=True,
-# help="Parse a kinship file. This is an nxn plain text file and can be computed with the pylmmKinship program.")
-option_parser.add_option("--prefix", dest="prefix",
- help="Output prefix for output file(s)")
-option_parser.add_option("-q", "--quiet",
- action="store_false", dest="verbose", default=True,
- help="don't print status messages to stdout")
-option_parser.add_option("-v", "--verbose",
- action="store_true", dest="verbose", default=False,
- help="Print extra info")
-
-(options, args) = option_parser.parse_args()
-
-writer = None
-num_inds = None
-snp_names = []
-ind_names = []
-
-def msg(s):
- sys.stderr.write("INFO: ")
- sys.stderr.write(s)
- sys.stderr.write("\n")
-
-def wr(s):
- if writer:
- writer.write(s)
- else:
- sys.stdout.write(s)
-
-def wrln(s):
- wr(s)
- wr("\n")
-
-
-if options.pheno:
- if not options.plink:
- raise Exception("Use --plink switch")
- # Because plink does not track size we need to read the whole thing first
- msg("Converting pheno "+options.pheno)
- phenos = []
- count = 0
- count_pheno = None
- for line in open(options.pheno,'r'):
- count += 1
- list = line.split()
- pcount = len(list)-2
- assert(pcount > 0)
- if count_pheno == None:
- count_pheno = pcount
- assert(count_pheno == pcount)
- row = [list[0]]+list[2:]
- phenos.append(row)
-
- writer = None
- if options.prefix:
- writer = open(options.prefix+".pheno","w")
- wrln("# Phenotype format version 1.0")
- wrln("# Individuals = "+str(count))
- wrln("# Phenotypes = "+str(count_pheno))
- for i in range(count_pheno):
- wr("\t"+str(i+1))
- wr("\n")
- for i in range(count):
- wr("\t".join(phenos[i]))
- wr("\n")
- num_inds = count
- msg(str(count)+" pheno lines written")
-
-if options.kinship:
- is_header = True
- count = 0
- msg("Converting kinship "+options.kinship)
- writer = None
- if options.prefix:
- writer = open(options.prefix+".kin","w")
- for line in open(options.kinship,'r'):
- count += 1
- if is_header:
- size = len(line.split())
- wrln("# Kinship format version 1.0")
- wrln("# Size="+str(size))
- for i in range(size):
- wr("\t"+str(i+1))
- wr("\n")
- is_header = False
- wr(str(count))
- wr("\t")
- wr("\t".join(line.split()))
- wr("\n")
- num_inds = count
- msg(str(count)+" kinship lines written")
-
-if options.geno:
- msg("Converting geno "+options.geno+'.bed')
- if not options.plink:
- raise Exception("Use --plink switch")
- if not num_inds:
- raise Exception("Can not figure out the number of individuals, use --pheno or --kinship")
- bim_snps = plink.readbim(options.geno+'.bim')
- num_snps = len(bim_snps)
- writer = None
- if options.prefix:
- writer = open(options.prefix+".geno","w")
- wrln("# Genotype format version 1.0")
- wrln("# Individuals = "+str(num_inds))
- wrln("# SNPs = "+str(num_snps))
- wrln("# Encoding = HAB")
- for i in range(num_inds):
- wr("\t"+str(i+1))
- wr("\n")
-
- m = []
- def out(i,x):
- # wr(str(i)+"\t")
- # wr("\t".join(x))
- # wr("\n")
- m.append(x)
-
- snps = plink.readbed(options.geno+'.bed',num_inds, ('A','H','B','-'), out)
-
- msg("Write transposed genotype matrix")
- for g in range(num_snps):
- wr(bim_snps[g][1]+"\t")
- for i in range(num_inds):
- wr(m[g][i])
- wr("\n")
-
- msg(str(count)+" geno lines written (with "+str(snps)+" snps)")
-
-msg("Converting done")
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/genotype.py b/wqflask/wqflask/my_pylmm/pyLMM/genotype.py
deleted file mode 100644
index 49f32e3a..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/genotype.py
+++ /dev/null
@@ -1,51 +0,0 @@
-# Genotype routines
-
-# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com)
-# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-import numpy as np
-from collections import Counter
-import operator
-
-def replace_missing_with_MAF(snp_g):
- """
- Replace the missing genotype with the minor allele frequency (MAF)
- in the snp row. It is rather slow!
- """
- cnt = Counter(snp_g)
- tuples = sorted(cnt.items(), key=operator.itemgetter(1))
- l2 = [t for t in tuples if not np.isnan(t[0])]
- maf = l2[0][0]
- res = np.array([maf if np.isnan(snp) else snp for snp in snp_g])
- return res
-
-def normalize(ind_g):
- """
- Run for every SNP list (for one individual) and return
- normalized SNP genotype values with missing data filled in
- """
- g = np.copy(ind_g) # copy to avoid side effects
- missing = np.isnan(g)
- values = g[True - missing]
- mean = values.mean() # Global mean value
- stddev = np.sqrt(values.var()) # Global stddev
- g[missing] = mean # Plug-in mean values for missing data
- if stddev == 0:
- g = g - mean # Subtract the mean
- else:
- g = (g - mean) / stddev # Normalize the deviation
- return g
-
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/gn2.py b/wqflask/wqflask/my_pylmm/pyLMM/gn2.py
deleted file mode 100644
index 821195c8..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/gn2.py
+++ /dev/null
@@ -1,110 +0,0 @@
-# Standalone specific methods and callback handler
-#
-# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# Set the log level with
-#
-# logging.basicConfig(level=logging.DEBUG)
-
-from __future__ import absolute_import, print_function, division
-
-import numpy as np
-import sys
-import logging
-
-# logger = logging.getLogger(__name__)
-logger = logging.getLogger('lmm2')
-logging.basicConfig(level=logging.DEBUG)
-np.set_printoptions(precision=3,suppress=True)
-
-progress_location = None
-progress_current = None
-progress_prev_perc = None
-
-def progress_default_func(location,count,total):
- global progress_current
- value = round(count*100.0/total)
- progress_current = value
-
-progress_func = progress_default_func
-
-def progress_set_func(func):
- global progress_func
- progress_func = func
-
-def progress(location, count, total):
- global progress_location
- global progress_prev_perc
-
- perc = round(count*100.0/total)
- if perc != progress_prev_perc and (location != progress_location or perc > 98 or perc > progress_prev_perc + 5):
- progress_func(location, count, total)
- logger.info("Progress: %s %d%%" % (location,perc))
- progress_location = location
- progress_prev_perc = perc
-
-def mprint(msg,data):
- """
- Array/matrix print function
- """
- m = np.array(data)
- if m.ndim == 1:
- print(msg,m.shape,"=\n",m[0:3]," ... ",m[-3:])
- if m.ndim == 2:
- print(msg,m.shape,"=\n[",
- m[0][0:3]," ... ",m[0][-3:],"\n ",
- m[1][0:3]," ... ",m[1][-3:],"\n ...\n ",
- m[-2][0:3]," ... ",m[-2][-3:],"\n ",
- m[-1][0:3]," ... ",m[-1][-3:],"]")
-
-def fatal(msg):
- logger.critical(msg)
- raise Exception(msg)
-
-def callbacks():
- return dict(
- write = sys.stdout.write,
- writeln = print,
- debug = logger.debug,
- info = logger.info,
- warning = logger.warning,
- error = logger.error,
- critical = logger.critical,
- fatal = fatal,
- progress = progress,
- mprint = mprint
- )
-
-def uses(*funcs):
- """
- Some sugar
- """
- return [callbacks()[func] for func in funcs]
-
-# ----- Minor test cases:
-
-if __name__ == '__main__':
- # logging.basicConfig(level=logging.DEBUG)
- logging.debug("Test %i" % (1))
- d = callbacks()['debug']
- d("TEST")
- wrln = callbacks()['writeln']
- wrln("Hello %i" % 34)
- progress = callbacks()['progress']
- progress("I am half way",50,100)
- list = [0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
- 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
- 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
- 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
- 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15]
- mprint("list",list)
- matrix = [[1,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
- [2,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
- [3,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
- [4,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
- [5,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
- [6,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15]]
- mprint("matrix",matrix)
- ix,dx = uses("info","debug")
- ix("ix")
- dx("dx")
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/gwas.py b/wqflask/wqflask/my_pylmm/pyLMM/gwas.py
deleted file mode 100644
index 247a8729..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/gwas.py
+++ /dev/null
@@ -1,165 +0,0 @@
-# pylmm-based GWAS calculation
-#
-# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com)
-# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-#!/usr/bin/python
-
-import pdb
-import time
-import sys
-import lmm2
-
-import os
-import numpy as np
-import input
-from optmatrix import matrix_initialize
-from lmm2 import LMM2
-
-import multiprocessing as mp # Multiprocessing is part of the Python stdlib
-import Queue
-
-# ---- A trick to decide on the environment:
-try:
- from wqflask.my_pylmm.pyLMM import chunks
- from gn2 import uses
-except ImportError:
- has_gn2=False
- from standalone import uses
-
-progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal')
-
-
-def formatResult(id,beta,betaSD,ts,ps):
- return "\t".join([str(x) for x in [id,beta,betaSD,ts,ps]]) + "\n"
-
-def compute_snp(j,n,snp_ids,lmm2,REML,q = None):
- result = []
- for snp_id in snp_ids:
- snp,id = snp_id
- x = snp.reshape((n,1)) # all the SNPs
- # if refit:
- # L.fit(X=snp,REML=REML)
- ts,ps,beta,betaVar = lmm2.association(x,REML=REML,returnBeta=True)
- # result.append(formatResult(id,beta,np.sqrt(betaVar).sum(),ts,ps))
- result.append( (ts,ps) )
- if not q:
- q = compute_snp.q
- q.put([j,result])
- return j
-
-def f_init(q):
- compute_snp.q = q
-
-def gwas(Y,G,K,restricted_max_likelihood=True,refit=False,verbose=True):
- """
- GWAS. The G matrix should be n inds (cols) x m snps (rows)
- """
- info("In gwas.gwas")
- matrix_initialize()
- cpu_num = mp.cpu_count()
- numThreads = None # for now use all available threads
- kfile2 = False
- reml = restricted_max_likelihood
-
- mprint("G",G)
- n = G.shape[1] # inds
- inds = n
- m = G.shape[0] # snps
- snps = m
- info("%s SNPs",snps)
- assert snps>=inds, "snps should be larger than inds (snps=%d,inds=%d)" % (snps,inds)
-
- # CREATE LMM object for association
- # if not kfile2: L = LMM(Y,K,Kva,Kve,X0,verbose=verbose)
- # else: L = LMM_withK2(Y,K,Kva,Kve,X0,verbose=verbose,K2=K2)
-
- lmm2 = LMM2(Y,K) # ,Kva,Kve,X0,verbose=verbose)
- if not refit:
- info("Computing fit for null model")
- lmm2.fit() # follow GN model in run_other
- info("heritability=%0.3f, sigma=%0.3f" % (lmm2.optH,lmm2.optSigma))
-
- res = []
-
- # Set up the pool
- # mp.set_start_method('spawn')
- q = mp.Queue()
- p = mp.Pool(numThreads, f_init, [q])
- collect = []
-
- count = 0
- job = 0
- jobs_running = 0
- jobs_completed = 0
- for snp in G:
- snp_id = (snp,'SNPID')
- count += 1
- if count % 1000 == 0:
- job += 1
- debug("Job %d At SNP %d" % (job,count))
- if numThreads == 1:
- debug("Running on 1 THREAD")
- compute_snp(job,n,collect,lmm2,reml,q)
- collect = []
- j,lst = q.get()
- debug("Job "+str(j)+" finished")
- jobs_completed += 1
- progress("GWAS2",jobs_completed,snps/1000)
- res.append((j,lst))
- else:
- p.apply_async(compute_snp,(job,n,collect,lmm2,reml))
- jobs_running += 1
- collect = []
- while jobs_running > cpu_num:
- try:
- j,lst = q.get_nowait()
- debug("Job "+str(j)+" finished")
- jobs_completed += 1
- progress("GWAS2",jobs_completed,snps/1000)
- res.append((j,lst))
- jobs_running -= 1
- except Queue.Empty:
- time.sleep(0.1)
- pass
- if jobs_running > cpu_num*2:
- time.sleep(1.0)
- else:
- break
-
- collect.append(snp_id)
-
- if numThreads==1 or count<1000 or len(collect)>0:
- job += 1
- debug("Collect final batch size %i job %i @%i: " % (len(collect), job, count))
- compute_snp(job,n,collect,lmm2,reml,q)
- collect = []
- j,lst = q.get()
- res.append((j,lst))
- debug("count=%i running=%i collect=%i" % (count,jobs_running,len(collect)))
- for job in range(jobs_running):
- j,lst = q.get(True,15) # time out
- debug("Job "+str(j)+" finished")
- jobs_completed += 1
- progress("GWAS2",jobs_completed,snps/1000)
- res.append((j,lst))
-
- mprint("Before sort",[res1[0] for res1 in res])
- res = sorted(res,key=lambda x: x[0])
- mprint("After sort",[res1[0] for res1 in res])
- info([len(res1[1]) for res1 in res])
- ts = [item[0] for j,res1 in res for item in res1]
- ps = [item[1] for j,res1 in res for item in res1]
- return ts,ps
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/input.py b/wqflask/wqflask/my_pylmm/pyLMM/input.py
deleted file mode 100644
index 7063fedf..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/input.py
+++ /dev/null
@@ -1,267 +0,0 @@
-# pylmm is a python-based linear mixed-model solver with applications to GWAS
-
-# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com)
-
-#The program is free for academic use. Please contact Nick Furlotte
-#<nick.furlotte@gmail.com> if you are interested in using the software for
-#commercial purposes.
-
-#The software must not be modified and distributed without prior
-#permission of the author.
-
-#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-#"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-#LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
-#A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
-#CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
-#EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
-#PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
-#PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
-#LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
-#NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-#SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
-import os
-import sys
-import numpy as np
-import struct
-import pdb
-
-class plink:
- def __init__(self,fbase,kFile=None,phenoFile=None,type='b',normGenotype=True,readKFile=False):
- self.fbase = fbase
- self.type = type
- self.indivs = self.getIndivs(self.fbase,type)
- self.kFile = kFile
- self.phenos = None
- self.normGenotype = normGenotype
- self.phenoFile = phenoFile
- # Originally I was using the fastLMM style that has indiv IDs embedded.
- # NOW I want to use this module to just read SNPs so I'm allowing
- # the programmer to turn off the kinship reading.
- self.readKFile = readKFile
-
- if self.kFile:
- self.K = self.readKinship(self.kFile)
- elif os.path.isfile("%s.kin" % fbase):
- self.kFile = "%s.kin" %fbase
- if self.readKFile:
- self.K = self.readKinship(self.kFile)
- else:
- self.kFile = None
- self.K = None
-
- self.getPhenos(self.phenoFile)
-
- self.fhandle = None
- self.snpFileHandle = None
-
- def __del__(self):
- if self.fhandle: self.fhandle.close()
- if self.snpFileHandle: self.snpFileHandle.close()
-
- def getSNPIterator(self):
- if not self.type == 'b':
- sys.stderr.write("Have only implemented this for binary plink files (bed)\n")
- return
-
- # get the number of snps
- file = self.fbase + '.bim'
- i = 0
- f = open(file,'r')
- for line in f: i += 1
- f.close()
- self.numSNPs = i
- self.have_read = 0
- self.snpFileHandle = open(file,'r')
-
- self.BytestoRead = self.N / 4 + (self.N % 4 and 1 or 0)
- self._formatStr = 'c'*self.BytestoRead
-
- file = self.fbase + '.bed'
- self.fhandle = open(file,'rb')
-
- magicNumber = self.fhandle.read(2)
- order = self.fhandle.read(1)
- if not order == '\x01':
- sys.stderr.write("This is not in SNP major order - you did not handle this case\n")
- raise StopIteration
-
- return self
-
- def __iter__(self):
- return self.getSNPIterator()
-
- def next(self):
- if self.have_read == self.numSNPs:
- raise StopIteration
- X = self.fhandle.read(self.BytestoRead)
- XX = [bin(ord(x)) for x in struct.unpack(self._formatStr,X)]
- self.have_read += 1
- return self.formatBinaryGenotypes(XX,self.normGenotype),self.snpFileHandle.readline().strip().split()[1]
-
- def formatBinaryGenotypes(self,X,norm=True):
- D = { \
- '00': 0.0, \
- '10': 0.5, \
- '11': 1.0, \
- '01': np.nan \
- }
-
- D_tped = { \
- '00': '1 1', \
- '10': '1 2', \
- '11': '2 2', \
- '01': '0 0' \
- }
-
- #D = D_tped
-
- G = []
- for x in X:
- if not len(x) == 10:
- xx = x[2:]
- x = '0b' + '0'*(8 - len(xx)) + xx
- a,b,c,d = (x[8:],x[6:8],x[4:6],x[2:4])
- L = [D[y] for y in [a,b,c,d]]
- G += L
- # only take the leading values because whatever is left should be null
- G = G[:self.N]
- G = np.array(G)
- if norm:
- G = self.normalizeGenotype(G)
- return G
-
- def normalizeGenotype(self,G):
- # print "Before",G
- # print G.shape
- print "call input.normalizeGenotype"
- raise "This should not be used"
- x = True - np.isnan(G)
- m = G[x].mean()
- s = np.sqrt(G[x].var())
- G[np.isnan(G)] = m
- if s == 0: G = G - m
- else: G = (G - m) / s
- # print "After",G
- return G
-
- def getPhenos(self,phenoFile=None):
- if not phenoFile:
- self.phenoFile = phenoFile = self.fbase+".phenos"
- if not os.path.isfile(phenoFile):
- sys.stderr.write("Could not find phenotype file: %s\n" % (phenoFile))
- return
- f = open(phenoFile,'r')
- keys = []
- P = []
- for line in f:
- v = line.strip().split()
- keys.append((v[0],v[1]))
- P.append([(x == 'NA' or x == '-9') and np.nan or float(x) for x in v[2:]])
- f.close()
- P = np.array(P)
-
- # reorder to match self.indivs
- D = {}
- L = []
- for i in range(len(keys)):
- D[keys[i]] = i
- for i in range(len(self.indivs)):
- if not D.has_key(self.indivs[i]):
- continue
- L.append(D[self.indivs[i]])
- P = P[L,:]
-
- self.phenos = P
- return P
-
- def getIndivs(self,base,type='b'):
- if type == 't':
- famFile = "%s.tfam" % base
- else:
- famFile = "%s.fam" % base
- keys = []
- i = 0
- f = open(famFile,'r')
- for line in f:
- v = line.strip().split()
- famId = v[0]
- indivId = v[1]
- k = (famId.strip(),indivId.strip())
- keys.append(k)
- i += 1
- f.close()
-
- self.N = len(keys)
- sys.stderr.write("Read %d individuals from %s\n" % (self.N, famFile))
-
- return keys
-
- def readKinship(self,kFile):
- # Assume the fastLMM style
- # This will read in the kinship matrix and then reorder it
- # according to self.indivs - additionally throwing out individuals
- # that are not in both sets
- if self.indivs == None or len(self.indivs) == 0:
- sys.stderr.write("Did not read any individuals so can't load kinship\n")
- return
-
- sys.stderr.write("Reading kinship matrix from %s\n" % (kFile) )
-
- f = open(kFile,'r')
- # read indivs
- v = f.readline().strip().split("\t")[1:]
- keys = [tuple(y.split()) for y in v]
- D = {}
- for i in range(len(keys)): D[keys[i]] = i
-
- # read matrix
- K = []
- for line in f:
- K.append([float(x) for x in line.strip().split("\t")[1:]])
- f.close()
- K = np.array(K)
-
- # reorder to match self.indivs
- L = []
- KK = []
- X = []
- for i in range(len(self.indivs)):
- if not D.has_key(self.indivs[i]):
- X.append(self.indivs[i])
- else:
- KK.append(self.indivs[i])
- L.append(D[self.indivs[i]])
- K = K[L,:][:,L]
- self.indivs = KK
- self.indivs_removed = X
- if len(self.indivs_removed):
- sys.stderr.write("Removed %d individuals that did not appear in Kinship\n" % (len(self.indivs_removed)))
- return K
-
- def getCovariates(self,covFile=None):
- if not os.path.isfile(covFile):
- sys.stderr.write("Could not find covariate file: %s\n" % (phenoFile))
- return
- f = open(covFile,'r')
- keys = []
- P = []
- for line in f:
- v = line.strip().split()
- keys.append((v[0],v[1]))
- P.append([x == 'NA' and np.nan or float(x) for x in v[2:]])
- f.close()
- P = np.array(P)
-
- # reorder to match self.indivs
- D = {}
- L = []
- for i in range(len(keys)):
- D[keys[i]] = i
- for i in range(len(self.indivs)):
- if not D.has_key(self.indivs[i]): continue
- L.append(D[self.indivs[i]])
- P = P[L,:]
-
- return P
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
deleted file mode 100644
index 1c157fd8..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
+++ /dev/null
@@ -1,168 +0,0 @@
-# pylmm kinship calculation
-#
-# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com)
-# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-# env PYTHONPATH=$pylmm_lib_path:./lib python $pylmm_lib_path/runlmm.py --pheno test.pheno --geno test9000.geno kinship --test
-
-import sys
-import os
-import numpy as np
-from scipy import linalg
-import multiprocessing as mp # Multiprocessing is part of the Python stdlib
-import Queue
-import time
-
-from optmatrix import matrix_initialize, matrixMultT
-
-# ---- A trick to decide on the environment:
-try:
- from wqflask.my_pylmm.pyLMM import chunks
- from gn2 import uses, progress_set_func
-except ImportError:
- has_gn2=False
- import standalone as handlers
- from standalone import uses, progress_set_func
-
-progress,debug,info,mprint = uses('progress','debug','info','mprint')
-
-def kinship_full(G):
- """
- Calculate the Kinship matrix using a full dot multiplication
- """
- # mprint("kinship_full G",G)
- m = G.shape[0] # snps
- n = G.shape[1] # inds
- info("%d SNPs",m)
- assert m>n, "n should be larger than m (%d snps > %d inds)" % (m,n)
- # m = np.dot(G.T,G)
- m = matrixMultT(G.T)
- m = m/G.shape[0]
- # mprint("kinship_full K",m)
- return m
-
-def compute_W(job,G,n,snps,compute_size):
- """
- Read 1000 SNPs at a time into matrix and return the result
- """
- m = compute_size
- W = np.ones((n,m)) * np.nan # W matrix has dimensions individuals x SNPs (initially all NaNs)
- for j in range(0,compute_size):
- pos = job*m + j # real position
- if pos >= snps:
- W = W[:,range(0,j)]
- break
- snp = G[job*compute_size+j]
- if snp.var() == 0:
- continue
- W[:,j] = snp # set row to list of SNPs
- return W
-
-def compute_matrixMult(job,W,q = None):
- """
- Compute Kinship(W)*j
-
- For every set of SNPs matrixMult is used to multiply matrices T(W)*W
- """
- res = matrixMultT(W)
- if not q: q=compute_matrixMult.q
- q.put([job,res])
- return job
-
-def f_init(q):
- compute_matrixMult.q = q
-
-# Calculate the kinship matrix from G (SNPs as rows!), returns K
-#
-def kinship(G,computeSize=1000,numThreads=None,useBLAS=False):
-
- matrix_initialize(useBLAS)
-
- mprint("G",G)
- n = G.shape[1] # inds
- inds = n
- m = G.shape[0] # snps
- snps = m
- info("%i SNPs" % (m))
- assert snps>=inds, "snps should be larger than inds (%i snps, %i inds)" % (snps,inds)
-
- q = mp.Queue()
- p = mp.Pool(numThreads, f_init, [q])
- cpu_num = mp.cpu_count()
- info("CPU cores: %i" % cpu_num)
- iterations = snps/computeSize+1
-
- results = []
- K = np.zeros((n,n)) # The Kinship matrix has dimension individuals x individuals
-
- completed = 0
- for job in range(iterations):
- info("Processing job %d first %d SNPs" % (job, ((job+1)*computeSize)))
- W = compute_W(job,G,n,snps,computeSize)
- if numThreads == 1:
- # Single-core
- compute_matrixMult(job,W,q)
- j,x = q.get()
- debug("Job "+str(j)+" finished")
- progress("kinship",j,iterations)
- K_j = x
- K = K + K_j
- else:
- # Multi-core
- results.append(p.apply_async(compute_matrixMult, (job,W)))
- # Do we have a result?
- while (len(results)-completed>cpu_num*2):
- time.sleep(0.1)
- try:
- j,x = q.get_nowait()
- debug("Job "+str(j)+" finished")
- K_j = x
- K = K + K_j
- completed += 1
- progress("kinship",completed,iterations)
- except Queue.Empty:
- pass
-
- if numThreads == None or numThreads > 1:
- for job in range(len(results)-completed):
- j,x = q.get(True,15)
- debug("Job "+str(j)+" finished")
- K_j = x
- K = K + K_j
- completed += 1
- progress("kinship",completed,iterations)
-
- K = K / float(snps)
- return K
-
-def kvakve(K):
- """
- Obtain eigendecomposition for K and return Kva,Kve where Kva is cleaned
- of small values < 1e-6 (notably smaller than zero)
- """
- info("Obtaining eigendecomposition for %dx%d matrix" % (K.shape[0],K.shape[1]) )
- Kva,Kve = linalg.eigh(K)
- mprint("Kva",Kva)
- mprint("Kve",Kve)
-
- if sum(Kva < 0):
- info("Cleaning %d eigen values (Kva<0)" % (sum(Kva < 0)))
- Kva[Kva < 1e-6] = 1e-6
- return Kva,Kve
-
-
-
-
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
deleted file mode 100644
index 2a0c7fdc..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
+++ /dev/null
@@ -1,995 +0,0 @@
-# pylmm is a python-based linear mixed-model solver with applications to GWAS
-
-# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-from __future__ import absolute_import, print_function, division
-
-import sys
-import time
-import uuid
-
-import numpy as np
-from scipy import linalg
-from scipy import optimize
-from scipy import stats
-# import pdb
-
-# import gzip
-# import zlib
-import datetime
-# import cPickle as pickle
-from pprint import pformat as pf
-
-# Add local dir to PYTHONPATH
-import os
-cwd = os.path.dirname(__file__)
-if sys.path[0] != cwd:
- sys.path.insert(1,cwd)
-
-# pylmm imports
-from kinship import kinship, kinship_full, kvakve
-import genotype
-import phenotype
-import gwas
-from benchmark import Bench
-
-# The following imports are for exchanging data with the webserver
-import simplejson as json
-from redis import Redis
-Redis = Redis()
-import temp_data
-
-has_gn2=None
-
-# sys.stderr.write("INFO: pylmm system path is "+":".join(sys.path)+"\n")
-sys.stderr.write("INFO: pylmm file is "+__file__+"\n")
-
-# ---- A trick to decide on the environment:
-try:
- sys.stderr.write("INFO: lmm try loading module\n")
- import utility.formatting # this is never used, just to check the environment
- sys.stderr.write("INFO: This is a genenetwork2 environment\n")
- from gn2 import uses, progress_set_func
- has_gn2=True
-except ImportError:
- # Failed to load gn2
- has_gn2=False
- import standalone as handlers
- from standalone import uses, progress_set_func
- sys.stderr.write("WARNING: LMM standalone version missing the Genenetwork2 environment\n")
-
-progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal')
-
-#np.seterr('raise')
-
-#def run_human(pheno_vector,
-# covariate_matrix,
-# plink_input_file,
-# kinship_matrix,
-# refit=False,
-# loading_progress=None):
-
-def run_human(pheno_vector,
- covariate_matrix,
- plink_input_file,
- kinship_matrix,
- refit=False):
-
- v = np.isnan(pheno_vector)
- keep = True - v
- keep = keep.reshape((len(keep),))
-
- identifier = str(uuid.uuid4())
-
- #print("pheno_vector: ", pf(pheno_vector))
- #print("kinship_matrix: ", pf(kinship_matrix))
- #print("kinship_matrix.shape: ", pf(kinship_matrix.shape))
-
- #lmm_vars = pickle.dumps(dict(
- # pheno_vector = pheno_vector,
- # covariate_matrix = covariate_matrix,
- # kinship_matrix = kinship_matrix
- #))
- #Redis.hset(identifier, "lmm_vars", lmm_vars)
- #Redis.expire(identifier, 60*60)
-
- if v.sum():
- pheno_vector = pheno_vector[keep]
- print("pheno_vector shape is now: ", pf(pheno_vector.shape))
- covariate_matrix = covariate_matrix[keep,:]
- print("kinship_matrix shape is: ", pf(kinship_matrix.shape))
- print("keep is: ", pf(keep.shape))
- kinship_matrix = kinship_matrix[keep,:][:,keep]
-
- print("kinship_matrix:", pf(kinship_matrix))
-
- n = kinship_matrix.shape[0]
- print("n is:", n)
- lmm_ob = LMM(pheno_vector,
- kinship_matrix,
- covariate_matrix)
- lmm_ob.fit()
-
-
- # Buffers for pvalues and t-stats
- p_values = []
- t_stats = []
-
- #print("input_file: ", plink_input_file)
-
- with Bench("Opening and loading pickle file"):
- with gzip.open(plink_input_file, "rb") as input_file:
- data = pickle.load(input_file)
-
- plink_input = data['plink_input']
-
- #plink_input.getSNPIterator()
- with Bench("Calculating numSNPs"):
- total_snps = data['numSNPs']
-
- with Bench("snp iterator loop"):
- count = 0
-
- with Bench("Create list of inputs"):
- inputs = list(plink_input)
-
- with Bench("Divide into chunks"):
- results = chunks.divide_into_chunks(inputs, 64)
-
- result_store = []
-
- key = "plink_inputs"
-
- # Todo: Delete below line when done testing
- Redis.delete(key)
-
- timestamp = datetime.datetime.utcnow().isoformat()
-
- # Pickle chunks of input SNPs (from Plink interator) and compress them
- #print("Starting adding loop")
- for part, result in enumerate(results):
- #data = pickle.dumps(result, pickle.HIGHEST_PROTOCOL)
- holder = pickle.dumps(dict(
- identifier = identifier,
- part = part,
- timestamp = timestamp,
- result = result
- ), pickle.HIGHEST_PROTOCOL)
-
- #print("Adding:", part)
- Redis.rpush(key, zlib.compress(holder))
- #print("End adding loop")
- #print("***** Added to {} queue *****".format(key))
- for snp, this_id in plink_input:
- #with Bench("part before association"):
- #if count > 1000:
- # break
- count += 1
- progress("human",count,total_snps)
-
- #with Bench("actual association"):
- ps, ts = human_association(snp,
- n,
- keep,
- lmm_ob,
- pheno_vector,
- covariate_matrix,
- kinship_matrix,
- refit)
-
- #with Bench("after association"):
- p_values.append(ps)
- t_stats.append(ts)
-
- return p_values, t_stats
-
-
-#class HumanAssociation(object):
-# def __init__(self):
-#
-
-def human_association(snp,
- n,
- keep,
- lmm_ob,
- pheno_vector,
- covariate_matrix,
- kinship_matrix,
- refit):
-
- x = snp[keep].reshape((n,1))
- #x[[1,50,100,200,3000],:] = np.nan
- v = np.isnan(x).reshape((-1,))
-
- # Check SNPs for missing values
- if v.sum():
- keeps = True - v
- xs = x[keeps,:]
- # If no variation at this snp or all genotypes missing
- if keeps.sum() <= 1 or xs.var() <= 1e-6:
- return np.nan, np.nan
- #p_values.append(np.nan)
- #t_stats.append(np.nan)
- #continue
-
- # Its ok to center the genotype - I used options.normalizeGenotype to
- # force the removal of missing genotypes as opposed to replacing them with MAF.
-
- #if not options.normalizeGenotype:
- # xs = (xs - xs.mean()) / np.sqrt(xs.var())
-
- filtered_pheno = pheno_vector[keeps]
- filtered_covariate_matrix = covariate_matrix[keeps,:]
-
- print("kinship_matrix shape is: ", pf(kinship_matrix.shape))
- print("keeps is: ", pf(keeps.shape))
- filtered_kinship_matrix = kinship_matrix[keeps,:][:,keeps]
- filtered_lmm_ob = lmm.LMM(filtered_pheno,filtered_kinship_matrix,X0=filtered_covariate_matrix)
- if refit:
- filtered_lmm_ob.fit(X=xs)
- else:
- #try:
- filtered_lmm_ob.fit()
- #except: pdb.set_trace()
- ts,ps,beta,betaVar = Ls.association(xs,returnBeta=True)
- else:
- if x.var() == 0:
- return np.nan, np.nan
- #p_values.append(np.nan)
- #t_stats.append(np.nan)
- #continue
- if refit:
- lmm_ob.fit(X=x)
- ts, ps, beta, betaVar = lmm_ob.association(x)
- return ps, ts
-
-
-#def run(pheno_vector,
-# genotype_matrix,
-# restricted_max_likelihood=True,
-# refit=False,
-# temp_data=None):
-
-def run_other_old(pheno_vector,
- genotype_matrix,
- restricted_max_likelihood=True,
- refit=False):
-
- """Takes the phenotype vector and genotype matrix and returns a set of p-values and t-statistics
-
- restricted_max_likelihood -- whether to use restricted max likelihood; True or False
- refit -- whether to refit the variance component for each marker
-
- """
-
- print("Running the original LMM engine in run_other (old)")
- print("REML=",restricted_max_likelihood," REFIT=",refit)
- with Bench("Calculate Kinship"):
- kinship_matrix,genotype_matrix = calculate_kinship_new(genotype_matrix)
-
- print("kinship_matrix: ", pf(kinship_matrix))
- print("kinship_matrix.shape: ", pf(kinship_matrix.shape))
-
- # with Bench("Create LMM object"):
- # lmm_ob = LMM(pheno_vector, kinship_matrix)
-
- # with Bench("LMM_ob fitting"):
- # lmm_ob.fit()
-
- print("run_other_old genotype_matrix: ", genotype_matrix.shape)
- print(genotype_matrix)
-
- with Bench("Doing GWAS"):
- t_stats, p_values = GWAS(pheno_vector,
- genotype_matrix.T,
- kinship_matrix,
- restricted_max_likelihood=True,
- refit=False)
- Bench().report()
- return p_values, t_stats
-
-def run_other_new(n,m,pheno_vector,
- geno,
- restricted_max_likelihood=True,
- refit=False):
-
- """Takes the phenotype vector and genotype matrix and returns a set of p-values and t-statistics
-
- restricted_max_likelihood -- whether to use restricted max likelihood; True or False
- refit -- whether to refit the variance component for each marker
-
- """
-
- print("Running the new LMM2 engine in run_other_new")
- print("REML=",restricted_max_likelihood," REFIT=",refit)
-
- # Adjust phenotypes
- n,Y,keep = phenotype.remove_missing_new(n,pheno_vector)
-
- # if options.maf_normalization:
- # G = np.apply_along_axis( genotype.replace_missing_with_MAF, axis=0, arr=g )
- # print "MAF replacements: \n",G
- # if not options.skip_genotype_normalization:
- # G = np.apply_along_axis( genotype.normalize, axis=1, arr=G)
-
- geno = geno[:,keep]
- with Bench("Calculate Kinship"):
- K,G = calculate_kinship_new(geno)
-
- print("kinship_matrix: ", pf(K))
- print("kinship_matrix.shape: ", pf(K.shape))
-
- # with Bench("Create LMM object"):
- # lmm_ob = lmm2.LMM2(Y,K)
- # with Bench("LMM_ob fitting"):
- # lmm_ob.fit()
-
- print("run_other_new genotype_matrix: ", G.shape)
- print(G)
-
- with Bench("Doing GWAS"):
- t_stats, p_values = gwas.gwas(Y,
- G,
- K,
- restricted_max_likelihood=True,
- refit=False,verbose=True)
- Bench().report()
- return p_values, t_stats
-
-# def matrixMult(A,B):
-# return np.dot(A,B)
-
-def matrixMult(A,B):
-
- # If there is no fblas then we will revert to np.dot()
-
- try:
- linalg.fblas
- except AttributeError:
- return np.dot(A,B)
-
- #print("A is:", pf(A.shape))
- #print("B is:", pf(B.shape))
-
- # If the matrices are in Fortran order then the computations will be faster
- # when using dgemm. Otherwise, the function will copy the matrix and that takes time.
- if not A.flags['F_CONTIGUOUS']:
- AA = A.T
- transA = True
- else:
- AA = A
- transA = False
-
- if not B.flags['F_CONTIGUOUS']:
- BB = B.T
- transB = True
- else:
- BB = B
- transB = False
-
- return linalg.fblas.dgemm(alpha=1.,a=AA,b=BB,trans_a=transA,trans_b=transB)
-
-def calculate_kinship_new(genotype_matrix):
- """
- Call the new kinship calculation where genotype_matrix contains
- inds (columns) by snps (rows).
- """
- assert type(genotype_matrix) is np.ndarray
- info("call genotype.normalize")
- G = np.apply_along_axis( genotype.normalize, axis=1, arr=genotype_matrix)
- mprint("G",genotype_matrix)
- info("call calculate_kinship_new")
- return kinship(G),G # G gets transposed, we'll turn this into an iterator (FIXME)
-
-def calculate_kinship_iter(geno):
- """
- Call the new kinship calculation where genotype_matrix contains
- inds (columns) by snps (rows).
- """
- assert type(genotype_matrix) is iter
- info("call genotype.normalize")
- G = np.apply_along_axis( genotype.normalize, axis=0, arr=genotype_matrix)
- info("call calculate_kinship_new")
- return kinship(G)
-
-def calculate_kinship_old(genotype_matrix):
- """
- genotype_matrix is an n x m matrix encoding SNP minor alleles.
-
- This function takes a matrix oF SNPs, imputes missing values with the maf,
- normalizes the resulting vectors and returns the RRM matrix.
-
- """
- info("call calculate_kinship_old")
- fatal("THE FUNCTION calculate_kinship_old IS OBSOLETE, use calculate_kinship_new instead - see Genotype Normalization Problem on Pjotr's blog")
- n = genotype_matrix.shape[0]
- m = genotype_matrix.shape[1]
- info("genotype 2D matrix n (inds) is: %d" % (n))
- info("genotype 2D matrix m (snps) is: %d" % (m))
- assert m>n, "n should be larger than m (snps>inds)"
- keep = []
- mprint("G (before old normalize)",genotype_matrix)
- for counter in range(m):
- #print("type of genotype_matrix[:,counter]:", pf(genotype_matrix[:,counter]))
- #Checks if any values in column are not numbers
- not_number = np.isnan(genotype_matrix[:,counter])
-
- #Gets vector of values for column (no values in vector if not all values in col are numbers)
- marker_values = genotype_matrix[True - not_number, counter]
- #print("marker_values is:", pf(marker_values))
-
- #Gets mean of values in vector
- values_mean = marker_values.mean()
-
- genotype_matrix[not_number,counter] = values_mean
- vr = genotype_matrix[:,counter].var()
- if vr == 0:
- continue
- keep.append(counter)
- genotype_matrix[:,counter] = (genotype_matrix[:,counter] - values_mean) / np.sqrt(vr)
- progress('kinship_old normalize genotype',counter,m)
-
- genotype_matrix = genotype_matrix[:,keep]
- mprint("G (after old normalize)",genotype_matrix.T)
- kinship_matrix = np.dot(genotype_matrix, genotype_matrix.T) * 1.0/float(m)
- return kinship_matrix,genotype_matrix
- # return kinship_full(genotype_matrix.T),genotype_matrix
-
-def GWAS(pheno_vector,
- genotype_matrix,
- kinship_matrix,
- kinship_eigen_vals=None,
- kinship_eigen_vectors=None,
- covariate_matrix=None,
- restricted_max_likelihood=True,
- refit=False,
- temp_data=None):
- """
- Performs a basic GWAS scan using the LMM. This function
- uses the LMM module to assess association at each SNP and
- does some simple cleanup, such as removing missing individuals
- per SNP and re-computing the eigen-decomp
-
- pheno_vector - n x 1 phenotype vector
- genotype_matrix - n x m SNP matrix
- kinship_matrix - n x n kinship matrix
- kinship_eigen_vals, kinship_eigen_vectors = linalg.eigh(K) - or the eigen vectors and values for K
- covariate_matrix - n x q covariate matrix
- restricted_max_likelihood - use restricted maximum likelihood
- refit - refit the variance component for each SNP
-
- """
- if kinship_eigen_vals is None:
- kinship_eigen_vals = []
- if kinship_eigen_vectors is None:
- kinship_eigen_vectors = []
-
- n = genotype_matrix.shape[0]
- m = genotype_matrix.shape[1]
-
- if covariate_matrix == None:
- covariate_matrix = np.ones((n,1))
-
- # Remove missing values in pheno_vector and adjust associated parameters
- v = np.isnan(pheno_vector)
- if v.sum():
- keep = True - v
- print(pheno_vector.shape,pheno_vector)
- print(keep.shape,keep)
- pheno_vector = pheno_vector[keep]
- #genotype_matrix = genotype_matrix[keep,:]
- #covariate_matrix = covariate_matrix[keep,:]
- #kinship_matrix = kinship_matrix[keep,:][:,keep]
- kinship_eigen_vals = []
- kinship_eigen_vectors = []
-
- lmm_ob = LMM(pheno_vector,
- kinship_matrix,
- kinship_eigen_vals,
- kinship_eigen_vectors,
- covariate_matrix)
- if not refit:
- lmm_ob.fit()
-
- p_values = []
- t_statistics = []
-
- n = genotype_matrix.shape[0]
- m = genotype_matrix.shape[1]
-
- for counter in range(m):
- x = genotype_matrix[:,counter].reshape((n, 1))
- v = np.isnan(x).reshape((-1,))
- if v.sum():
- keep = True - v
- xs = x[keep,:]
- if xs.var() == 0:
- p_values.append(0)
- t_statistics.append(np.nan)
- continue
-
- print(genotype_matrix.shape,pheno_vector.shape,keep.shape)
-
- pheno_vector = pheno_vector[keep]
- covariate_matrix = covariate_matrix[keep,:]
- kinship_matrix = kinship_matrix[keep,:][:,keep]
- lmm_ob_2 = LMM(pheno_vector,
- kinship_matrix,
- X0=covariate_matrix)
- if refit:
- lmm_ob_2.fit(X=xs)
- else:
- lmm_ob_2.fit()
- ts, ps, beta, betaVar = lmm_ob_2.association(xs, REML=restricted_max_likelihood)
- else:
- if x.var() == 0:
- p_values.append(0)
- t_statistics.append(np.nan)
- continue
-
- if refit:
- lmm_ob.fit(X=x)
- ts, ps, beta, betaVar = lmm_ob.association(x, REML=restricted_max_likelihood)
-
- progress("gwas_old",counter,m)
-
- p_values.append(ps)
- t_statistics.append(ts)
-
- return t_statistics, p_values
-
-
-class LMM:
-
- """
- This is a simple version of EMMA/fastLMM.
- The main purpose of this module is to take a phenotype vector (Y), a set of covariates (X) and a kinship matrix (K)
- and to optimize this model by finding the maximum-likelihood estimates for the model parameters.
- There are three model parameters: heritability (h), covariate coefficients (beta) and the total
- phenotypic variance (sigma).
- Heritability as defined here is the proportion of the total variance (sigma) that is attributed to
- the kinship matrix.
-
- For simplicity, we assume that everything being input is a numpy array.
- If this is not the case, the module may throw an error as conversion from list to numpy array
- is not done consistently.
-
- """
- def __init__(self,Y,K,Kva=[],Kve=[],X0=None,verbose=True):
-
- """
- The constructor takes a phenotype vector or array of size n.
- It takes a kinship matrix of size n x n. Kva and Kve can be computed as Kva,Kve = linalg.eigh(K) and cached.
- If they are not provided, the constructor will calculate them.
- X0 is an optional covariate matrix of size n x q, where there are q covariates.
- When this parameter is not provided, the constructor will set X0 to an n x 1 matrix of all ones to represent a mean effect.
- """
-
- if X0 is None: X0 = np.ones(len(Y)).reshape(len(Y),1)
- self.verbose = verbose
-
- #x = Y != -9
- x = True - np.isnan(Y)
- #pdb.set_trace()
- if not x.sum() == len(Y):
- print("Removing %d missing values from Y\n" % ((True - x).sum()))
- if self.verbose: sys.stderr.write("Removing %d missing values from Y\n" % ((True - x).sum()))
- Y = Y[x]
- print("x: ", len(x))
- print("K: ", K.shape)
- #K = K[x,:][:,x]
- X0 = X0[x,:]
- Kva = []
- Kve = []
- self.nonmissing = x
-
- print("this K is:", K.shape, pf(K))
-
- if len(Kva) == 0 or len(Kve) == 0:
- # if self.verbose: sys.stderr.write("Obtaining eigendecomposition for %dx%d matrix\n" % (K.shape[0],K.shape[1]) )
- begin = time.time()
- # Kva,Kve = linalg.eigh(K)
- Kva,Kve = kvakve(K)
- end = time.time()
- if self.verbose: sys.stderr.write("Total time: %0.3f\n" % (end - begin))
- print("sum(Kva),sum(Kve)=",sum(Kva),sum(Kve))
-
- self.K = K
- self.Kva = Kva
- self.Kve = Kve
- print("self.Kva is: ", self.Kva.shape, pf(self.Kva))
- print("self.Kve is: ", self.Kve.shape, pf(self.Kve))
- self.Y = Y
- self.X0 = X0
- self.N = self.K.shape[0]
-
- # ----> Below moved to kinship.kvakve(K)
- # if sum(self.Kva < 1e-6):
- # if self.verbose: sys.stderr.write("Cleaning %d eigen values\n" % (sum(self.Kva < 0)))
- # self.Kva[self.Kva < 1e-6] = 1e-6
-
- self.transform()
-
- def transform(self):
-
- """
- Computes a transformation on the phenotype vector and the covariate matrix.
- The transformation is obtained by left multiplying each parameter by the transpose of the
- eigenvector matrix of K (the kinship).
- """
-
- self.Yt = matrixMult(self.Kve.T, self.Y)
- self.X0t = matrixMult(self.Kve.T, self.X0)
- self.X0t_stack = np.hstack([self.X0t, np.ones((self.N,1))])
- self.q = self.X0t.shape[1]
-
- def getMLSoln(self,h,X):
-
- """
- Obtains the maximum-likelihood estimates for the covariate coefficients (beta),
- the total variance of the trait (sigma) and also passes intermediates that can
- be utilized in other functions. The input parameter h is a value between 0 and 1 and represents
- the heritability or the proportion of the total variance attributed to genetics. The X is the
- covariate matrix.
- """
-
- S = 1.0/(h*self.Kva + (1.0 - h))
- Xt = X.T*S
- XX = matrixMult(Xt,X)
- XX_i = linalg.inv(XX)
- beta = matrixMult(matrixMult(XX_i,Xt),self.Yt)
- Yt = self.Yt - matrixMult(X,beta)
- Q = np.dot(Yt.T*S,Yt)
- sigma = Q * 1.0 / (float(self.N) - float(X.shape[1]))
- return beta,sigma,Q,XX_i,XX
-
- def LL_brent(self,h,X=None,REML=False):
- #brent will not be bounded by the specified bracket.
- # I return a large number if we encounter h < 0 to avoid errors in LL computation during the search.
- if h < 0: return 1e6
- return -self.LL(h,X,stack=False,REML=REML)[0]
-
- def LL(self,h,X=None,stack=True,REML=False):
-
- """
- Computes the log-likelihood for a given heritability (h). If X==None, then the
- default X0t will be used. If X is set and stack=True, then X0t will be matrix concatenated with
- the input X. If stack is false, then X is used in place of X0t in the LL calculation.
- REML is computed by adding additional terms to the standard LL and can be computed by setting REML=True.
- """
-
- if X is None:
- X = self.X0t
- elif stack:
- self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0]
- X = self.X0t_stack
-
- n = float(self.N)
- q = float(X.shape[1])
- beta,sigma,Q,XX_i,XX = self.getMLSoln(h,X)
- LL = n*np.log(2*np.pi) + np.log(h*self.Kva + (1.0-h)).sum() + n + n*np.log(1.0/n * Q)
- LL = -0.5 * LL
-
- if REML:
- LL_REML_part = q*np.log(2.0*np.pi*sigma) + np.log(linalg.det(matrixMult(X.T,X))) - np.log(linalg.det(XX))
- LL = LL + 0.5*LL_REML_part
-
- return LL,beta,sigma,XX_i
-
- def getMax(self,H, X=None,REML=False):
-
- """
- Helper functions for .fit(...).
- This function takes a set of LLs computed over a grid and finds possible regions
- containing a maximum. Within these regions, a Brent search is performed to find the
- optimum.
-
- """
- n = len(self.LLs)
- HOpt = []
- for i in range(1,n-2):
- if self.LLs[i-1] < self.LLs[i] and self.LLs[i] > self.LLs[i+1]:
- HOpt.append(optimize.brent(self.LL_brent,args=(X,REML),brack=(H[i-1],H[i+1])))
- if np.isnan(HOpt[-1][0]):
- HOpt[-1][0] = [self.LLs[i-1]]
-
- if len(HOpt) > 1:
- if self.verbose:
- sys.stderr.write("NOTE: Found multiple optima. Returning first...\n")
- return HOpt[0]
- elif len(HOpt) == 1:
- return HOpt[0]
- elif self.LLs[0] > self.LLs[n-1]:
- return H[0]
- else:
- return H[n-1]
-
- def fit(self,X=None,ngrids=100,REML=True):
-
- """
- Finds the maximum-likelihood solution for the heritability (h) given the current parameters.
- X can be passed and will transformed and concatenated to X0t. Otherwise, X0t is used as
- the covariate matrix.
-
- This function calculates the LLs over a grid and then uses .getMax(...) to find the optimum.
- Given this optimum, the function computes the LL and associated ML solutions.
- """
-
- if X == None:
- X = self.X0t
- else:
- #X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)])
- self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0]
- X = self.X0t_stack
-
- H = np.array(range(ngrids)) / float(ngrids)
- L = np.array([self.LL(h,X,stack=False,REML=REML)[0] for h in H])
- self.LLs = L
-
- hmax = self.getMax(H,X,REML)
- L,beta,sigma,betaSTDERR = self.LL(hmax,X,stack=False,REML=REML)
-
- self.H = H
- self.optH = hmax
- self.optLL = L
- self.optBeta = beta
- self.optSigma = sigma
-
- return hmax,beta,sigma,L
-
- def association(self,X, h = None, stack=True,REML=True, returnBeta=True):
-
- """
- Calculates association statitics for the SNPs encoded in the vector X of size n.
- If h == None, the optimal h stored in optH is used.
-
- """
- if stack:
- #X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)])
- self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0]
- X = self.X0t_stack
-
- if h == None:
- h = self.optH
-
- L,beta,sigma,betaVAR = self.LL(h,X,stack=False,REML=REML)
- q = len(beta)
- ts,ps = self.tstat(beta[q-1],betaVAR[q-1,q-1],sigma,q)
-
- if returnBeta:
- return ts,ps,beta[q-1].sum(),betaVAR[q-1,q-1].sum()*sigma
- return ts,ps
-
- def tstat(self,beta,var,sigma,q):
-
- """
- Calculates a t-statistic and associated p-value given the estimate of beta and its standard error.
- This is actually an F-test, but when only one hypothesis is being performed, it reduces to a t-test.
- """
-
- ts = beta / np.sqrt(var * sigma)
- ps = 2.0*(1.0 - stats.t.cdf(np.abs(ts), self.N-q))
- if not len(ts) == 1 or not len(ps) == 1:
- print("ts=",ts)
- print("ps=",ps)
- raise Exception("Something bad happened :(")
- return ts.sum(),ps.sum()
-
- def plotFit(self,color='b-',title=''):
-
- """
- Simple function to visualize the likelihood space. It takes the LLs
- calcualted over a grid and normalizes them by subtracting off the mean and exponentiating.
- The resulting "probabilities" are normalized to one and plotted against heritability.
- This can be seen as an approximation to the posterior distribuiton of heritability.
-
- For diagnostic purposes this lets you see if there is one distinct maximum or multiple
- and what the variance of the parameter looks like.
- """
- import matplotlib.pyplot as pl
-
- mx = self.LLs.max()
- p = np.exp(self.LLs - mx)
- p = p/p.sum()
-
- pl.plot(self.H,p,color)
- pl.xlabel("Heritability")
- pl.ylabel("Probability of data")
- pl.title(title)
-
-def run_gwas(species,n,m,k,y,geno,cov=None,reml=True,refit=False,inputfn=None,new_code=True):
- """
- Invoke pylmm using genotype as a matrix or as a (SNP) iterator.
- """
- info("run_gwas")
- print('pheno', y)
-
- if species == "human" :
- print('kinship', k )
- ps, ts = run_human(pheno_vector = y,
- covariate_matrix = cov,
- plink_input_file = inputfn,
- kinship_matrix = k,
- refit = refit)
- else:
- print('geno', geno.shape, geno)
-
- if new_code:
- ps, ts = run_other_new(n,m,pheno_vector = y,
- geno = geno,
- restricted_max_likelihood = reml,
- refit = refit)
- else:
- ps, ts = run_other_old(pheno_vector = y,
- genotype_matrix = geno,
- restricted_max_likelihood = reml,
- refit = refit)
- return ps,ts
-
-def gwas_with_redis(key,species,new_code=True):
- """
- Invoke pylmm using Redis as a container. new_code runs the new
- version. All the Redis code goes here!
- """
- json_params = Redis.get(key)
-
- params = json.loads(json_params)
-
- tempdata = temp_data.TempData(params['temp_uuid'])
- def update_tempdata(loc,i,total):
- """
- This is the single method that updates Redis for percentage complete!
- """
- tempdata.store("percent_complete",round(i*100.0/total))
- debug("Updating REDIS percent_complete=%d" % (round(i*100.0/total)))
- progress_set_func(update_tempdata)
-
- def narray(t):
- info("Type is "+t)
- v = params.get(t)
- if v is not None:
- # Note input values can be array of string or float
- v1 = [x if x != 'NA' else 'nan' for x in v]
- v = np.array(v1).astype(np.float)
- return v
-
- def marray(t):
- info("Type is "+t)
- v = params.get(t)
- if v is not None:
- m = []
- for r in v:
- # Note input values can be array of string or float
- r1 = [x if x != 'NA' else 'nan' for x in r]
- m.append(np.array(r1).astype(np.float))
- return np.array(m)
- return np.array(v)
-
- def marrayT(t):
- m = marray(t)
- if m is not None:
- return m.T
- return m
-
- # We are transposing before we enter run_gwas - this should happen on the webserver
- # side (or when reading data from file)
- k = marray('kinship_matrix')
- g = marrayT('genotype_matrix')
- mprint("geno",g)
- y = narray('pheno_vector')
- n = len(y)
- m = params.get('num_genotypes')
- if m is None:
- m = g.shape[0]
- info("m=%d,n=%d" % (m,n))
- ps,ts = run_gwas(species,n,m,k,y,g,narray('covariate_matrix'),params['restricted_max_likelihood'],params['refit'],params.get('input_file_name'),new_code)
-
- results_key = "pylmm:results:" + params['temp_uuid']
-
- # fatal(results_key)
- json_results = json.dumps(dict(p_values = ps,
- t_stats = ts))
-
- #Pushing json_results into a list where it is the only item because blpop needs a list
- Redis.rpush(results_key, json_results)
- Redis.expire(results_key, 60*60)
- return ps, ts
-
-def gn2_load_redis(key,species,kinship,pheno,geno,new_code=True):
- """
- This function emulates current GN2 behaviour by pre-loading Redis (note the input
- genotype is transposed to emulate GN2 (FIXME!)
- """
- info("Loading Redis from parsed data")
- if kinship == None:
- k = None
- else:
- k = kinship.tolist()
- params = dict(pheno_vector = pheno.tolist(),
- genotype_matrix = geno.T.tolist(),
- num_genotypes = geno.shape[0],
- kinship_matrix = k,
- covariate_matrix = None,
- input_file_name = None,
- restricted_max_likelihood = True,
- refit = False,
- temp_uuid = "testrun_temp_uuid",
-
- # meta data
- timestamp = datetime.datetime.now().isoformat())
-
- json_params = json.dumps(params)
- Redis.set(key, json_params)
- Redis.expire(key, 60*60)
-
- return gwas_with_redis(key,species,new_code)
-
-def gn2_load_redis_iter(key,species,kinship,pheno,geno_iterator):
- """
- This function emulates GN2 behaviour by pre-loading Redis with
- a SNP iterator, for this it sets a key for every genotype (SNP)
- """
- print("Loading Redis using a SNP iterator")
- for i,genotypes in enumerate(geno_iterator):
- gkey = key+'_geno_'+str(i)
- Redis.set(gkey, genotypes)
- Redis.expire(gkey, 60*60)
-
- if kinship == None:
- k = None
- else:
- k = kinship.tolist()
- params = dict(pheno_vector = pheno.tolist(),
- genotype_matrix = "iterator",
- num_genotypes = i,
- kinship_matrix = k,
- covariate_matrix = None,
- input_file_name = None,
- restricted_max_likelihood = True,
- refit = False,
- temp_uuid = "testrun_temp_uuid",
-
- # meta data
- timestamp = datetime.datetime.now().isoformat(),
- )
-
- json_params = json.dumps(params)
- Redis.set(key, json_params)
- Redis.expire(key, 60*60)
- return gwas_with_redis(key,species)
-
-# This is the main function used by Genenetwork2 (with environment)
-#
-# Note that this calling route will become OBSOLETE (we should use runlmm.py
-# instead)
-def gn2_main():
- import argparse
- parser = argparse.ArgumentParser(description='Run pyLMM')
- parser.add_argument('-k', '--key')
- parser.add_argument('-s', '--species')
-
- opts = parser.parse_args()
-
- key = opts.key
- species = opts.species
-
- gwas_with_redis(key,species)
-
-
-if __name__ == '__main__':
- print("WARNING: Calling pylmm from lmm.py will become OBSOLETE, use runlmm.py instead!")
- gn2_main()
-
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py
deleted file mode 100644
index d871d8d2..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py
+++ /dev/null
@@ -1,433 +0,0 @@
-# pylmm is a python-based linear mixed-model solver with applications to GWAS
-
-# Copyright (C) 2013,2014 Nicholas A. Furlotte (nick.furlotte@gmail.com)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-import sys
-import time
-import numpy as np
-from scipy.linalg import eigh, inv, det
-import scipy.stats as stats # t-tests
-from scipy import optimize
-from optmatrix import matrixMult
-import kinship
-
-sys.stderr.write("INFO: pylmm (lmm2) system path is "+":".join(sys.path)+"\n")
-sys.stderr.write("INFO: pylmm (lmm2) file is "+__file__+"\n")
-
-# ---- A trick to decide on the environment:
-try:
- sys.stderr.write("INFO: lmm2 try loading module\n")
- import utility.formatting # this is never used, just to check the environment
- sys.stderr.write("INFO: This is a genenetwork2 environment (lmm2)\n")
- from gn2 import uses, progress_set_func
-except ImportError:
- # Failed to load gn2
- has_gn2=False
- import standalone as handlers
- from standalone import uses, progress_set_func
- sys.stderr.write("WARNING: LMM2 standalone version missing the Genenetwork2 environment\n")
-
-progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal')
-
-
-def calculateKinship(W,center=False):
- """
- W is an n x m matrix encoding SNP minor alleles.
-
- This function takes a matrix oF SNPs, imputes missing values with the maf,
- normalizes the resulting vectors and returns the RRM matrix.
- """
- n = W.shape[0]
- m = W.shape[1]
- keep = []
- for i in range(m):
- mn = W[True - np.isnan(W[:,i]),i].mean()
- W[np.isnan(W[:,i]),i] = mn
- vr = W[:,i].var()
- if vr == 0: continue
-
- keep.append(i)
- W[:,i] = (W[:,i] - mn) / np.sqrt(vr)
-
- W = W[:,keep]
- K = matrixMult(W,W.T) * 1.0/float(m)
- if center:
- P = np.diag(np.repeat(1,n)) - 1/float(n) * np.ones((n,n))
- S = np.trace(matrixMult(matrixMult(P,K),P))
- K_n = (n - 1)*K / S
- return K_n
- return K
-
-def GWAS(Y, X, K, Kva=[], Kve=[], X0=None, REML=True, refit=False):
- """
-
- Performs a basic GWAS scan using the LMM. This function
- uses the LMM module to assess association at each SNP and
- does some simple cleanup, such as removing missing individuals
- per SNP and re-computing the eigen-decomp
-
- Y - n x 1 phenotype vector
- X - n x m SNP matrix (genotype matrix)
- K - n x n kinship matrix
- Kva,Kve = linalg.eigh(K) - or the eigen vectors and values for K
- X0 - n x q covariate matrix
- REML - use restricted maximum likelihood
- refit - refit the variance component for each SNP
-
- """
- n = X.shape[0]
- m = X.shape[1]
- prins("Initialize GWAS")
- print("genotype matrix n is:", n)
- print("genotype matrix m is:", m)
-
- if X0 is None:
- X0 = np.ones((n,1))
-
- # Remove missing values in Y and adjust associated parameters
- v = np.isnan(Y)
- if v.sum():
- keep = True - v
- keep = keep.reshape((-1,))
- Y = Y[keep]
- X = X[keep,:]
- X0 = X0[keep,:]
- K = K[keep,:][:,keep]
- Kva = []
- Kve = []
-
- if len(Y) == 0:
- return np.ones(m)*np.nan,np.ones(m)*np.nan
-
- L = LMM(Y,K,Kva,Kve,X0)
- if not refit: L.fit()
-
- PS = []
- TS = []
-
- n = X.shape[0]
- m = X.shape[1]
-
- for i in range(m):
- x = X[:,i].reshape((n,1))
- v = np.isnan(x).reshape((-1,))
- if v.sum():
- keep = True - v
- xs = x[keep,:]
- if xs.var() == 0:
- PS.append(np.nan)
- TS.append(np.nan)
- continue
-
- Ys = Y[keep]
- X0s = X0[keep,:]
- Ks = K[keep,:][:,keep]
- Ls = LMM(Ys,Ks,X0=X0s)
- if refit:
- Ls.fit(X=xs)
- else:
- Ls.fit()
- ts,ps = Ls.association(xs,REML=REML)
- else:
- if x.var() == 0:
- PS.append(np.nan)
- TS.append(np.nan)
- continue
-
- if refit:
- L.fit(X=x)
- ts,ps = L.association(x,REML=REML)
-
- PS.append(ps)
- TS.append(ts)
-
- return TS,PS
-
-class LMM2:
-
- """This is a simple version of EMMA/fastLMM.
-
- The main purpose of this module is to take a phenotype vector (Y),
- a set of covariates (X) and a kinship matrix (K) and to optimize
- this model by finding the maximum-likelihood estimates for the
- model parameters. There are three model parameters: heritability
- (h), covariate coefficients (beta) and the total phenotypic
- variance (sigma). Heritability as defined here is the proportion
- of the total variance (sigma) that is attributed to the kinship
- matrix.
-
- For simplicity, we assume that everything being input is a numpy
- array. If this is not the case, the module may throw an error as
- conversion from list to numpy array is not done consistently.
-
- """
- def __init__(self,Y,K,Kva=[],Kve=[],X0=None,verbose=False):
-
- """The constructor takes a phenotype vector or array Y of size n. It
- takes a kinship matrix K of size n x n. Kva and Kve can be
- computed as Kva,Kve = linalg.eigh(K) and cached. If they are
- not provided, the constructor will calculate them. X0 is an
- optional covariate matrix of size n x q, where there are q
- covariates. When this parameter is not provided, the
- constructor will set X0 to an n x 1 matrix of all ones to
- represent a mean effect.
- """
-
- if X0 is None:
- X0 = np.ones(len(Y)).reshape(len(Y),1)
- self.verbose = verbose
-
- x = True - np.isnan(Y)
- x = x.reshape(-1,)
- if not x.sum() == len(Y):
- if self.verbose: sys.stderr.write("Removing %d missing values from Y\n" % ((True - x).sum()))
- Y = Y[x]
- K = K[x,:][:,x]
- X0 = X0[x,:]
- Kva = []
- Kve = []
- self.nonmissing = x
-
- print("this K is:", K.shape, K)
-
- if len(Kva) == 0 or len(Kve) == 0:
- # if self.verbose: sys.stderr.write("Obtaining eigendecomposition for %dx%d matrix\n" % (K.shape[0],K.shape[1]) )
- begin = time.time()
- # Kva,Kve = linalg.eigh(K)
- Kva,Kve = kinship.kvakve(K)
- end = time.time()
- if self.verbose: sys.stderr.write("Total time: %0.3f\n" % (end - begin))
- print("sum(Kva),sum(Kve)=",sum(Kva),sum(Kve))
-
- self.K = K
- self.Kva = Kva
- self.Kve = Kve
- self.N = self.K.shape[0]
- self.Y = Y.reshape((self.N,1))
- self.X0 = X0
-
- if sum(self.Kva < 1e-6):
- if self.verbose: sys.stderr.write("Cleaning %d eigen values\n" % (sum(self.Kva < 0)))
- self.Kva[self.Kva < 1e-6] = 1e-6
-
- self.transform()
-
- def transform(self):
-
- """
- Computes a transformation on the phenotype vector and the covariate matrix.
- The transformation is obtained by left multiplying each parameter by the transpose of the
- eigenvector matrix of K (the kinship).
- """
-
- self.Yt = matrixMult(self.Kve.T, self.Y)
- self.X0t = matrixMult(self.Kve.T, self.X0)
- self.X0t_stack = np.hstack([self.X0t, np.ones((self.N,1))])
- self.q = self.X0t.shape[1]
-
- def getMLSoln(self,h,X):
-
- """
- Obtains the maximum-likelihood estimates for the covariate coefficients (beta),
- the total variance of the trait (sigma) and also passes intermediates that can
- be utilized in other functions. The input parameter h is a value between 0 and 1 and represents
- the heritability or the proportion of the total variance attributed to genetics. The X is the
- covariate matrix.
- """
-
- S = 1.0/(h*self.Kva + (1.0 - h))
- Xt = X.T*S
- XX = matrixMult(Xt,X)
- XX_i = inv(XX)
- beta = matrixMult(matrixMult(XX_i,Xt),self.Yt)
- Yt = self.Yt - matrixMult(X,beta)
- Q = np.dot(Yt.T*S,Yt)
- sigma = Q * 1.0 / (float(self.N) - float(X.shape[1]))
- return beta,sigma,Q,XX_i,XX
-
- def LL_brent(self,h,X=None,REML=False):
- #brent will not be bounded by the specified bracket.
- # I return a large number if we encounter h < 0 to avoid errors in LL computation during the search.
- if h < 0: return 1e6
- return -self.LL(h,X,stack=False,REML=REML)[0]
-
- def LL(self,h,X=None,stack=True,REML=False):
-
- """
- Computes the log-likelihood for a given heritability (h). If X==None, then the
- default X0t will be used. If X is set and stack=True, then X0t will be matrix concatenated with
- the input X. If stack is false, then X is used in place of X0t in the LL calculation.
- REML is computed by adding additional terms to the standard LL and can be computed by setting REML=True.
- """
-
- if X is None: X = self.X0t
- elif stack:
- self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0]
- X = self.X0t_stack
-
- n = float(self.N)
- q = float(X.shape[1])
- beta,sigma,Q,XX_i,XX = self.getMLSoln(h,X)
- LL = n*np.log(2*np.pi) + np.log(h*self.Kva + (1.0-h)).sum() + n + n*np.log(1.0/n * Q)
- LL = -0.5 * LL
-
- if REML:
- LL_REML_part = q*np.log(2.0*np.pi*sigma) + np.log(det(matrixMult(X.T,X))) - np.log(det(XX))
- LL = LL + 0.5*LL_REML_part
-
-
- LL = LL.sum()
- return LL,beta,sigma,XX_i
-
- def getMax(self,H, X=None,REML=False):
-
- """
- Helper functions for .fit(...).
- This function takes a set of LLs computed over a grid and finds possible regions
- containing a maximum. Within these regions, a Brent search is performed to find the
- optimum.
-
- """
- n = len(self.LLs)
- HOpt = []
- for i in range(1,n-2):
- if self.LLs[i-1] < self.LLs[i] and self.LLs[i] > self.LLs[i+1]:
- HOpt.append(optimize.brent(self.LL_brent,args=(X,REML),brack=(H[i-1],H[i+1])))
- if np.isnan(HOpt[-1]): HOpt[-1] = H[i-1]
- #if np.isnan(HOpt[-1]): HOpt[-1] = self.LLs[i-1]
- #if np.isnan(HOpt[-1][0]): HOpt[-1][0] = [self.LLs[i-1]]
-
- if len(HOpt) > 1:
- if self.verbose: sys.stderr.write("NOTE: Found multiple optima. Returning first...\n")
- return HOpt[0]
- elif len(HOpt) == 1: return HOpt[0]
- elif self.LLs[0] > self.LLs[n-1]: return H[0]
- else: return H[n-1]
-
-
- def fit(self,X=None,ngrids=100,REML=True):
-
- """
- Finds the maximum-likelihood solution for the heritability (h) given the current parameters.
- X can be passed and will transformed and concatenated to X0t. Otherwise, X0t is used as
- the covariate matrix.
-
- This function calculates the LLs over a grid and then uses .getMax(...) to find the optimum.
- Given this optimum, the function computes the LL and associated ML solutions.
- """
-
- if X is None: X = self.X0t
- else:
- #X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)])
- self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0]
- X = self.X0t_stack
-
- H = np.array(range(ngrids)) / float(ngrids)
- L = np.array([self.LL(h,X,stack=False,REML=REML)[0] for h in H])
- self.LLs = L
-
- hmax = self.getMax(H,X,REML)
- L,beta,sigma,betaSTDERR = self.LL(hmax,X,stack=False,REML=REML)
-
- self.H = H
- self.optH = hmax.sum()
- self.optLL = L
- self.optBeta = beta
- self.optSigma = sigma.sum()
-
- return hmax,beta,sigma,L
-
- def association(self,X,h=None,stack=True,REML=True,returnBeta=False):
- """
- Calculates association statitics for the SNPs encoded in the vector X of size n.
- If h is None, the optimal h stored in optH is used.
-
- """
- if False:
- print "X=",X
- print "h=",h
- print "q=",self.q
- print "self.Kve=",self.Kve
- print "X0t_stack=",self.X0t_stack.shape,self.X0t_stack
-
- if stack:
- # X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)])
- m = matrixMult(self.Kve.T,X)
- # print "m=",m
- m = m[:,0]
- self.X0t_stack[:,(self.q)] = m
- X = self.X0t_stack
-
- if h is None: h = self.optH
-
- L,beta,sigma,betaVAR = self.LL(h,X,stack=False,REML=REML)
- q = len(beta)
- ts,ps = self.tstat(beta[q-1],betaVAR[q-1,q-1],sigma,q)
-
- if returnBeta: return ts,ps,beta[q-1].sum(),betaVAR[q-1,q-1].sum()*sigma
- return ts,ps
-
- def tstat(self,beta,var,sigma,q,log=False):
-
- """
- Calculates a t-statistic and associated p-value given the estimate of beta and its standard error.
- This is actually an F-test, but when only one hypothesis is being performed, it reduces to a t-test.
- """
-
- ts = beta / np.sqrt(var * sigma)
- #ps = 2.0*(1.0 - stats.t.cdf(np.abs(ts), self.N-q))
- # sf == survival function - this is more accurate -- could also use logsf if the precision is not good enough
- if log:
- ps = 2.0 + (stats.t.logsf(np.abs(ts), self.N-q))
- else:
- ps = 2.0*(stats.t.sf(np.abs(ts), self.N-q))
- if not len(ts) == 1 or not len(ps) == 1:
- raise Exception("Something bad happened :(")
- return ts.sum(),ps.sum()
-
- def plotFit(self,color='b-',title=''):
-
- """
- Simple function to visualize the likelihood space. It takes the LLs
- calcualted over a grid and normalizes them by subtracting off the mean and exponentiating.
- The resulting "probabilities" are normalized to one and plotted against heritability.
- This can be seen as an approximation to the posterior distribuiton of heritability.
-
- For diagnostic purposes this lets you see if there is one distinct maximum or multiple
- and what the variance of the parameter looks like.
- """
- import matplotlib.pyplot as pl
-
- mx = self.LLs.max()
- p = np.exp(self.LLs - mx)
- p = p/p.sum()
-
- pl.plot(self.H,p,color)
- pl.xlabel("Heritability")
- pl.ylabel("Probability of data")
- pl.title(title)
-
- def meanAndVar(self):
-
- mx = self.LLs.max()
- p = np.exp(self.LLs - mx)
- p = p/p.sum()
-
- mn = (self.H * p).sum()
- vx = ((self.H - mn)**2 * p).sum()
-
- return mn,vx
-
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/optmatrix.py b/wqflask/wqflask/my_pylmm/pyLMM/optmatrix.py
deleted file mode 100644
index 5c71db6a..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/optmatrix.py
+++ /dev/null
@@ -1,55 +0,0 @@
-import sys
-import time
-import numpy as np
-from numpy.distutils.system_info import get_info;
-from scipy import linalg
-from scipy import optimize
-from scipy import stats
-
-useNumpy = None
-hasBLAS = None
-
-def matrix_initialize(useBLAS=True):
- global useNumpy # module based variable
- if useBLAS and useNumpy == None:
- print get_info('blas_opt')
- try:
- linalg.fblas
- sys.stderr.write("INFO: using linalg.fblas\n")
- useNumpy = False
- hasBLAS = True
- except AttributeError:
- sys.stderr.write("WARNING: linalg.fblas not found, using numpy.dot instead!\n")
- useNumpy = True
- else:
- sys.stderr.write("INFO: using numpy.dot\n")
- useNumpy=True
-
-def matrixMult(A,B):
- global useNumpy # module based variable
-
- if useNumpy:
- return np.dot(A,B)
-
- # If the matrices are in Fortran order then the computations will be faster
- # when using dgemm. Otherwise, the function will copy the matrix and that takes time.
- if not A.flags['F_CONTIGUOUS']:
- AA = A.T
- transA = True
- else:
- AA = A
- transA = False
-
- if not B.flags['F_CONTIGUOUS']:
- BB = B.T
- transB = True
- else:
- BB = B
- transB = False
-
- return linalg.fblas.dgemm(alpha=1.,a=AA,b=BB,trans_a=transA,trans_b=transB)
-
-def matrixMultT(M):
- # res = np.dot(W,W.T)
- # return linalg.fblas.dgemm(alpha=1.,a=M.T,b=M.T,trans_a=True,trans_b=False)
- return matrixMult(M,M.T)
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py b/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py
deleted file mode 100644
index 7b652515..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py
+++ /dev/null
@@ -1,65 +0,0 @@
-# Phenotype routines
-
-# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com)
-# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-import sys
-import numpy as np
-
-# ---- A trick to decide on the environment:
-try:
- from wqflask.my_pylmm.pyLMM import chunks
- from gn2 import uses, progress_set_func
-except ImportError:
- has_gn2=False
- import standalone as handlers
- from standalone import uses, progress_set_func
-
-progress,debug,info,mprint = uses('progress','debug','info','mprint')
-
-def remove_missing(n,y,g):
- """
- Remove missing data from matrices, make sure the genotype data has
- individuals as rows
- """
- assert(y is not None)
- assert y.shape[0] == g.shape[0],"y (n) %d, g (n,m) %s" % (y.shape[0],g.shape)
-
- y1 = y
- g1 = g
- v = np.isnan(y)
- keep = True - v
- if v.sum():
- info("runlmm.py: Cleaning the phenotype vector and genotype matrix by removing %d individuals...\n" % (v.sum()))
- y1 = y[keep]
- g1 = g[keep,:]
- n = y1.shape[0]
- return n,y1,g1,keep
-
-def remove_missing_new(n,y):
- """
- Remove missing data. Returns new n,y,keep
- """
- assert(y is not None)
- y1 = y
- v = np.isnan(y)
- keep = True - v
- if v.sum():
- info("runlmm.py: Cleaning the phenotype vector by removing %d individuals" % (v.sum()))
- y1 = y[keep]
- n = y1.shape[0]
- return n,y1,keep
-
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/plink.py b/wqflask/wqflask/my_pylmm/pyLMM/plink.py
deleted file mode 100644
index 7bd2df91..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/plink.py
+++ /dev/null
@@ -1,107 +0,0 @@
-# Plink module
-#
-# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl)
-# Some of the BED file parsing came from pylmm:
-# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com)
-
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-# According to the PLINK information
-
-# Parse a textual BIM file and return the contents as a list of tuples
-#
-# Extended variant information file accompanying a .bed binary genotype table.
-#
-# A text file with no header line, and one line per variant with the following six fields:
-#
-# Chromosome code (either an integer, or 'X'/'Y'/'XY'/'MT'; '0' indicates unknown) or name
-# Variant identifier
-# Position in morgans or centimorgans (safe to use dummy value of '0')
-# Base-pair coordinate (normally 1-based, but 0 ok; limited to 231-2)
-# Allele 1 (corresponding to clear bits in .bed; usually minor)
-# Allele 2 (corresponding to set bits in .bed; usually major)
-#
-# Allele codes can contain more than one character. Variants with negative bp coordinates are ignored by PLINK. Example
-#
-# 1 mm37-1-3125499 0 3125499 1 2
-# 1 mm37-1-3125701 0 3125701 1 2
-# 1 mm37-1-3187481 0 3187481 1 2
-
-import struct
-# import numpy as np
-
-def readbim(fn):
- res = []
- for line in open(fn):
- list = line.split()
- if len([True for e in list if e == 'nan']) == 0:
- res.append( (list[0],list[1],int(list[2]),int(list[3]),int(list[4]),int(list[5])) )
- else:
- res.append( (list[0],list[1],list[2],float('nan'),float('nan'),float('nan')) )
- return res
-
-# .bed (PLINK binary biallelic genotype table)
-#
-# Primary representation of genotype calls at biallelic variants. Must
-# be accompanied by .bim and .fam files. Basically contains num SNP
-# blocks containing IND (compressed 4 IND into a byte)
-#
-# Since it is a biallelic format it supports for every individual
-# whether the first allele is homozygous (b00), the second allele is
-# homozygous (b11), it is heterozygous (b10) or that it is missing
-# (b01).
-
-# http://pngu.mgh.harvard.edu/~purcell/plink2/formats.html#bed
-# http://pngu.mgh.harvard.edu/~purcell/plink2/formats.html#fam
-# http://pngu.mgh.harvard.edu/~purcell/plink2/formats.html#bim
-
-def readbed(fn,inds,encoding,func=None):
-
- # For every SNP block fetch the individual genotypes using values
- # 0.0 and 1.0 for homozygous and 0.5 for heterozygous alleles
- def fetchGenotypes(X):
- # D = { \
- # '00': 0.0, \
- # '10': 0.5, \
- # '11': 1.0, \
- # '01': float('nan') \
- # }
-
- Didx = { '00': 0, '10': 1, '11': 2, '01': 3 }
- G = []
- for x in X:
- if not len(x) == 10:
- xx = x[2:]
- x = '0b' + '0'*(8 - len(xx)) + xx
- a,b,c,d = (x[8:],x[6:8],x[4:6],x[2:4])
- L = [encoding[Didx[y]] for y in [a,b,c,d]]
- G += L
- G = G[:inds]
- # G = np.array(G)
- return G
-
- bytes = inds / 4 + (inds % 4 and 1 or 0)
- format = 'c'*bytes
- count = 0
- with open(fn,'rb') as f:
- magic = f.read(3)
- assert( ":".join("{:02x}".format(ord(c)) for c in magic) == "6c:1b:01")
- while True:
- count += 1
- X = f.read(bytes)
- if not X:
- return(count-1)
- XX = [bin(ord(x)) for x in struct.unpack(format,X)]
- xs = fetchGenotypes(XX)
- func(count,xs)
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
deleted file mode 100644
index 6b241cd6..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ /dev/null
@@ -1,229 +0,0 @@
-# This is the LMM runner that calls the possible methods using command line
-# switches. It acts as a multiplexer where all the invocation complexity
-# is kept outside the main LMM routines.
-#
-# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-from optparse import OptionParser
-import sys
-import tsvreader
-import numpy as np
-
-# Add local dir to PYTHONPATH
-import os
-cwd = os.path.dirname(__file__)
-if sys.path[0] != cwd:
- sys.path.insert(1,cwd)
-
-# pylmm modules
-from lmm import gn2_load_redis, gn2_load_redis_iter, calculate_kinship_new, run_gwas
-from kinship import kinship, kinship_full
-import genotype
-import phenotype
-from standalone import uses
-
-progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal')
-
-usage = """
-python runlmm.py [options] command
-
- runlmm.py processing multiplexer reads standardised input formats
- and calls the different routines (writes to stdout)
-
- Current commands are:
-
- parse : only parse input files
- redis : use Redis to call into GN2
- kinship : calculate (new) kinship matrix
-
- try --help for more information
-"""
-
-
-parser = OptionParser(usage=usage)
-# parser.add_option("-f", "--file", dest="input file",
-# help="In", metavar="FILE")
-parser.add_option("--kinship",dest="kinship",
- help="Kinship file format 1.0")
-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("--genotype-normalization",
- action="store_true", dest="genotype_normalization", default=False,
- help="Force genotype normalization")
-parser.add_option("--remove-missing-phenotypes",
- action="store_true", dest="remove_missing_phenotypes", default=False,
- help="Remove missing phenotypes")
-parser.add_option("-q", "--quiet",
- action="store_false", dest="verbose", default=True,
- help="don't print status messages to stdout")
-parser.add_option("--blas", action="store_true", default=False, dest="useBLAS", help="Use BLAS instead of numpy matrix multiplication")
-parser.add_option("-t", "--threads",
- type="int", dest="numThreads",
- help="Threads to use")
-parser.add_option("--saveKvaKve",
- action="store_true", dest="saveKvaKve", default=False,
- help="Testing mode")
-parser.add_option("--test",
- action="store_true", dest="testing", default=False,
- help="Testing mode")
-parser.add_option("--test-kinship",
- action="store_true", dest="test_kinship", default=False,
- help="Testing mode for Kinship calculation")
-
-(options, args) = parser.parse_args()
-
-if len(args) != 1:
- print usage
- sys.exit(1)
-
-cmd = args[0]
-print "Command: ",cmd
-
-k = None
-y = None
-g = None
-
-if options.kinship:
- k = tsvreader.kinship(options.kinship)
- print k.shape
-
-if options.pheno:
- y = tsvreader.pheno(options.pheno)
- print y.shape
-
-if options.geno and cmd != 'iterator':
- g = tsvreader.geno(options.geno)
- print g.shape
-
-def check_results(ps,ts):
- print np.array(ps)
- print len(ps),sum(ps)
- p1 = round(ps[0],4)
- p2 = round(ps[-1],4)
- if options.geno == 'data/small.geno':
- info("Validating results for "+options.geno)
- assert p1==0.7387, "p1=%f" % p1
- assert p2==0.7387, "p2=%f" % p2
- if options.geno == 'data/small_na.geno':
- info("Validating results for "+options.geno)
- assert p1==0.062, "p1=%f" % p1
- assert p2==0.062, "p2=%f" % p2
- if options.geno == 'data/test8000.geno':
- info("Validating results for "+options.geno)
- assert round(sum(ps)) == 4070
- assert len(ps) == 8000
- info("Run completed")
-
-if y is not None:
- n = y.shape[0]
-
-if cmd == 'run':
- if options.remove_missing_phenotypes:
- raise Exception('Can not use --remove-missing-phenotypes with LMM2')
- n = len(y)
- m = g.shape[1]
- ps, ts = run_gwas('other',n,m,k,y,g) # <--- pass in geno by SNP
- check_results(ps,ts)
-elif cmd == 'iterator':
- if options.remove_missing_phenotypes:
- raise Exception('Can not use --remove-missing-phenotypes with LMM2')
- geno_iterator = tsvreader.geno_iter(options.geno)
- ps, ts = gn2_load_redis_iter('testrun_iter','other',k,y,geno_iterator)
- check_results(ps,ts)
-elif cmd == 'redis_new':
- # The main difference between redis_new and redis is that missing
- # phenotypes are handled by the first
- if options.remove_missing_phenotypes:
- raise Exception('Can not use --remove-missing-phenotypes with LMM2')
- Y = y
- G = g
- print "Original G",G.shape, "\n", G
- # gt = G.T
- # G = None
- ps, ts = gn2_load_redis('testrun','other',k,Y,G,new_code=True)
- check_results(ps,ts)
-elif cmd == 'redis':
- # Emulating the redis setup of GN2
- G = g
- print "Original G",G.shape, "\n", G
- if y is not None and options.remove_missing_phenotypes:
- gnt = np.array(g).T
- n,Y,g,keep = phenotype.remove_missing(n,y,gnt)
- G = g.T
- print "Removed missing phenotypes",G.shape, "\n", G
- else:
- Y = y
- if options.maf_normalization:
- G = np.apply_along_axis( genotype.replace_missing_with_MAF, axis=0, arr=g )
- print "MAF replacements: \n",G
- if options.genotype_normalization:
- G = np.apply_along_axis( genotype.normalize, axis=1, arr=G)
- g = None
- gnt = None
-
- # gt = G.T
- # G = None
- ps, ts = gn2_load_redis('testrun','other',k,Y,G, new_code=False)
- check_results(ps,ts)
-elif cmd == 'kinship':
- G = g
- print "Original G",G.shape, "\n", G
- if y != None and options.remove_missing_phenotypes:
- gnt = np.array(g).T
- n,Y,g,keep = phenotype.remove_missing(n,y,g.T)
- G = g.T
- print "Removed missing phenotypes",G.shape, "\n", G
- if options.maf_normalization:
- G = np.apply_along_axis( genotype.replace_missing_with_MAF, axis=0, arr=g )
- print "MAF replacements: \n",G
- if options.genotype_normalization:
- G = np.apply_along_axis( genotype.normalize, axis=1, arr=G)
- g = None
- gnt = None
-
- if options.test_kinship:
- K = kinship_full(np.copy(G))
- print "Genotype",G.shape, "\n", G
- print "first Kinship method",K.shape,"\n",K
- k1 = round(K[0][0],4)
- K2,G = calculate_kinship_new(np.copy(G))
- print "Genotype",G.shape, "\n", G
- print "GN2 Kinship method",K2.shape,"\n",K2
- k2 = round(K2[0][0],4)
-
- print "Genotype",G.shape, "\n", G
- K3 = kinship(G)
- print "third Kinship method",K3.shape,"\n",K3
- sys.stderr.write(options.geno+"\n")
- k3 = round(K3[0][0],4)
- if options.geno == 'data/small.geno':
- assert k1==0.8333, "k1=%f" % k1
- assert k2==0.9375, "k2=%f" % k2
- assert k3==0.9375, "k3=%f" % k3
- if options.geno == 'data/small_na.geno':
- assert k1==0.8333, "k1=%f" % k1
- assert k2==0.7172, "k2=%f" % k2
- assert k3==0.7172, "k3=%f" % k3
- if options.geno == 'data/test8000.geno':
- assert k3==1.4352, "k3=%f" % k3
-
-else:
- fatal("Doing nothing")
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/standalone.py b/wqflask/wqflask/my_pylmm/pyLMM/standalone.py
deleted file mode 100644
index 40b2021d..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/standalone.py
+++ /dev/null
@@ -1,110 +0,0 @@
-# Standalone specific methods and callback handler
-#
-# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# Set the log level with
-#
-# logging.basicConfig(level=logging.DEBUG)
-
-from __future__ import absolute_import, print_function, division
-
-import numpy as np
-import sys
-import logging
-
-# logger = logging.getLogger(__name__)
-logger = logging.getLogger('lmm2')
-logging.basicConfig(level=logging.DEBUG)
-np.set_printoptions(precision=3,suppress=True)
-
-progress_location = None
-progress_current = None
-progress_prev_perc = None
-
-def progress_default_func(location,count,total):
- global progress_current
- value = round(count*100.0/total)
- progress_current = value
-
-progress_func = progress_default_func
-
-def progress_set_func(func):
- global progress_func
- progress_func = func
-
-def progress(location, count, total):
- global progress_location
- global progress_prev_perc
-
- perc = round(count*100.0/total)
- if perc != progress_prev_perc and (location != progress_location or perc > 98 or perc > progress_prev_perc + 5):
- progress_func(location, count, total)
- logger.info("Progress: %s %d%%" % (location,perc))
- progress_location = location
- progress_prev_perc = perc
-
-def mprint(msg,data):
- """
- Array/matrix print function
- """
- m = np.array(data)
- if m.ndim == 1:
- print(msg,m.shape,"=\n",m[0:3]," ... ",m[-3:])
- if m.ndim == 2:
- print(msg,m.shape,"=\n[",
- m[0][0:3]," ... ",m[0][-3:],"\n ",
- m[1][0:3]," ... ",m[1][-3:],"\n ...\n ",
- m[-2][0:3]," ... ",m[-2][-3:],"\n ",
- m[-1][0:3]," ... ",m[-1][-3:],"]")
-
-def fatal(msg):
- logger.critical(msg)
- raise Exception(msg)
-
-def callbacks():
- return dict(
- write = sys.stdout.write,
- writeln = print,
- debug = logger.debug,
- info = logger.info,
- warning = logger.warning,
- error = logger.error,
- critical = logger.critical,
- fatal = fatal,
- progress = progress,
- mprint = mprint
- )
-
-def uses(*funcs):
- """
- Some sugar
- """
- return [callbacks()[func] for func in funcs]
-
-# ----- Minor test cases:
-
-if __name__ == '__main__':
- # logging.basicConfig(level=logging.DEBUG)
- logging.debug("Test %i" % (1))
- d = callbacks()['debug']
- d("TEST")
- wrln = callbacks()['writeln']
- wrln("Hello %i" % 34)
- progress = callbacks()['progress']
- progress("I am half way",50,100)
- list = [0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
- 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
- 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
- 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15,
- 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15]
- mprint("list",list)
- matrix = [[1,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
- [2,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
- [3,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
- [4,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
- [5,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15],
- [6,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15]]
- mprint("matrix",matrix)
- ix,dx = uses("info","debug")
- ix("ix")
- dx("dx")
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/temp_data.py b/wqflask/wqflask/my_pylmm/pyLMM/temp_data.py
deleted file mode 100644
index 004d45c6..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/temp_data.py
+++ /dev/null
@@ -1,25 +0,0 @@
-from __future__ import print_function, division, absolute_import
-from redis import Redis
-
-import simplejson as json
-
-class TempData(object):
-
- def __init__(self, temp_uuid):
- self.temp_uuid = temp_uuid
- self.redis = Redis()
- self.key = "tempdata:{}".format(self.temp_uuid)
-
- def store(self, field, value):
- self.redis.hset(self.key, field, value)
- self.redis.expire(self.key, 60*15) # Expire in 15 minutes
-
- def get_all(self):
- return self.redis.hgetall(self.key)
-
-
-if __name__ == "__main__":
- redis = Redis()
- for key in redis.keys():
- for field in redis.hkeys(key):
- print("{}.{}={}".format(key, field, redis.hget(key, field)))
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py b/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py
deleted file mode 100644
index 66b34ee2..00000000
--- a/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py
+++ /dev/null
@@ -1,122 +0,0 @@
-# Standard file readers
-#
-# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Affero General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
-
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Affero General Public License for more details.
-
-# You should have received a copy of the GNU Affero General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-import sys
-import os
-import numpy as np
-import csv
-
-def kinship(fn):
- K1 = []
- print fn
- with open(fn,'r') as tsvin:
- assert(tsvin.readline().strip() == "# Kinship format version 1.0")
- tsvin.readline()
- tsvin.readline()
- tsv = csv.reader(tsvin, delimiter='\t')
- for row in tsv:
- ns = np.genfromtxt(row[1:])
- K1.append(ns) # <--- slow
- K = np.array(K1)
- return K
-
-def pheno(fn):
- Y1 = []
- print fn
- with open(fn,'r') as tsvin:
- assert(tsvin.readline().strip() == "# Phenotype format version 1.0")
- tsvin.readline()
- tsvin.readline()
- tsvin.readline()
- tsv = csv.reader(tsvin, delimiter='\t')
- for row in tsv:
- ns = np.genfromtxt(row[1:])
- Y1.append(ns) # <--- slow
- Y = np.array(Y1)
- return Y
-
-def geno(fn):
- G1 = []
- hab_mapper = {'A':0,'H':1,'B':2,'-':3}
- pylmm_mapper = [ 0.0, 0.5, 1.0, float('nan') ]
-
- print fn
- with open(fn,'r') as tsvin:
- line = tsvin.readline().strip()
- assert line == "# Genotype format version 1.0", line
- tsvin.readline()
- tsvin.readline()
- tsvin.readline()
- tsvin.readline()
- tsv = csv.reader(tsvin, delimiter='\t')
- for row in tsv:
- # print(row)
- id = row[0]
- gs = list(row[1])
- # print id,gs
- gs2 = [pylmm_mapper[hab_mapper[g]] for g in gs]
- # print id,gs2
- # ns = np.genfromtxt(row[1:])
- G1.append(gs2) # <--- slow
- G = np.array(G1)
- return G
-
-def geno(fn):
- G1 = []
- for id,values in geno_iter(fn):
- G1.append(values) # <--- slow
- G = np.array(G1)
- return G
-
-def geno_callback(fn,func):
- hab_mapper = {'A':0,'H':1,'B':2,'-':3}
- pylmm_mapper = [ 0.0, 0.5, 1.0, float('nan') ]
-
- print fn
- with open(fn,'r') as tsvin:
- assert(tsvin.readline().strip() == "# Genotype format version 1.0")
- tsvin.readline()
- tsvin.readline()
- tsvin.readline()
- tsvin.readline()
- tsv = csv.reader(tsvin, delimiter='\t')
- for row in tsv:
- id = row[0]
- gs = list(row[1])
- gs2 = [pylmm_mapper[hab_mapper[g]] for g in gs]
- func(id,gs2)
-
-def geno_iter(fn):
- """
- Yield a tuple of snpid and values
- """
- hab_mapper = {'A':0,'H':1,'B':2,'-':3}
- pylmm_mapper = [ 0.0, 0.5, 1.0, float('nan') ]
-
- print fn
- with open(fn,'r') as tsvin:
- assert(tsvin.readline().strip() == "# Genotype format version 1.0")
- tsvin.readline()
- tsvin.readline()
- tsvin.readline()
- tsvin.readline()
- tsv = csv.reader(tsvin, delimiter='\t')
- for row in tsv:
- id = row[0]
- gs = list(row[1])
- gs2 = [pylmm_mapper[hab_mapper[g]] for g in gs]
- yield (id,gs2)
diff --git a/wqflask/wqflask/templates/collections/view.html b/wqflask/wqflask/templates/collections/view.html
index fc1edf2a..e41ade7c 100755
--- a/wqflask/wqflask/templates/collections/view.html
+++ b/wqflask/wqflask/templates/collections/view.html
@@ -59,6 +59,9 @@
<button class="btn btn-default" id="select_all"><span class="glyphicon glyphicon-ok"></span> Select All</button>
<button class="btn btn-default" id="deselect_all"><span class="glyphicon glyphicon-remove"></span> Deselect All</button>
<button class="btn btn-default" id="invert"><span class="glyphicon glyphicon-resize-vertical"></span> Invert</button>
+ <button class="btn" id="add" disabled="disabled"><i class="icon-plus-sign"></i> Add Record to Other Collection</button>
+ <button class="btn" id="remove" disabled="disabled"><i class="icon-minus-sign"></i> Remove Record</button>
+ <button class="btn btn-primary pull-right"><i class="icon-download icon-white"></i> Download Table</button>
<br />
<br />
<table class="table table-hover table-striped" id='trait_table'>
generated by cgit v1.2.3 (git 2.44.0) at 2024-12-11 12:41:00 +0000