diff options
Diffstat (limited to 'gn3')
| -rw-r--r-- | gn3/api/rqtl.py | 16 | ||||
| -rw-r--r-- | gn3/computations/rqtl.py | 64 |
2 files changed, 75 insertions, 5 deletions
diff --git a/gn3/api/rqtl.py b/gn3/api/rqtl.py index de620f7..0194b6f 100644 --- a/gn3/api/rqtl.py +++ b/gn3/api/rqtl.py @@ -1,10 +1,13 @@ """Endpoints for running the rqtl cmd""" +import os + from flask import Blueprint from flask import current_app from flask import jsonify from flask import request -from gn3.computations.rqtl import generate_rqtl_cmd +from gn3.commands import run_cmd +from gn3.computations.rqtl import generate_rqtl_cmd, process_rqtl_output, process_perm_output from gn3.computations.gemma import do_paths_exist rqtl = Blueprint("rqtl", __name__) @@ -35,10 +38,17 @@ run the rqtl_wrapper script and return the results as JSON if kwarg in boolean_kwargs: rqtl_bool_kwargs.append(kwarg) - results = generate_rqtl_cmd( + rqtl_cmd = generate_rqtl_cmd( rqtl_wrapper_cmd=current_app.config.get("RQTL_WRAPPER_CMD"), rqtl_wrapper_kwargs=rqtl_kwargs, rqtl_wrapper_bool_kwargs=boolean_kwargs ) - return jsonify(results) + os.system(rqtl_cmd.get('rqtl_cmd')) + + rqtl_output = {} + rqtl_output['results'] = process_rqtl_output(rqtl_cmd.get('output_file')) + if int(rqtl_kwargs['nperm']) > 0: + rqtl_output['perm_results'], rqtl_output['suggestive'], rqtl_output['significant'] = process_perm_output(rqtl_cmd.get('output_file')) + + return jsonify(rqtl_output) diff --git a/gn3/computations/rqtl.py b/gn3/computations/rqtl.py index 605e0e1..22d9faf 100644 --- a/gn3/computations/rqtl.py +++ b/gn3/computations/rqtl.py @@ -1,6 +1,11 @@ """Procedures related rqtl computations""" - +import os +import numpy as np from typing import Dict +from typing import List + +from flask import current_app + from gn3.commands import compose_rqtl_cmd from gn3.computations.gemma import generate_hash_of_string from gn3.fs_helpers import get_hash_of_files @@ -29,7 +34,9 @@ output filename generated from a hash of the genotype and phenotype files # Temporarily substitute forward-slashes in hash with underscores _hash = _hash.replace("/", "_") - _output_filename = f"{_hash}-output.json" + _output_filename = f"{_hash}-output.csv" + rqtl_wrapper_kwargs["filename"] = _output_filename + return { "output_file": _output_filename, @@ -38,3 +45,56 @@ output filename generated from a hash of the genotype and phenotype files rqtl_wrapper_kwargs=rqtl_wrapper_kwargs, rqtl_wrapper_bool_kwargs=rqtl_wrapper_bool_kwargs) } + + +def process_rqtl_output(file_name: str) -> List: + """Given an output file name, read in R/qtl results and return + a List of marker objects + + """ + marker_obs = [] + # Later I should probably redo this using csv.read to avoid the + # awkwardness with removing quotes with [1:-1] + with open(os.path.join(current_app.config.get("TMPDIR"), "output", file_name), "r") as the_file: + 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 + try: + the_chr = int(line_items[1][1:-1]) + except: + the_chr = line_items[1][1:-1] + this_marker = { + "name": line_items[0][1:-1], + "chr": the_chr, + "cM": float(line_items[2]), + "Mb": float(line_items[2]), + "lod_score": float(line_items[3]) + } + marker_obs.append(this_marker) + + return marker_obs + + +def process_perm_output(file_name: str): + """Given base filename, read in R/qtl permutation output and calculate + suggestive and significant thresholds + + """ + perm_results = [] + with open(os.path.join(current_app.config.get("TMPDIR"), "output", "PERM_" + file_name), "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])) + + suggestive = np.percentile(np.array(perm_results), 67) + significant = np.percentile(np.array(perm_results), 95) + + return perm_results, suggestive, significant |