From 5849fe30c792a45d872c505b01d4705ce3711710 Mon Sep 17 00:00:00 2001 From: Zachary Sloan Date: Fri, 8 Mar 2013 21:15:07 +0000 Subject: Refactored GWAS function in lmm.py --- .../wqflask/marker_regression/marker_regression.py | 4 +- wqflask/wqflask/my_pylmm/pyLMM/lmm.py | 223 +++++++-------------- .../new/javascript/show_trait_mapping_tools.js | 70 ++++++- 3 files changed, 145 insertions(+), 152 deletions(-) diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py index 4ddc89c6..4f6e8e1a 100755 --- a/wqflask/wqflask/marker_regression/marker_regression.py +++ b/wqflask/wqflask/marker_regression/marker_regression.py @@ -491,12 +491,12 @@ class MarkerRegression(object): with Bench("LMM_ob fitting"): lmm_ob.fit() - + with Bench("Doing gwas"): t_stats, p_values = lmm.GWAS(pheno_vector, genotype_matrix, kinship_matrix, - REML=True, + restricted_max_likelihood=True, refit=False, temp_data=self.temp_data) diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py index 12f7c2ea..6c101ba1 100644 --- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py +++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py @@ -95,179 +95,104 @@ def calculate_kinship(genotype_matrix, temp_data): kinship_matrix = np.dot(genotype_matrix,genotype_matrix.T) * 1.0/float(m) return kinship_matrix -def GWAS(Y, X, K, Kva=[], Kve=[], X0=None, REML=True, refit=False, temp_data=None): +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 - - Y - n x 1 phenotype vector - X - n x m SNP 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 + 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 + """ - n = X.shape[0] - m = X.shape[1] + if kinship_eigen_vals == None: + kinship_eigen_vals = [] + if kinship_eigen_vectors= == None: + kinship_eigen_vectors = [] + + n = genotype_matrix.shape[0] + m = genotype_matrix.shape[1] - if X0 == None: X0 = np.ones((n,1)) + if covariate_matrix == None: + covariate_matrix = np.ones((n,1)) - # Remove missing values in Y and adjust associated parameters - v = np.isnan(Y) + # Remove missing values in pheno_vector and adjust associated parameters + v = np.isnan(pheno_vector) if v.sum(): keep = True - v - Y = Y[keep] - X = X[keep,:] - X0 = X0[keep,:] - K = K[keep,:][:,keep] - Kva = [] - Kve = [] - - L = LMM(Y,K,Kva,Kve,X0) - if not refit: L.fit() - - PS = [] - TS = [] + 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 = [] for counter in range(m): - x = X[:,counter].reshape((n,1)) + 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: - PS.append(np.nan) - TS.append(np.nan) + p_values.append(np.nan) + t_statistics.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) + 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 = lmm_ob_2.association(xs, REML=restricted_max_likelihood) else: if x.var() == 0: - PS.append(np.nan) - TS.append(np.nan) + p_values.append(np.nan) + t_statistics.append(np.nan) continue - if refit: L.fit(X=x) - ts,ps = L.association(x,REML=REML) + if refit: + lmm_ob.fit(X=x) + ts,ps = lmm_ob.association(x,REML=restricted_max_likelihood) percent_complete = 45 + int(round((counter/m)*55)) print("Percent complete: ", percent_complete) temp_data.store("percent_complete", percent_complete) - PS.append(ps) - TS.append(ts) - - return TS,PS - -#def GWAS(pheno_vector, -# genotype_matrix, -# kinship_matrix, -# kinship_eigenvals=None, -# kinship_eigenvectors=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 -# -# Y - n x 1 phenotype vector -# X - n x m SNP 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 -# -# """ -# -# assert temp_data, "You forgot to pass in temp_data" -# -# if kinship_eigenvals == None: -# kinship_eigenvals = [] -# if kinship_eigenvectors == None: -# kinship_eigenvectors = [] -# -# n = genotype_matrix.shape[0] -# m = genotype_matrix.shape[1] -# -# if covariate_matrix == None: -# covariate_matrix = np.ones((n,1)) -# -# # Remove missing values in Y and adjust associated parameters -# pheno_not_number = np.isnan(pheno_vector) -# if pheno_not_number.sum(): -# keep = True - pheno_not_number -# pheno_vector = pheno_vector[keep] -# genotype_matrix = genotype_matrix[keep,:] -# covariate_matrix = covariate_matrix[keep,:] -# kinship_matrix = kinship_matrix[keep,:][:,keep] -# kinship_eigenvals = [] -# kinship_eigenvectors = [] -# -# lmm_ob = LMM(pheno_vector, -# kinship_matrix, -# kinship_eigenvals, -# kinship_eigenvectors, -# covariate_matrix) -# if not refit: -# lmm_ob.fit() -# -# p_value_matrix = [] -# t_stats_matrix = [] -# -# for counter in range(m): -# #pheno_vector_2 = geno_vector[:, counter] -# #x = pheno_vector_2.reshape((n,1)) -# 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_value_matrix.append(np.nan) -# t_stats_matrix.append(np.nan) -# continue -# -# pheno_vector_2 = pheno_vector[keep] -# covariate_matrix_2 = covariate_matrix[keep,:] -# kinship_matrix_2 = kinship_matrix[keep,:][:,keep] -# lmm_ob_2 = LMM(pheno_vector, kinship_matrix, covariate_matrix=covariate_matrix) -# if refit: -# lmm_ob_2.fit(X=xs) -# else: -# lmm_ob_2.fit() -# t_stats, p_values = lmm_ob_2.association(xs, REML=restricted_max_likelihood) -# else: -# if x.var() == 0: -# p_value_matrix.append(np.nan) -# t_stats_matrix.append(np.nan) -# continue -# -# if refit: -# lmm_ob.fit(X=x) -# t_stats,p_values = lmm_ob.association(x, REML=restricted_max_likelihood) -# -# p_value_matrix.append(p_values) -# t_stats_matrix.append(t_stats) -# -# percent_complete = 45 + int(round((counter/m)*55)) -# print("Percent complete: ", percent_complete) -# temp_data.store("percent_complete", percent_complete) -# -# return p_value_matrix, t_stats_matrix + p_values.append(ps) + t_statistics.append(ts) + + return t_statistics,p_values + class LMM: diff --git a/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js b/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js index c8b0aa7b..78459692 100644 --- a/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js +++ b/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js @@ -1,6 +1,74 @@ // Generated by CoffeeScript 1.4.0 (function() { - + $(function() { + var composite_mapping_fields, get_progress, submit_special, toggle_enable_disable, + _this = this; + submit_special = function() { + var url; + console.log("In submit_special"); + console.log("this is:", this); + console.log("$(this) is:", $(this)); + url = $(this).data("url"); + console.log("url is:", url); + $("#trait_data_form").attr("action", url); + return $("#trait_data_form").submit(); + }; + get_progress = function() { + var params, params_str, temp_uuid, url, + _this = this; + console.log("temp_uuid:", $("#temp_uuid").val()); + temp_uuid = $("#temp_uuid").val(); + params = { + key: temp_uuid + }; + params_str = $.param(params); + url = "/get_temp_data?" + params_str; + console.log("url:", url); + $.ajax({ + type: "GET", + url: url, + success: function(progress_data) { + console.log("in get_progress data:", progress_data); + console.log(progress_data['percent_complete'] + "%"); + return $('#marker_regression_progress').css("width", progress_data['percent_complete'] + "%"); + } + }); + return false; + }; + $("#marker_regression").click(function() { + var form_data, url; + $("#progress_bar_container").modal(); + url = "/marker_regression"; + form_data = $('#trait_data_form').serialize(); + console.log("form_data is:", form_data); + $.ajax({ + type: "POST", + url: url, + data: form_data, + success: function(data) { + clearInterval(_this.my_timer); + $('#progress_bar_container').modal('hide'); + return $("body").html(data); + } + }); + console.log("settingInterval"); + _this.my_timer = setInterval(get_progress, 1000); + return false; + }); + composite_mapping_fields = function() { + return $(".composite_fields").toggle(); + }; + $("#use_composite_choice").change(composite_mapping_fields); + toggle_enable_disable = function(elem) { + return $(elem).prop("disabled", !$(elem).prop("disabled")); + }; + $("#choose_closet_control").change(function() { + return toggle_enable_disable("#control_locus"); + }); + return $("#display_all_lrs").change(function() { + return toggle_enable_disable("#suggestive_lrs"); + }); + }); }).call(this); -- cgit v1.2.3