"""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 r_qtl.fileerrors import InvalidValue 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 qc_phenocovar_file(phenocovarfile: Path, separator: str, comment_char: str): """Check that `phenocovar` files are structured correctly.""" _csvfile = rqtl2.read_csv_file(phenocovarfile, separator, comment_char) _headings = tuple(heading.lower() for heading in next(_csvfile)) _errors = tuple() for heading in ("description", "units"): if heading not in _headings: _errors = (InvalidValue( phenocovarfile.name, "header row", "-", "-", (f"File {phenocovarfile.name} is missing the {heading} heading " "in the header line.")),) def collect_errors(errors_and_linecount, line): _errs, _lc = errors_and_linecount if len(line) != len(_headings): _errs = _errs + (InvalidValue( phenocovarfile.name, line[0], "-", "-", (f"Record {_lc} in file {phenocovarfile.name} has a different " "number of columns than the number of headings")),) _line = dict(zip(_headings, line)) if not bool(_line["description"]): _errs = _errs + ( InvalidValue(phenocovarfile.name, _line[_headings[0]], "description", _line["description"], "The description is not provided!"),) return _errs, _lc+1 return { phenocovarfile.name: dict(zip( ("errors", "linecount"), reduce(collect_errors, _csvfile, (_errors, 1)))) } def merge_dicts(*dicts): """Merge multiple dicts into a single one.""" return reduce(lambda merged, dct: {**merged, **dct}, dicts, {}) 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 with mproc.Pool(mproc.cpu_count() - 1) as pool: # This call is way too busy. Maybe just return the errors? qc_results = merge_dicts(*pool.starmap(qc_phenocovar_file, tuple( (extractiondir.joinpath(_file), cdata["sep"], cdata["comment.char"]) for _file in cdata.get("phenocovar", [])))) # - 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())