From 3325184b1dd310619626dd31852ab84cae6dc7fc Mon Sep 17 00:00:00 2001 From: Zachary Sloan Date: Tue, 8 Oct 2013 12:13:56 -0500 Subject: Did some work with interval mapping page; will use qtl reaper to do the mapping, since it is reliable/fast and avoids us having to rewrite from scratch using something like r/qtl --- wqflask/base/data_set.py | 20 ++++---- wqflask/maintenance/quick_search_table.py | 4 +- .../wqflask/interval_mapping/interval_mapping.py | 53 ++++++++++++++++++++++ 3 files changed, 67 insertions(+), 10 deletions(-) (limited to 'wqflask') diff --git a/wqflask/base/data_set.py b/wqflask/base/data_set.py index 96e04df0..befbd518 100755 --- a/wqflask/base/data_set.py +++ b/wqflask/base/data_set.py @@ -168,13 +168,13 @@ class Markers(object): for marker, p_value in itertools.izip(self.markers, p_values): marker['p_value'] = p_value - if marker['p_value'] == 0: - marker['lod_score'] = 0 - marker['lrs_value'] = 0 - else: - marker['lod_score'] = -math.log10(marker['p_value']) - #Using -log(p) for the LRS; need to ask Rob how he wants to get LRS from p-values - marker['lrs_value'] = -math.log10(marker['p_value']) * 4.61 + if math.isnan(marker['p_value']): + print("p_value is:", marker['p_value']) + marker['lod_score'] = -math.log10(marker['p_value']) + #Using -log(p) for the LRS; need to ask Rob how he wants to get LRS from p-values + marker['lrs_value'] = -math.log10(marker['p_value']) * 4.61 + + class HumanMarkers(Markers): @@ -189,6 +189,8 @@ class HumanMarkers(Markers): marker['name'] = splat[1] marker['Mb'] = float(splat[3]) / 1000000 self.markers.append(marker) + + #print("markers is: ", pf(self.markers)) def add_pvalues(self, p_values): @@ -315,10 +317,10 @@ class DatasetGroup(object): #determine default genotype object if self.incparentsf1 and genotype_1.type != "intercross": - genotype = genotype_2 + self.genotype = genotype_2 else: self.incparentsf1 = 0 - genotype = genotype_1 + self.genotype = genotype_1 self.samplelist = list(genotype.prgy) diff --git a/wqflask/maintenance/quick_search_table.py b/wqflask/maintenance/quick_search_table.py index 9cd792ef..23bd505c 100644 --- a/wqflask/maintenance/quick_search_table.py +++ b/wqflask/maintenance/quick_search_table.py @@ -11,8 +11,10 @@ each trait, its dataset, and several columns determined by its trait type (pheno """ -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function +# We do this here so we can use zach_settings +# Not to avoid other absoulte_imports import sys sys.path.append("../../..") diff --git a/wqflask/wqflask/interval_mapping/interval_mapping.py b/wqflask/wqflask/interval_mapping/interval_mapping.py index 48e8018e..5d660224 100644 --- a/wqflask/wqflask/interval_mapping/interval_mapping.py +++ b/wqflask/wqflask/interval_mapping/interval_mapping.py @@ -12,6 +12,7 @@ import collections import numpy as np from scipy import linalg +import rpy2.robjects import simplejson as json @@ -83,6 +84,28 @@ class IntervalMapping(object): """Generates qtl results for plotting interval map""" self.dataset.group.get_markers() + self.dataset.read_genotype_file() + + samples, values, variances = self.trait.export_informative() + if self.control_locus: + if self.weighted_regression: + qtl_result = self.dataset.genotype.regression(strains = samples, + trait = values, + variance = variances, + control = self.control_locus) + else: + qtl_result = self.dataset.genotype.regression(strains = samples, + trait = values, + control = self.control_locus) + else: + if self.weighted_regression: + qtl_result = self.dataset.genotype.regression(strains = samples, + trait = values, + variance = variances) + else: + qtl_result = self.dataset.genotype.regression(strains = samples, + trait = values) + pheno_vector = np.array([val == "x" and np.nan or float(val) for val in self.vals]) @@ -108,6 +131,36 @@ class IntervalMapping(object): self.qtl_results = self.dataset.group.markers.markers + #def gen_qtl_results_2(self, tempdata): + # """Generates qtl results for plotting interval map""" + # + # self.dataset.group.get_markers() + # self.dataset.read_genotype_file() + # + # pheno_vector = np.array([val == "x" and np.nan or float(val) for val in self.vals]) + # + # #if self.dataset.group.species == "human": + # # p_values, t_stats = self.gen_human_results(pheno_vector, tempdata) + # #else: + # 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 + # + # t_stats, p_values = lmm.run( + # pheno_vector, + # genotype_matrix, + # restricted_max_likelihood=True, + # refit=False, + # temp_data=tempdata + # ) + # + # self.dataset.group.markers.add_pvalues(p_values) + # + # self.qtl_results = self.dataset.group.markers.markers + def identify_empty_samples(self): no_val_samples = [] -- cgit v1.2.3