run_rqtl now reads in permutation output

1 files changed, 32 insertions, 2 deletions

diff --git a/wqflask/wqflask/marker_regression/rqtl_mapping.py b/wqflask/wqflask/marker_regression/rqtl_mapping.py
index 3f4899b0..c3f047fc 100644
--- a/wqflask/wqflask/marker_regression/rqtl_mapping.py
+++ b/wqflask/wqflask/marker_regression/rqtl_mapping.py
@@ -8,6 +8,8 @@ from typing import List
 from typing import Optional
 from typing import TextIO
 
+import numpy as np
+
 from base.webqtlConfig import TMPDIR
 from base.trait import create_trait
 from utility.tools import locate
@@ -46,7 +48,11 @@ def run_rqtl(trait_name, vals, samples, dataset, mapping_scale, method, model, p
 
     out_file = requests.post(GN3_RQTL_URL, data=post_data).json()['output_file']
 
-    return process_rqtl_results(out_file)
+    if num_perm > 0:
+        perm_results, suggestive, significant = process_perm_results(out_file)
+        return perm_results, suggestive, significant, process_rqtl_results(out_file)
+    else:
+        return process_rqtl_results(out_file)
 
 
 def process_rqtl_results(out_file: str) -> List:
@@ -63,6 +69,7 @@ def process_rqtl_results(out_file: str) -> List:
         for line in the_file:
             line_items = line.split(",")
             if line_items[1][1:-1] == "chr" or not line_items:
+                # Skip header line
                 continue
             else:
                 # Convert chr to int if possible
@@ -81,6 +88,30 @@ def process_rqtl_results(out_file: str) -> List:
 
     return marker_obs
 
+
+def process_perm_results(out_file: str):
+    """Given base filename, read in R/qtl permutation output and calculate
+    suggestive and significant thresholds
+
+    """
+    perm_results = []
+    with open(GN3_TMP_PATH + "/output/PERM_" + out_file, "r") as the_file:
+        for i, line in enumerate(the_file):
+            if i == 0:
+                # Skip header line
+                continue
+            else:
+                line_items = line.split(",")
+                perm_results.append(float(line_items[1]))
+
+    logger.debug("PERM RESULTS:", perm_results)
+
+    suggestive = np.percentile(np.array(perm_results), 67)
+    significant = np.percentile(np.array(perm_results), 95)
+
+    return perm_results, suggestive, significant
+
+
 def get_hash_of_textio(the_file: TextIO) -> str:
     """Given a StringIO, return the hash of its contents"""
 
@@ -146,7 +177,6 @@ def cofactors_to_dict(cofactors: str, dataset_ob, samples) -> Dict:
                                         name=cofactor_name)
                 sample_data = trait_ob.data
                 for index, sample in enumerate(samples):
-                    #if (sample in sample_list) and (sample in sample_data):
                     if sample in sample_data:
                         sample_value = sample_data[sample].value
                         cofactors[cofactor_name].append(sample_value)