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 |