aboutsummaryrefslogtreecommitdiff
path: root/scripts/qc_on_rqtl2_bundle2.py
blob: 913624300b7670bdd01a4b429e87529be64caf20 (plain)
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
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
"""Run Quality Control checks on R/qtl2 bundle."""
import os
import sys
import json
from time import sleep
from pathlib import Path
from zipfile import ZipFile
from argparse import Namespace
from datetime import timedelta
import multiprocessing as mproc
from functools import reduce, partial
from logging import Logger, getLogger, StreamHandler
from typing import Union, Sequence, Callable, Iterator

import MySQLdb as mdb
from redis import Redis

from quality_control.errors import InvalidValue
from quality_control.checks import decimal_points_error

from uploader import jobs
from uploader.db_utils import database_connection
from uploader.check_connections import check_db, check_redis

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 import fileerrors as rqfe

from scripts.process_rqtl2_bundle import parse_job
from scripts.redis_logger import setup_redis_logger
from scripts.cli_parser import init_cli_parser, add_global_data_arguments


def build_line_splitter(cdata: dict) -> Callable[[str], tuple[Union[str, None], ...]]:
    """Build and return a function to use to split data in the files.

    Parameters
    ----------
    cdata: A dict holding the control information included with the R/qtl2
        bundle.

    Returns
    -------
    A function that takes a string and return a tuple of strings.
    """
    separator = cdata["sep"]
    na_strings = cdata["na.strings"]
    def __splitter__(line: str) -> tuple[Union[str, None], ...]:
        return tuple(
            item if item not in na_strings else None
            for item in
            (field.strip() for field in line.strip().split(separator)))
    return __splitter__


def build_line_joiner(cdata: dict) -> Callable[[tuple[Union[str, None], ...]], str]:
    """Build and return a function to use to split data in the files.

    Parameters
    ----------
    cdata: A dict holding the control information included with the R/qtl2
        bundle.

    Returns
    -------
    A function that takes a string and return a tuple of strings.
    """
    separator = cdata["sep"]
    na_strings = cdata["na.strings"]
    def __joiner__(row: tuple[Union[str, None], ...]) -> str:
        return separator.join(
            (na_strings[0] if item is None else item)
            for item in row)
    return __joiner__


def check_for_missing_files(
        rconn: Redis, fqjobid: str, extractpath: Path, logger: Logger) -> bool:
    """Check that all files listed in the control file do actually exist."""
    logger.info("Checking for missing files.")
    missing = rqc.missing_files(extractpath)
    # add_to_errors(rconn, fqjobid, "errors-generic", tuple(
    #     rqfe.MissingFile(
    #         mfile[0], mfile[1], (
    #             f"File '{mfile[1]}' is listed in the control file under "
    #             f"the '{mfile[0]}' key, but it does not actually exist in "
    #             "the bundle."))
    #     for mfile in missing))
    if len(missing) > 0:
        logger.error(f"Missing files in the bundle!")
        return True
    return False


def open_file(file_: Path) -> Iterator:
    """Open file and return one line at a time."""
    with open(file_, "r", encoding="utf8") as infile:
        for line in infile:
            yield line


def check_markers(
        filename: str,
        row: tuple[str, ...],
        save_error: lambda val: val
) -> tuple[rqfe.InvalidValue]:
    """Check that the markers are okay"""
    errors = tuple()
    counts = {}
    for marker in row:
        counts = {**counts, marker: counts.get(marker, 0) + 1}
        if marker is None or marker == "":
            errors = errors + (save_error(rqfe.InvalidValue(
                filename,
                "markers"
                "-",
                marker,
                "A marker MUST be a valid value.")),)

    return errors + tuple(
        save_error(rqfe.InvalidValue(
            filename,
            "markers",
            key,
            f"Marker '{key}' was repeated {value} times"))
        for key,value in counts.items() if value > 1)


def check_geno_line(
        filename: str,
        headers: tuple[str, ...],
        row: tuple[Union[str, None]],
        cdata: dict,
        save_error: lambda val: val
) -> tuple[rqfe.InvalidValue]:
    """Check that the geno line is correct."""
    errors = tuple()
    # Verify that line has same number of columns as headers
    if len(headers) != len(row):
        errors = errors + (save_error(rqfe.InvalidValue(
            filename,
            headers[0],
            row[0],
            row[0],
            "Every line MUST have the same number of columns.")),)

    # First column is the individuals/cases/samples
    if not bool(row[0]):
        errors = errors + (save_error(rqfe.InvalidValue(
            filename,
            headers[0],
            row[0],
            row[0],
            "The sample/case MUST be a valid value.")),)

    def __process_value__(val):
        if val in cdata["na.strings"]:
            return None
        if val in cdata["alleles"]:
            return cdata["genotypes"][val]

    genocode = cdata.get("genotypes", {})
    for coltitle, cellvalue in zip(headers[1:],row[1:]):
        if (
                bool(genocode) and
                cellvalue is not None and
                cellvalue not in genocode.keys()
        ):
            errors = errors + (save_error(rqfe.InvalidValue(
                filename, row[0], coltitle, cellvalue,
                f"Value '{cellvalue}' is invalid. Expected one of "
                f"'{', '.join(genocode.keys())}'.")),)

    return errors


def push_file_error_to_redis(rconn: Redis, key: str, error: InvalidValue) -> InvalidValue:
    """Push the file error to redis a json string

    Parameters
    ----------
    rconn: Connection to redis
    key: The name of the list where we push the errors
    error: The file error to save

    Returns
    -------
    Returns the file error it saved
    """
    if bool(error):
        rconn.rpush(key, json.dumps(error._asdict()))
    return error


def file_errors_and_details(
        redisargs: dict[str, str],
        file_: Path,
        filetype: str,
        cdata: dict,
        linesplitterfn: Callable,
        linejoinerfn: Callable,
        headercheckers: tuple[Callable, ...],
        bodycheckers: tuple[Callable, ...]
) -> dict:
    """Compute errors, and other file metadata."""
    errors = tuple()
    if cdata[f"{filetype}_transposed"]:
        rqtl2.transpose_csv_with_rename(file_, linesplitterfn, linejoinerfn)

    with Redis.from_url(redisargs["redisuri"], decode_responses=True) as rconn:
        save_error_fn = partial(push_file_error_to_redis,
                                rconn,
                                error_list_name(filetype, file_.name))
        for lineno, line in enumerate(open_file(file_), start=1):
            row = linesplitterfn(line)
            if lineno == 1:
                headers = tuple(row)
                errors = errors + reduce(
                    lambda errs, fnct: errs + fnct(
                        file_.name, row[1:], save_error_fn),
                    headercheckers,
                    tuple())
                continue

            errors = errors + reduce(
                lambda errs, fnct: errs + fnct(
                    file_.name, headers, row, cdata, save_error_fn),
                bodycheckers,
                tuple())

        filedetails = {
            "filename": file_.name,
            "filesize": os.stat(file_).st_size,
            "linecount": lineno
        }
        rconn.hset(redisargs["fqjobid"],
                   f"file-details:{filetype}:{file_.name}",
                   json.dumps(filedetails))
        return {**filedetails, "errors": errors}


def error_list_name(filetype: str, filename: str):
    """Compute the name of the list where the errors will be pushed.

    Parameters
    ----------
    filetype: The type of file. One of `r_qtl.r_qtl2.FILE_TYPES`
    filename: The name of the file.
    """
    return f"errors:{filetype}:{filename}"


def check_for_geno_errors(
        redisargs: dict[str, str],
        extractdir: Path,
        cdata: dict,
        linesplitterfn: Callable[[str], tuple[Union[str, None]]],
        linejoinerfn: Callable[[tuple[Union[str, None], ...]], str],
        logger: Logger
) -> bool:
    """Check for errors in genotype files."""
    if "geno" in cdata or "founder_geno" in cdata:
        genofiles = tuple(
            extractdir.joinpath(fname) for fname in cdata.get("geno", []))
        fgenofiles = tuple(
            extractdir.joinpath(fname) for fname in cdata.get("founder_geno", []))
        allgenofiles = genofiles + fgenofiles
        with Redis.from_url(redisargs["redisuri"], decode_responses=True) as rconn:
            error_list_names = [
                error_list_name("geno", file_.name) for file_ in allgenofiles]
            for list_name in error_list_names:
                rconn.delete(list_name)
            rconn.hset(
                redisargs["fqjobid"],
                "geno-errors-lists",
                json.dumps(error_list_names))
            processes = [
                mproc.Process(target=file_errors_and_details,
                              args=(
                                  redisargs,
                                  file_,
                                  ftype,
                                  cdata,
                                  linesplitterfn,
                                  linejoinerfn,
                                  (check_markers,),
                                  (check_geno_line,))
                              )
                for ftype, file_ in (
                        tuple(("geno", file_) for file_ in genofiles) +
                        tuple(("founder_geno", file_) for file_ in fgenofiles))
            ]
            for process in processes:
                process.start()
            # Set expiry for any created error lists
            for key in error_list_names:
                rconn.expire(name=key,
                             time=timedelta(seconds=redisargs["redisexpiry"]))

            # TOD0: Add the errors to redis
            if any(rconn.llen(errlst) > 0 for errlst in error_list_names):
                logger.error("At least one of the 'geno' files has (an) error(s).")
                return True
            logger.info("No error(s) found in any of the 'geno' files.")

    else:
        logger.info("No 'geno' files to check.")

    return False


# def check_for_pheno_errors(...):
#     """Check for errors in phenotype files."""
#     pass


# def check_for_phenose_errors(...):
#     """Check for errors in phenotype, standard-error files."""
#     pass


# def check_for_phenocovar_errors(...):
#     """Check for errors in phenotype-covariates files."""
#     pass


def run_qc(rconn: Redis, args: Namespace, fqjobid: str, logger: Logger) -> int:
    """Run quality control checks on R/qtl2 bundles."""
    thejob = parse_job(rconn, args.redisprefix, args.jobid)
    print(f"THE JOB =================> {thejob}")
    jobmeta = thejob["job-metadata"]
    inpath = Path(jobmeta["rqtl2-bundle-file"])
    extractdir = inpath.parent.joinpath(f"{inpath.name}__extraction_dir")
    with ZipFile(inpath, "r") as zfile:
        rqtl2.extract(zfile, extractdir)

    ### BEGIN: The quality control checks ###
    cdata = rqtl2.control_data(extractdir)
    splitter = build_line_splitter(cdata)
    joiner = build_line_joiner(cdata)

    redisargs = {
        "fqjobid": fqjobid,
        "redisuri": args.redisuri,
        "redisexpiry": args.redisexpiry
    }
    check_for_missing_files(rconn, fqjobid, extractdir, logger)
    # check_for_pheno_errors(...)
    check_for_geno_errors(redisargs, extractdir, cdata, splitter, joiner, logger)
    # check_for_phenose_errors(...)
    # check_for_phenocovar_errors(...)
    ### END: The quality control checks ###

    def __fetch_errors__(rkey: str) -> tuple:
        return tuple(json.loads(rconn.hget(fqjobid, rkey) or "[]"))

    return (1 if any((
        bool(__fetch_errors__(key))
        for key in
        ("errors-geno", "errors-pheno", "errors-phenos", "errors-phenocovar")))
            else 0)


if __name__ == "__main__":
    def main():
        """Enter R/qtl2 bundle QC runner."""
        args = add_global_data_arguments(init_cli_parser(
            "qc-on-rqtl2-bundle", "Run QC on R/qtl2 bundle.")).parse_args()
        check_redis(args.redisuri)
        check_db(args.databaseuri)

        logger = getLogger("qc-on-rqtl2-bundle")
        logger.addHandler(StreamHandler(stream=sys.stderr))
        logger.setLevel("DEBUG")

        fqjobid = jobs.job_key(args.redisprefix, args.jobid)
        with Redis.from_url(args.redisuri, decode_responses=True) as rconn:
            logger.addHandler(setup_redis_logger(
                rconn, fqjobid, f"{fqjobid}:log-messages",
                args.redisexpiry))

            exitcode = run_qc(rconn, args, fqjobid, logger)
            rconn.hset(
                jobs.job_key(args.redisprefix, args.jobid), "exitcode", exitcode)
            return exitcode

    sys.exit(main())