This should fix an issue where the new qtlreaper wouldn't work properlty

with genotypes that only have cM positions
* wqflask/wqflask/marker_regression/qtlreaper_mapping.py -
parse_reaper_output is changed to check if cM and Mb exist in output
when creating marker obs
* wqflask/wqflask/marker_regression/run_mapping.py - Changed
export_mapping_results to properly deal with all combinations of cM and
Mb positions (both and either/or)

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