diff options
Diffstat (limited to 'wqflask')
23 files changed, 325 insertions, 4031 deletions
diff --git a/wqflask/base/data_set.py b/wqflask/base/data_set.py index 489bd374..1cd57b4b 100755 --- a/wqflask/base/data_set.py +++ b/wqflask/base/data_set.py @@ -42,7 +42,7 @@ from base import species from dbFunction import webqtlDatabaseFunction from utility import webqtlUtil from utility.benchmark import Bench -from wqflask.my_pylmm.pyLMM import chunks +from utility import chunks from maintenance import get_group_samplelists diff --git a/wqflask/wqflask/my_pylmm/pyLMM/chunks.py b/wqflask/utility/chunks.py index 9565fb96..9565fb96 100644 --- a/wqflask/wqflask/my_pylmm/pyLMM/chunks.py +++ b/wqflask/utility/chunks.py 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/interval_mapping/interval_mapping.py b/wqflask/wqflask/interval_mapping/interval_mapping.py index 5511826a..1cd3fc80 100755 --- a/wqflask/wqflask/interval_mapping/interval_mapping.py +++ b/wqflask/wqflask/interval_mapping/interval_mapping.py @@ -24,9 +24,9 @@ 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 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 diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py index 49521bd6..76d05bd8 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/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) |