diff options
Diffstat (limited to 'wqflask')
| -rw-r--r-- | wqflask/wqflask/marker_regression/qtlreaper_mapping.py | 42 | ||||
| -rw-r--r-- | wqflask/wqflask/marker_regression/run_mapping.py | 39 |
2 files changed, 60 insertions, 21 deletions
diff --git a/wqflask/wqflask/marker_regression/qtlreaper_mapping.py b/wqflask/wqflask/marker_regression/qtlreaper_mapping.py index 0c560582..6b4c05ea 100644 --- a/wqflask/wqflask/marker_regression/qtlreaper_mapping.py +++ b/wqflask/wqflask/marker_regression/qtlreaper_mapping.py @@ -37,6 +37,8 @@ def run_reaper(this_trait, this_dataset, samples, vals, json_data, num_perm, boo opt_list.append("--permu_output " + webqtlConfig.GENERATED_IMAGE_DIR + permu_filename + ".txt") if control_marker != "" and do_control == "true": opt_list.append("-c " + control_marker) + if manhattan_plot != True: + opt_list.append("--interval 1") reaper_command = REAPER_COMMAND + ' --geno {0}/{1}.geno --traits {2}/gn2/{3}.txt {4} -o {5}{6}.txt'.format(flat_files('genotype'), genofile_name, @@ -85,9 +87,17 @@ def parse_reaper_output(gwa_filename, permu_filename, bootstrap_filename): p_values = [] marker_obs = [] + only_cm = False + only_mb = False + with open("{}{}.txt".format(webqtlConfig.GENERATED_IMAGE_DIR, gwa_filename)) as output_file: for line in output_file: if line.startswith("ID\t"): + if len(line.split("\t")) < 8: + if 'cM' in line.split("\t"): + only_cm = True + else: + only_mb = True continue else: marker = {} @@ -96,16 +106,30 @@ def parse_reaper_output(gwa_filename, permu_filename, bootstrap_filename): marker['chr'] = int(line.split("\t")[2]) except: marker['chr'] = line.split("\t")[2] - marker['cM'] = float(line.split("\t")[3]) - if float(line.split("\t")[4]) > 1000: - marker['Mb'] = float(line.split("\t")[4])/1000000 + if only_cm or only_mb: + if only_cm: + marker['cM'] = float(line.split("\t")[3]) + else: + if float(line.split("\t")[3]) > 1000: + marker['Mb'] = float(line.split("\t")[3])/1000000 + else: + marker['Mb'] = float(line.split("\t")[3]) + if float(line.split("\t")[6]) != 1: + marker['p_value'] = float(line.split("\t")[6]) + marker['lrs_value'] = float(line.split("\t")[4]) + marker['lod_score'] = marker['lrs_value'] / 4.61 + marker['additive'] = float(line.split("\t")[5]) else: - marker['Mb'] = float(line.split("\t")[4]) - if float(line.split("\t")[7]) != 1: - marker['p_value'] = float(line.split("\t")[7]) - marker['lrs_value'] = float(line.split("\t")[5]) - marker['lod_score'] = marker['lrs_value'] / 4.61 - marker['additive'] = float(line.split("\t")[6]) + marker['cM'] = float(line.split("\t")[3]) + if float(line.split("\t")[4]) > 1000: + marker['Mb'] = float(line.split("\t")[4])/1000000 + else: + marker['Mb'] = float(line.split("\t")[4]) + if float(line.split("\t")[7]) != 1: + marker['p_value'] = float(line.split("\t")[7]) + marker['lrs_value'] = float(line.split("\t")[5]) + marker['lod_score'] = marker['lrs_value'] / 4.61 + marker['additive'] = float(line.split("\t")[6]) marker_obs.append(marker) #ZS: Results have to be reordered because the new reaper returns results sorted alphabetically by chr for some reason, resulting in chr 1 being followed by 10, etc diff --git a/wqflask/wqflask/marker_regression/run_mapping.py b/wqflask/wqflask/marker_regression/run_mapping.py index c9d10f7c..9bde343c 100644 --- a/wqflask/wqflask/marker_regression/run_mapping.py +++ b/wqflask/wqflask/marker_regression/run_mapping.py @@ -374,10 +374,15 @@ class RunMapping(object): self.annotations_for_browser = [] highest_chr = 1 #This is needed in order to convert the highest chr to X/Y for marker in results: + if 'Mb' in marker: + this_ps = marker['Mb']*1000000 + else: + this_ps = marker['cM']*1000000 + browser_marker = dict( chr = str(marker['chr']), rs = marker['name'], - ps = marker['Mb']*1000000, + ps = this_ps, url = "/show_trait?trait_id=" + marker['name'] + "&dataset=" + self.dataset.group.name + "Geno" ) @@ -386,7 +391,7 @@ class RunMapping(object): name = str(marker['name']), chr = str(marker['chr']), rs = marker['name'], - pos = marker['Mb']*1000000, + pos = this_ps, url = "/show_trait?trait_id=" + marker['name'] + "&dataset=" + self.dataset.group.name + "Geno" ) else: @@ -394,7 +399,7 @@ class RunMapping(object): name = str(marker['name']), chr = str(marker['chr']), rs = marker['name'], - pos = marker['Mb']*1000000 + pos = this_ps ) #if 'p_value' in marker: # logger.debug("P EXISTS:", marker['p_value']) @@ -533,10 +538,10 @@ def export_mapping_results(dataset, trait, markers, results_path, mapping_scale, output_file.write("\n") output_file.write("Name,Chr,") if score_type.lower() == "-log(p)": - score_type = "'-log(p)" - if mapping_scale == "physic": + score_type = "-log(p)" + if 'Mb' in markers[0]: output_file.write("Mb," + score_type) - else: + if 'cM' in markers[0]: output_file.write("Cm," + score_type) if "additive" in markers[0].keys(): output_file.write(",Additive") @@ -544,7 +549,11 @@ def export_mapping_results(dataset, trait, markers, results_path, mapping_scale, output_file.write(",Dominance") output_file.write("\n") for i, marker in enumerate(markers): - output_file.write(marker['name'] + "," + str(marker['chr']) + "," + str(marker['Mb']) + ",") + output_file.write(marker['name'] + "," + str(marker['chr']) + ",") + if 'Mb' in marker: + output_file.write(str(marker['Mb']) + ",") + if 'cM' in marker: + output_file.write(str(marker['cM']) + ",") if "lod_score" in marker.keys(): output_file.write(str(marker['lod_score'])) else: @@ -669,16 +678,22 @@ def get_chr_lengths(mapping_scale, mapping_method, dataset, qtl_results): except: chr_as_num = 20 if chr_as_num > this_chr or i == (len(qtl_results) - 1): - chr_lengths.append({ "chr": str(this_chr), "size": str(highest_pos)}) - this_chr = chr_as_num - highest_pos = 0 + if i == (len(qtl_results) - 1): + if mapping_method == "reaper": + highest_pos = float(result['cM']) * 1000000 + else: + highest_pos = float(result['Mb']) * 1000000 + chr_lengths.append({ "chr": str(this_chr), "size": str(highest_pos)}) + else: + chr_lengths.append({ "chr": str(this_chr), "size": str(highest_pos)}) + this_chr = chr_as_num else: if mapping_method == "reaper": if float(result['cM']) > highest_pos: - highest_pos = float(result['cM']) + highest_pos = float(result['cM']) * 1000000 else: if float(result['Mb']) > highest_pos: - highest_pos = float(result['Mb']) + highest_pos = float(result['Mb']) * 1000000 return chr_lengths |