diff options
| -rw-r--r-- | wqflask/wqflask/marker_regression/marker_regression.py | 16 | ||||
| -rw-r--r-- | wqflask/wqflask/marker_regression/qtlreaper_mapping.py | 186 |
2 files changed, 101 insertions, 101 deletions
diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py index f340e309..200f2207 100644 --- a/wqflask/wqflask/marker_regression/marker_regression.py +++ b/wqflask/wqflask/marker_regression/marker_regression.py @@ -201,14 +201,14 @@ class MarkerRegression(object): self.control_marker = start_vars['control_marker'] self.do_control = start_vars['do_control'] logger.info("Running qtlreaper") - results, self.json_data, self.perm_output, self.suggestive, self.significant, self.bootstrap_results = qtlreaper_mapping.gen_reaper_results(self.this_trait, - self.dataset, - self.samples, - self.json_data, - self.num_perm, - self.bootCheck, - self.num_bootstrap, - self.do_control, + results, self.json_data, self.perm_output, self.suggestive, self.significant, self.bootstrap_results = qtlreaper_mapping.gen_reaper_results(self.this_trait, + self.dataset, + self.samples, + self.json_data, + self.num_perm, + self.bootCheck, + self.num_bootstrap, + self.do_control, self.control_marker, self.manhattan_plot) elif self.mapping_method == "plink": diff --git a/wqflask/wqflask/marker_regression/qtlreaper_mapping.py b/wqflask/wqflask/marker_regression/qtlreaper_mapping.py index 568991f7..b072584c 100644 --- a/wqflask/wqflask/marker_regression/qtlreaper_mapping.py +++ b/wqflask/wqflask/marker_regression/qtlreaper_mapping.py @@ -1,93 +1,93 @@ -def gen_reaper_results(this_trait, dataset, samples_before, json_data, num_perm, bootCheck, num_bootstrap, do_control, control_marker, manhattan_plot):
- genotype = dataset.group.read_genotype_file()
-
- if manhattan_plot != True:
- genotype = genotype.addinterval()
-
- samples, values, variances, sample_aliases = this_trait.export_informative()
-
- trimmed_samples = []
- trimmed_values = []
- for i in range(0, len(samples)):
- if this_trait.data[samples[i]].name in samples_before:
- trimmed_samples.append(samples[i])
- trimmed_values.append(values[i])
-
- perm_output = []
- bootstrap_results = []
-
- if num_perm < 100:
- suggestive = 0
- significant = 0
- else:
- perm_output = genotype.permutation(strains = trimmed_samples, trait = trimmed_values, nperm=num_perm)
- suggestive = perm_output[int(num_perm*0.37-1)]
- significant = perm_output[int(num_perm*0.95-1)]
- highly_significant = perm_output[int(num_perm*0.99-1)]
-
- json_data['suggestive'] = suggestive
- json_data['significant'] = significant
-
- if control_marker != "" and do_control == "true":
- reaper_results = genotype.regression(strains = trimmed_samples,
- trait = trimmed_values,
- control = str(control_marker))
- if bootCheck:
- control_geno = []
- control_geno2 = []
- _FIND = 0
- for _chr in genotype:
- for _locus in _chr:
- if _locus.name == control_marker:
- control_geno2 = _locus.genotype
- _FIND = 1
- break
- if _FIND:
- break
- if control_geno2:
- _prgy = list(genotype.prgy)
- for _strain in trimmed_samples:
- _idx = _prgy.index(_strain)
- control_geno.append(control_geno2[_idx])
-
- bootstrap_results = genotype.bootstrap(strains = trimmed_samples,
- trait = trimmed_values,
- control = control_geno,
- nboot = num_bootstrap)
- else:
- reaper_results = genotype.regression(strains = trimmed_samples,
- trait = trimmed_values)
-
- if bootCheck:
- bootstrap_results = genotype.bootstrap(strains = trimmed_samples,
- trait = trimmed_values,
- nboot = num_bootstrap)
-
- json_data['chr'] = []
- json_data['pos'] = []
- json_data['lod.hk'] = []
- json_data['markernames'] = []
- #if self.additive:
- # self.json_data['additive'] = []
-
- #Need to convert the QTL objects that qtl reaper returns into a json serializable dictionary
- qtl_results = []
- for qtl in reaper_results:
- reaper_locus = qtl.locus
- #ZS: Convert chr to int
- converted_chr = reaper_locus.chr
- if reaper_locus.chr != "X" and reaper_locus.chr != "X/Y":
- converted_chr = int(reaper_locus.chr)
- json_data['chr'].append(converted_chr)
- json_data['pos'].append(reaper_locus.Mb)
- json_data['lod.hk'].append(qtl.lrs)
- json_data['markernames'].append(reaper_locus.name)
- #if self.additive:
- # self.json_data['additive'].append(qtl.additive)
- locus = {"name":reaper_locus.name, "chr":reaper_locus.chr, "cM":reaper_locus.cM, "Mb":reaper_locus.Mb}
- qtl = {"lrs_value": qtl.lrs, "chr":converted_chr, "Mb":reaper_locus.Mb,
- "cM":reaper_locus.cM, "name":reaper_locus.name, "additive":qtl.additive, "dominance":qtl.dominance}
- qtl_results.append(qtl)
-
-
- return qtl_results, json_data, perm_output, suggestive, significant, bootstrap_results
+def gen_reaper_results(this_trait, dataset, samples_before, json_data, num_perm, bootCheck, num_bootstrap, do_control, control_marker, manhattan_plot): + genotype = dataset.group.read_genotype_file() + + if manhattan_plot != True: + genotype = genotype.addinterval() + + samples, values, variances, sample_aliases = this_trait.export_informative() + + trimmed_samples = [] + trimmed_values = [] + for i in range(0, len(samples)): + if this_trait.data[samples[i]].name in samples_before: + trimmed_samples.append(samples[i]) + trimmed_values.append(values[i]) + + perm_output = [] + bootstrap_results = [] + + if num_perm < 100: + suggestive = 0 + significant = 0 + else: + perm_output = genotype.permutation(strains = trimmed_samples, trait = trimmed_values, nperm=num_perm) + suggestive = perm_output[int(num_perm*0.37-1)] + significant = perm_output[int(num_perm*0.95-1)] + highly_significant = perm_output[int(num_perm*0.99-1)] + + json_data['suggestive'] = suggestive + json_data['significant'] = significant + + if control_marker != "" and do_control == "true": + reaper_results = genotype.regression(strains = trimmed_samples, + trait = trimmed_values, + control = str(control_marker)) + if bootCheck: + control_geno = [] + control_geno2 = [] + _FIND = 0 + for _chr in genotype: + for _locus in _chr: + if _locus.name == control_marker: + control_geno2 = _locus.genotype + _FIND = 1 + break + if _FIND: + break + if control_geno2: + _prgy = list(genotype.prgy) + for _strain in trimmed_samples: + _idx = _prgy.index(_strain) + control_geno.append(control_geno2[_idx]) + + bootstrap_results = genotype.bootstrap(strains = trimmed_samples, + trait = trimmed_values, + control = control_geno, + nboot = num_bootstrap) + else: + reaper_results = genotype.regression(strains = trimmed_samples, + trait = trimmed_values) + + if bootCheck: + bootstrap_results = genotype.bootstrap(strains = trimmed_samples, + trait = trimmed_values, + nboot = num_bootstrap) + + json_data['chr'] = [] + json_data['pos'] = [] + json_data['lod.hk'] = [] + json_data['markernames'] = [] + #if self.additive: + # self.json_data['additive'] = [] + + #Need to convert the QTL objects that qtl reaper returns into a json serializable dictionary + qtl_results = [] + for qtl in reaper_results: + reaper_locus = qtl.locus + #ZS: Convert chr to int + converted_chr = reaper_locus.chr + if reaper_locus.chr != "X" and reaper_locus.chr != "X/Y": + converted_chr = int(reaper_locus.chr) + json_data['chr'].append(converted_chr) + json_data['pos'].append(reaper_locus.Mb) + json_data['lod.hk'].append(qtl.lrs) + json_data['markernames'].append(reaper_locus.name) + #if self.additive: + # self.json_data['additive'].append(qtl.additive) + locus = {"name":reaper_locus.name, "chr":reaper_locus.chr, "cM":reaper_locus.cM, "Mb":reaper_locus.Mb} + qtl = {"lrs_value": qtl.lrs, "chr":converted_chr, "Mb":reaper_locus.Mb, + "cM":reaper_locus.cM, "name":reaper_locus.name, "additive":qtl.additive, "dominance":qtl.dominance} + qtl_results.append(qtl) + + + return qtl_results, json_data, perm_output, suggestive, significant, bootstrap_results |