1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
|
"""Run quality control on phenotypes-specific files in the bundle."""
import sys
import shutil
import tempfile
from pathlib import Path
from zipfile import ZipFile
from functools import reduce
import multiprocessing as mproc
from logging import Logger, getLogger, StreamHandler
import MySQLdb as mdb
from r_qtl import r_qtl2 as rqtl2
from r_qtl import r_qtl2_qc as rqc
from r_qtl import exceptions as rqe
from uploader.files import sha256_digest_over_file
from uploader.samples.models import samples_by_species_and_population
from scripts.rqtl2.entry import build_main
from scripts.rqtl2.cli_parser import add_bundle_argument
from scripts.cli_parser import init_cli_parser, add_global_data_arguments
from scripts.rqtl2.bundleutils import build_line_joiner, build_line_splitter
def validate(phenobundle: Path, logger: Logger) -> dict:
"""Check that the bundle is generally valid"""
try:
rqc.validate_bundle(phenobundle)
except rqe.RQTLError as rqtlerr:
# logger.error("Bundle file validation failed!", exc_info=True)
return {
"skip": True,
"logger": logger,
"phenobundle": phenobundle,
"errors": (" ".join(rqtlerr.args),)
}
return {
"errors": tuple(),
"skip": False,
"phenobundle": phenobundle,
"logger": logger
}
def check_for_mandatory_pheno_keys(
phenobundle: Path,
logger: Logger,
**kwargs
) -> dict:
"""Check that the mandatory keys exist for phenotypes."""
if kwargs.get("skip", False):
return {
**kwargs,
"logger": logger,
"phenobundle": phenobundle
}
_mandatory_keys = ("pheno", "phenocovar")
_cdata = rqtl2.read_control_file(phenobundle)
_errors = kwargs.get("errors", tuple()) + tuple(
f"Expected '{key}' file(s) are not declared in the bundle."
for key in _mandatory_keys if key not in _cdata.keys())
return {
**kwargs,
"logger": logger,
"phenobundle": phenobundle,
"errors": _errors,
"skip": len(_errors) > 0
}
def check_for_averages_files(
phenobundle: Path,
logger: Logger,
**kwargs
) -> dict:
"""Check that averages files appear together"""
if kwargs.get("skip", False):
return {
**kwargs,
"logger": logger,
"phenobundle": phenobundle
}
_together = (("phenose", "phenonum"), ("phenonum", "phenose"))
_cdata = rqtl2.read_control_file(phenobundle)
_errors = kwargs.get("errors", tuple()) + tuple(
f"'{first}' is defined in the control file but there is no "
f"corresponding '{second}'"
for first, second in _together
if ((first in _cdata.keys()) and (second not in _cdata.keys())))
return {
**kwargs,
"logger": logger,
"phenobundle": phenobundle,
"errors": _errors,
"skip": len(_errors) > 0
}
def extract_bundle(bundle: Path, workdir: Path) -> tuple[Path, tuple[Path, ...]]:
"""Extract the bundle."""
with ZipFile(bundle) as zfile:
extractiondir = workdir.joinpath(
f"{sha256_digest_over_file(bundle)}_phenotype_qc_{bundle.name}")
return extractiondir, rqtl2.extract(zfile, extractiondir)
def undo_transpose(filetype: str, cdata: dict, extractiondir):
"""Undo transposition of all files of type `filetype` in thebundle."""
if len(cdata.get(filetype, [])) > 0 and cdata.get(f"{filetype}_transposed", False):
files = (extractiondir.joinpath(_file) for _file in cdata[filetype])
for _file in files:
rqtl2.transpose_csv_with_rename(
_file,
build_line_splitter(cdata),
build_line_joiner(cdata))
def run_qc(# pylint: disable=[too-many-arguments]
dbconn: mdb.Connection,
phenobundle: Path,
workingdir: Path,
speciesid: int,
populationid: int,
logger: Logger
) -> int:
"""Run quality control checks on the bundle."""
results = check_for_averages_files(
**check_for_mandatory_pheno_keys(
**validate(phenobundle, logger)))
errors = results.get("errors", tuple())
if len(errors) > 0:
logger.error("We found the following errors:\n%s",
"\n".join(f" - {error}" for error in errors))
return 1
# TODO: Run QC on actual values
# Steps:
# - Extract file to specific directory
extractiondir, *_bundlefiles = extract_bundle(phenobundle, workingdir)
# - For every pheno, phenocovar, phenose, phenonum file, undo
# transposition where relevant
cdata = rqtl2.control_data(extractiondir)
with mproc.Pool(mproc.cpu_count() - 1) as pool:
pool.starmap(
undo_transpose,
((ftype, cdata, extractiondir)
for ftype in ("pheno", "phenocovar", "phenose", "phenonum")))
# - Fetch samples/individuals from database.
samples = tuple(
item for item in set(reduce(
lambda acc, item: acc + (
item["Name"], item["Name2"], item["Symbol"], item["Alias"]),
samples_by_species_and_population(dbconn, speciesid, populationid),
tuple()))
if bool(item))
# - Check that `description` and `units` is present in phenocovar for
# all phenotypes
# qc_phenocovar_files(…)
# - Check all samples in pheno files exist in database
# - Check all phenotypes in pheno files exist in phenocovar files
# - Check all numeric values in pheno files
# qc_pheno_files(…)
# - Check the 3 checks above for phenose and phenonum values too
# qc_phenose_files(…)
# qc_phenonum_files(…)
# - Delete all extracted files
shutil.rmtree(extractiondir)
raise NotImplementedError("WIP!")
if __name__ == "__main__":
def cli_args():
"""Process command-line arguments for `install_phenos`"""
parser = add_bundle_argument(add_global_data_arguments(init_cli_parser(
program="PhenotypesQC",
description=(
"Perform Quality Control checks on a phenotypes bundle file"))))
parser.add_argument(
"--working-dir",
default=f"{tempfile.gettempdir()}/phenotypes_qc",
help=("The directory where this script will put its intermediate "
"files."),
type=Path)
return parser.parse_args()
_logger = getLogger("phenotypes_qc")
_logger.addHandler(StreamHandler(stream=sys.stderr))
main = build_main(
cli_args(),
lambda dbconn, args: run_qc(dbconn,
args.rqtl2bundle,
args.working_dir,
args.speciesid,
args.populationid,
_logger),
_logger,
"DEBUG")
sys.exit(main())
|
