aboutsummaryrefslogtreecommitdiff
path: root/wqflask
diff options
context:
space:
mode:
Diffstat (limited to 'wqflask')
-rw-r--r--wqflask/wqflask/marker_regression/qtlreaper_mapping.py42
-rw-r--r--wqflask/wqflask/marker_regression/run_mapping.py39
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