aboutsummaryrefslogtreecommitdiff
path: root/gn3/computations/rqtl.py
diff options
context:
space:
mode:
authorBonfaceKilz2022-03-23 15:52:25 +0300
committerBonfaceKilz2022-03-23 16:08:26 +0300
commitf8897da6e4069b592ca5846156e472b76a235724 (patch)
tree2acd96fbd1dfb62c0bf5fb949b5294d02ede58eb /gn3/computations/rqtl.py
parentf4b02281cdf0a29d82caf8c06ce96f947e5cf623 (diff)
downloadgenenetwork3-f8897da6e4069b592ca5846156e472b76a235724.tar.gz
Run python-black on file and fix other pylint issues.
See: <https://ci.genenetwork.org/jobs/genenetwork3-pylint/126> * gn3/computations/rqtl.py: Run `black gn3/computations/rqtl.py`. Also, manually fix other pylint issues.
Diffstat (limited to 'gn3/computations/rqtl.py')
-rw-r--r--gn3/computations/rqtl.py222
1 files changed, 146 insertions, 76 deletions
diff --git a/gn3/computations/rqtl.py b/gn3/computations/rqtl.py
index 65ee6de..2f9c4a2 100644
--- a/gn3/computations/rqtl.py
+++ b/gn3/computations/rqtl.py
@@ -11,24 +11,34 @@ 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
-def generate_rqtl_cmd(rqtl_wrapper_cmd: str,
- rqtl_wrapper_kwargs: Dict,
- rqtl_wrapper_bool_kwargs: list) -> Dict:
+
+def generate_rqtl_cmd(
+ rqtl_wrapper_cmd: str,
+ rqtl_wrapper_kwargs: Dict,
+ rqtl_wrapper_bool_kwargs: list,
+) -> Dict:
"""Given the base rqtl_wrapper command and
-dict of keyword arguments, return the full rqtl_wrapper command and an
-output filename generated from a hash of the genotype and phenotype files"""
+ dict of keyword arguments, return the full rqtl_wrapper command and an
+ output filename generated from a hash of the genotype and phenotype files"""
# Generate a hash from contents of the genotype and phenotype files
_hash = get_hash_of_files(
- [v for k, v in rqtl_wrapper_kwargs.items() if k in ["g", "p"]])
+ [v for k, v in rqtl_wrapper_kwargs.items() if k in ["g", "p"]]
+ )
# Append to hash a hash of keyword arguments
_hash += generate_hash_of_string(
- ",".join([f"{k}:{v}" for k, v in rqtl_wrapper_kwargs.items() if k not in ["g", "p"]]))
+ ",".join(
+ [
+ f"{k}:{v}"
+ for k, v in rqtl_wrapper_kwargs.items()
+ if k not in ["g", "p"]
+ ]
+ )
+ )
# Append to hash a hash of boolean keyword arguments
- _hash += generate_hash_of_string(
- ",".join(rqtl_wrapper_bool_kwargs))
+ _hash += generate_hash_of_string(",".join(rqtl_wrapper_bool_kwargs))
# Temporarily substitute forward-slashes in hash with underscores
_hash = _hash.replace("/", "_")
@@ -37,12 +47,12 @@ output filename generated from a hash of the genotype and phenotype files"""
rqtl_wrapper_kwargs["filename"] = _output_filename
return {
- "output_file":
- _output_filename,
- "rqtl_cmd":
- compose_rqtl_cmd(rqtl_wrapper_cmd=rqtl_wrapper_cmd,
- rqtl_wrapper_kwargs=rqtl_wrapper_kwargs,
- rqtl_wrapper_bool_kwargs=rqtl_wrapper_bool_kwargs)
+ "output_file": _output_filename,
+ "rqtl_cmd": compose_rqtl_cmd(
+ rqtl_wrapper_cmd=rqtl_wrapper_cmd,
+ rqtl_wrapper_kwargs=rqtl_wrapper_kwargs,
+ rqtl_wrapper_bool_kwargs=rqtl_wrapper_bool_kwargs,
+ ),
}
@@ -52,8 +62,13 @@ def process_rqtl_mapping(file_name: str) -> List:
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", "/tmp"),
- "output", file_name), "r", encoding="utf-8") as the_file:
+ with open(
+ os.path.join(
+ current_app.config.get("TMPDIR", "/tmp"), "output", file_name
+ ),
+ "r",
+ encoding="utf-8",
+ ) as the_file:
for line in the_file:
line_items = line.split(",")
if line_items[1][1:-1] == "chr" or not line_items:
@@ -71,53 +86,70 @@ def process_rqtl_mapping(file_name: str) -> List:
"chr": the_chr,
"cM": float(line_items[2]),
"Mb": float(line_items[2]),
- "lod_score": float(line_items[3])
+ "lod_score": float(line_items[3]),
}
marker_obs.append(this_marker)
return marker_obs
+
def process_rqtl_pairscan(file_name: str, geno_file: str) -> List:
"""Given an output file name, read in R/qtl pair-scan results and return
-a list of both the JSON needed for the d3panels figure and a list of results
-to be used when generating the results table (which will include marker names)"""
+ a list of both the JSON needed for the d3panels figure and a list of results
+ to be used when generating the results table (which will include marker names)"""
figure_data = pairscan_for_figure(file_name)
table_data = pairscan_for_table(file_name, geno_file)
return [figure_data, table_data]
+
def pairscan_for_figure(file_name: str) -> Dict:
"""Given an output file name, read in R/qtl pair-scan results and return
the JSON needed for the d3panels figure"""
figure_data = {}
# Open the file with the actual results, written as a list of lists
- with open(os.path.join(current_app.config.get("TMPDIR", "/tmp"),
- "output", file_name), "r") as the_file:
+ with open(
+ os.path.join(
+ current_app.config.get("TMPDIR", "/tmp"), "output", file_name
+ ),
+ "r",
+ encoding="utf8",
+ ) as the_file:
lod_results = []
for i, line in enumerate(the_file):
- if i == 0: # Skip first line
+ if i == 0: # Skip first line
continue
- line_items = [item.rstrip('\n') for item in line.split(",")]
- lod_results.append(line_items[1:]) # Append all but first item in line
- figure_data['lod'] = lod_results
-
- # Open the map file with the list of markers/pseudomarkers and their positions
- with open(os.path.join(current_app.config.get("TMPDIR", "/tmp"),
- "output", "MAP_" + file_name), "r") as the_file:
- chr_list = [] # type: List
- pos_list = [] # type: List
+ line_items = [item.rstrip("\n") for item in line.split(",")]
+ # Append all but first item in line
+ lod_results.append(line_items[1:])
+ figure_data["lod"] = lod_results
+
+ # Open the map file with the list of markers/pseudomarkers and their
+ # positions
+ with open(
+ os.path.join(
+ current_app.config.get("TMPDIR", "/tmp"),
+ "output",
+ "MAP_" + file_name,
+ ),
+ "r",
+ encoding="utf8",
+ ) as the_file:
+ chr_list = [] # type: List
+ pos_list = [] # type: List
for i, line in enumerate(the_file):
- if i == 0: # Skip first line
+ if i == 0: # Skip first line
continue
- line_items = [item.rstrip('\n') for item in line.split(",")]
+ line_items = [item.rstrip("\n") for item in line.split(",")]
chr_list.append(line_items[1][1:-1])
pos_list.append(line_items[2])
- figure_data['chr'] = chr_list
- figure_data['pos'] = pos_list
+ figure_data["chr"] = chr_list
+ figure_data["pos"] = pos_list
return figure_data
+
def get_marker_list(map_file: str) -> List:
"""
Open the map file with the list of markers/pseudomarkers and create list of marker obs
@@ -127,21 +159,31 @@ def get_marker_list(map_file: str) -> List:
"""
marker_list = []
- with open(os.path.join(current_app.config.get("TMPDIR", "/tmp"),
- "output", map_file), "r") as map_fh:
+ with open(
+ os.path.join(
+ current_app.config.get("TMPDIR", "/tmp"), "output", map_file
+ ),
+ "r",
+ encoding="utf8",
+ ) as map_fh:
for line in map_fh.readlines()[1:]:
- line_items = [item.rstrip('\n') for item in line.split(",")]
+ line_items = [item.rstrip("\n") for item in line.split(",")]
this_marker = {
- 'name': line_items[0],
- 'chr': line_items[1][1:-1], # Strip quotes from beginning and end of chr string
- 'pos': line_items[2]
+ "name": line_items[0],
+ "chr": line_items[1][
+ 1:-1
+ ], # Strip quotes from beginning and end of chr string
+ "pos": line_items[2],
}
marker_list.append(this_marker)
return marker_list
-def generate_table_rows(results_file: str, marker_list: List, original_markers: Dict) -> List:
+
+def generate_table_rows(
+ results_file: str, marker_list: List, original_markers: Dict
+) -> List:
"""
Open the file with the actual R/qtl pair-scan results and write them as
they will be displayed in the results table
@@ -152,38 +194,45 @@ def generate_table_rows(results_file: str, marker_list: List, original_markers:
original_markers: Dict of markers from the .geno file, for finding proximal/distal
markers to each pseudomarker
"""
-
table_data = []
- with open(os.path.join(current_app.config.get("TMPDIR", "/tmp"),
- "output", results_file), "r") as the_file:
+ with open(
+ os.path.join(
+ current_app.config.get("TMPDIR", "/tmp"), "output", results_file
+ ),
+ "r",
+ encoding="utf8",
+ ) as the_file:
for i, line in enumerate(the_file.readlines()[1:]):
marker_1 = marker_list[i]
- marker_1['proximal'], marker_1['distal'] = find_nearest_marker(marker_1['chr'],
- marker_1['pos'],
- original_markers)
- line_items = [item.rstrip('\n') for item in line.split(",")]
+ marker_1["proximal"], marker_1["distal"] = find_nearest_marker(
+ marker_1["chr"], marker_1["pos"], original_markers
+ )
+ line_items = [item.rstrip("\n") for item in line.split(",")]
for j, item in enumerate(line_items[1:]):
marker_2 = marker_list[j]
- marker_2['proximal'], marker_2['distal'] = find_nearest_marker(marker_2['chr'],
- marker_2['pos'],
- original_markers)
+ marker_2["proximal"], marker_2["distal"] = find_nearest_marker(
+ marker_2["chr"], marker_2["pos"], original_markers
+ )
try:
lod_score = f"{float(item):.3f}"
except ValueError:
lod_score = f"{item}"
- table_data.append({
- 'proximal1': marker_1['proximal'],
- 'distal1': marker_1['distal'],
- 'pos1': f"Chr {marker_1['chr']} @ {float(marker_1['pos']):.1f} cM",
- 'lod': lod_score,
- 'proximal2': marker_2['proximal'],
- 'distal2': marker_2['distal'],
- 'pos2': f"Chr {marker_2['chr']} @ {float(marker_2['pos']):.1f} cM"
- })
+ table_data.append(
+ {
+ "proximal1": marker_1["proximal"],
+ "distal1": marker_1["distal"],
+ "pos1": f"Chr {marker_1['chr']} @ {float(marker_1['pos']):.1f} cM",
+ "lod": lod_score,
+ "proximal2": marker_2["proximal"],
+ "distal2": marker_2["distal"],
+ "pos2": f"Chr {marker_2['chr']} @ {float(marker_2['pos']):.1f} cM",
+ }
+ )
return table_data
+
def pairscan_for_table(file_name: str, geno_file: str) -> List:
"""
Read in R/qtl pair-scan results and return as List representing
@@ -204,37 +253,50 @@ def pairscan_for_table(file_name: str, geno_file: str) -> List:
# they will be displayed in the results table
table_data = generate_table_rows(file_name, marker_list, original_markers)
- return sorted(table_data, key=lambda i: float(i['lod']), reverse=True)[:500]
+ return sorted(table_data, key=lambda i: float(i["lod"]), reverse=True)[:500]
+
def build_marker_pos_dict(genotype_file: str) -> Dict:
"""Gets list of markers and their positions from .geno file
Basically a pared-down version of parse_genotype_file for R/qtl pair-scan"""
- with open(genotype_file, "r") as infile:
+ with open(genotype_file, "r", encoding="utf8") as infile:
contents = infile.readlines()
# Get all lines after the metadata
- lines = tuple(line for line in contents if
- ((not line.strip().startswith("#")) and
- (not line.strip().startswith("@")) and
- (not line.strip() == "")))
+ lines = tuple(
+ line
+ for line in contents
+ if (
+ (not line.strip().startswith("#"))
+ and (not line.strip().startswith("@"))
+ and (not line.strip() == "")
+ )
+ )
header_items = lines[0].split("\t")
mb_exists = "Mb" in header_items
- pos_column = header_items.index("Mb") if mb_exists else header_items.index("cM")
+ pos_column = (
+ header_items.index("Mb") if mb_exists else header_items.index("cM")
+ )
- the_markers = {"1": {}} # type: Dict[str, Dict]
- for line in lines[1:]: # The lines with markers
+ the_markers = {"1": {}} # type: Dict[str, Dict]
+ for line in lines[1:]: # The lines with markers
line_items = line.split("\t")
this_chr = line_items[0]
if this_chr not in the_markers:
the_markers[this_chr] = {}
- the_markers[this_chr][str(float(line_items[pos_column]))] = line_items[1]
+ the_markers[this_chr][str(float(line_items[pos_column]))] = line_items[
+ 1
+ ]
return the_markers
-def find_nearest_marker(the_chr: str, the_pos: str, marker_list: Dict) -> Tuple[str, str]:
+
+def find_nearest_marker(
+ the_chr: str, the_pos: str, marker_list: Dict
+) -> Tuple[str, str]:
"""Given a chromosome and position of a pseudomarker (from R/qtl pair-scan results),
return the nearest real marker"""
@@ -243,18 +305,26 @@ def find_nearest_marker(the_chr: str, the_pos: str, marker_list: Dict) -> Tuple[
# Get the position of the pseudomarker in the list of markers for the chr
the_pos_index = bisect(pos_list, float(the_pos))
- proximal_marker = marker_list[the_chr][str(pos_list[the_pos_index-1])]
+ proximal_marker = marker_list[the_chr][str(pos_list[the_pos_index - 1])]
distal_marker = marker_list[the_chr][str(pos_list[the_pos_index])]
return proximal_marker, distal_marker
+
def process_perm_output(file_name: str) -> Tuple[List, float, float]:
"""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", "/tmp"),
- "output", "PERM_" + file_name), "r", encoding="utf-8") as the_file:
+ with open(
+ os.path.join(
+ current_app.config.get("TMPDIR", "/tmp"),
+ "output",
+ "PERM_" + file_name,
+ ),
+ "r",
+ encoding="utf-8",
+ ) as the_file:
for i, line in enumerate(the_file):
if i == 0:
# Skip header line