flowcraft/templates/integrity_coverage.py
def main(sample_id, fastq_pair, gsize, minimum_coverage, opts):
    """ Main executor of the integrity_coverage template.

    Parameters
    ----------
    sample_id : str
        Sample Identification string.
    fastq_pair : list
        Two element list containing the paired FastQ files.
    gsize : float or int
        Estimate of genome size in Mb.
    minimum_coverage : float or int
        Minimum coverage required for a sample to pass the coverage check
    opts : list
        List of arbitrary options. See `Expected input`_.

    """

    logger.info("Starting integrity coverage main")

    # Check for runtime options
    if "-e" in opts:
        skip_encoding = True
    else:
        skip_encoding = False

    # Information for encoding guess
    gmin, gmax = 99, 0
    encoding = []
    phred = None

    # Information for coverage estimation
    chars = 0
    nreads = 0

    # Information on maximum read length
    max_read_length = 0

    # Get compression of each FastQ pair file
    file_objects = []
    for fastq in fastq_pair:

        logger.info("Processing file {}".format(fastq))

        logger.info("[{}] Guessing file compression".format(fastq))
        ftype = guess_file_compression(fastq)

        # This can guess the compression of gz, bz2 and zip. If it cannot
        # find the compression type, it tries to open a regular file
        if ftype:
            logger.info("[{}] Found file compression: {}".format(
                fastq, ftype))
            file_objects.append(COPEN[ftype](fastq, "rt"))
        else:
            logger.info("[{}] File compression not found. Assuming an "
                        "uncompressed file".format(fastq))
            file_objects.append(open(fastq))

    logger.info("Starting FastQ file parsing")

    # The '*_encoding' file stores a string with the encoding ('Sanger')
    # If no encoding is guessed, 'None' should be stored
    # The '*_phred' file stores a string with the phred score ('33')
    # If no phred is guessed, 'None' should be stored
    # The '*_coverage' file stores the estimated coverage ('88')
    # The '*_report' file stores a csv report of the file
    # The '*_max_len' file stores a string with the maximum contig len ('155')
    with open("{}_encoding".format(sample_id), "w") as enc_fh, \
            open("{}_phred".format(sample_id), "w") as phred_fh, \
            open("{}_coverage".format(sample_id), "w") as cov_fh, \
            open("{}_report".format(sample_id), "w") as cov_rep, \
            open("{}_max_len".format(sample_id), "w") as len_fh, \
            open(".report.json", "w") as json_report, \
            open(".status", "w") as status_fh, \
            open(".fail", "w") as fail_fh:

        try:
            # Iterate over both pair files sequentially using itertools.chain
            for i, line in enumerate(chain(*file_objects)):

                # Parse only every 4th line of the file for the encoding
                # e.g.: AAAA/EEEEEEEEEEE<EEEEEEEEEEEEEEEEEEEEEEEEE (...)
                if (i + 1) % 4 == 0 and not skip_encoding:
                    # It is important to strip() the line so that any newline
                    # character is removed and not accounted for in the
                    # encoding guess
                    lmin, lmax = get_qual_range(line.strip())

                    # Guess new encoding if the range expands the previously
                    # set boundaries of gmin and gmax
                    if lmin < gmin or lmax > gmax:
                        gmin, gmax = min(lmin, gmin), max(lmax, gmax)
                        encoding, phred = get_encodings_in_range(gmin, gmax)
                        logger.debug(
                            "Updating estimates at line {} with range {} to"
                            " '{}' (encoding) and '{}' (phred)".format(
                                i, [lmin, lmax], encoding, phred))

                # Parse only every 2nd line of the file for the coverage
                # e.g.: GGATAATCTACCTTGACGATTTGTACTGGCGTTGGTTTCTTA (...)
                if (i + 3) % 4 == 0:
                    read_len = len(line.strip())
                    chars += read_len
                    nreads += 1

                    # Evaluate maximum read length for sample
                    if read_len > max_read_length:
                        logger.debug("Updating maximum read length at line "
                                     "{} to {}".format(i, read_len))
                        max_read_length = read_len

            # End of FastQ parsing
            logger.info("Finished FastQ file parsing")

            # The minimum expected coverage for a sample to pass
            exp_coverage = round(chars / (gsize * 1e6), 2)

            # Set json report
            if "-e" not in opts:

                json_dic = {
                    "tableRow": [{
                        "sample": sample_id,
                        "data": [
                            {"header": "Raw BP",
                             "value": chars,
                             "table": "qc",
                             "columnBar": True},
                            {"header": "Reads",
                             "value": nreads,
                             "table": "qc",
                             "columnBar": True},
                            {"header": "Coverage",
                             "value": exp_coverage,
                             "table": "qc",
                             "columnBar": True,
                             "failThreshold": minimum_coverage
                             }
                        ]
                    }],
                    "plotData": [{
                        "sample": sample_id,
                        "data": {
                            "sparkline": chars
                        }
                    }],
                }
            else:
                json_dic = {
                    "tableRow": [{
                        "sample": sample_id,
                        "data": [
                            {"header": "Coverage",
                             "value": exp_coverage,
                             "table": "qc",
                             "columnBar": True,
                             "failThreshold": minimum_coverage
                             }
                        ],
                    }],
                }

            # Get encoding
            if len(encoding) > 0:
                encoding = set(encoding)
                phred = set(phred)
                # Get encoding and phred as strings
                # e.g. enc: Sanger, Illumina-1.8
                # e.g. phred: 64
                enc = "{}".format(",".join([x for x in encoding]))
                phred = "{}".format(",".join(str(x) for x in phred))
                logger.info("Encoding set to {}".format(enc))
                logger.info("Phred set to {}".format(enc))

                enc_fh.write(enc)
                phred_fh.write(phred)
            # Encoding not found
            else:
                if not skip_encoding:
                    encoding_msg = "Could not guess encoding and phred from " \
                                   "FastQ"
                    logger.warning(encoding_msg)
                    json_dic["warnings"] = [{
                        "sample": sample_id,
                        "table": "qc",
                        "value": [encoding_msg]
                    }]
                    enc_fh.write("None")
                    phred_fh.write("None")

            # Estimate coverage
            logger.info("Estimating coverage based on a genome size of "
                        "{}".format(gsize))
            logger.info("Expected coverage is {}".format(exp_coverage))

            if exp_coverage >= minimum_coverage:
                cov_rep.write("{},{},{}\\n".format(
                    sample_id, str(exp_coverage), "PASS"))
                cov_fh.write(str(exp_coverage))
                status_fh.write("pass")
            # Estimated coverage does not pass minimum threshold
            else:
                fail_msg = "Sample with low coverage ({}), below the {} " \
                           "threshold".format(exp_coverage, minimum_coverage)
                logger.error(fail_msg)
                fail_fh.write(fail_msg)
                cov_fh.write("fail")
                status_fh.write("fail")
                cov_rep.write("{},{},{}\\n".format(
                    sample_id, str(exp_coverage), "FAIL"))
                json_dic["fail"] = [{
                    "sample": sample_id,
                    "table": "qc",
                    "value": [fail_msg]
                }]

            json_report.write(json.dumps(json_dic, separators=(",", ":")))
            # Maximum read length
            len_fh.write("{}".format(max_read_length))

        # This exception is raised when the input FastQ files are corrupted
        except EOFError:
            logger.error("The FastQ files could not be correctly "
                         "parsed. They may be corrupt")
            for fh in [enc_fh, phred_fh, cov_fh, cov_rep, len_fh]:
                fh.write("corrupt")
                status_fh.write("fail")
                fail_fh.write("Could not read/parse FastQ. "
                              "Possibly corrupt file")
def main(sample_id, fastq_pair, gsize, minimum_coverage, opts):
    """ Main executor of the integrity_coverage template.

    Parameters
    ----------
    sample_id : str
        Sample Identification string.
    fastq_pair : list
        Two element list containing the paired FastQ files.
    gsize : float or int
        Estimate of genome size in Mb.
    minimum_coverage : float or int
        Minimum coverage required for a sample to pass the coverage check
    opts : list
        List of arbitrary options. See `Expected input`_.

    """

    logger.info("Starting integrity coverage main")

    # Check for runtime options
    if "-e" in opts:
        skip_encoding = True
    else:
        skip_encoding = False

    # Information for encoding guess
    gmin, gmax = 99, 0
    encoding = []
    phred = None

    # Information for coverage estimation
    chars = 0
    nreads = 0

    # Information on maximum read length
    max_read_length = 0

    # Get compression of each FastQ pair file
    file_objects = []
    for fastq in fastq_pair:

        logger.info("Processing file {}".format(fastq))

        logger.info("[{}] Guessing file compression".format(fastq))
        ftype = guess_file_compression(fastq)

        # This can guess the compression of gz, bz2 and zip. If it cannot
        # find the compression type, it tries to open a regular file
        if ftype:
            logger.info("[{}] Found file compression: {}".format(
                fastq, ftype))
            file_objects.append(COPEN[ftype](fastq, "rt"))
        else:
            logger.info("[{}] File compression not found. Assuming an "
                        "uncompressed file".format(fastq))
            file_objects.append(open(fastq))

    logger.info("Starting FastQ file parsing")

    # The '*_encoding' file stores a string with the encoding ('Sanger')
    # If no encoding is guessed, 'None' should be stored
    # The '*_phred' file stores a string with the phred score ('33')
    # If no phred is guessed, 'None' should be stored
    # The '*_coverage' file stores the estimated coverage ('88')
    # The '*_report' file stores a csv report of the file
    # The '*_max_len' file stores a string with the maximum contig len ('155')
    with open("{}_encoding".format(sample_id), "w") as enc_fh, \
            open("{}_phred".format(sample_id), "w") as phred_fh, \
            open("{}_coverage".format(sample_id), "w") as cov_fh, \
            open("{}_report".format(sample_id), "w") as cov_rep, \
            open("{}_max_len".format(sample_id), "w") as len_fh, \
            open(".report.json", "w") as json_report, \
            open(".status", "w") as status_fh, \
            open(".fail", "w") as fail_fh:

        try:
            # Iterate over both pair files sequentially using itertools.chain
            for i, line in enumerate(chain(*file_objects)):

                # Parse only every 4th line of the file for the encoding
                # e.g.: AAAA/EEEEEEEEEEE<EEEEEEEEEEEEEEEEEEEEEEEEE (...)
                if (i + 1) % 4 == 0 and not skip_encoding:
                    # It is important to strip() the line so that any newline
                    # character is removed and not accounted for in the
                    # encoding guess
                    lmin, lmax = get_qual_range(line.strip())

                    # Guess new encoding if the range expands the previously
                    # set boundaries of gmin and gmax
                    if lmin < gmin or lmax > gmax:
                        gmin, gmax = min(lmin, gmin), max(lmax, gmax)
                        encoding, phred = get_encodings_in_range(gmin, gmax)
                        logger.debug(
                            "Updating estimates at line {} with range {} to"
                            " '{}' (encoding) and '{}' (phred)".format(
                                i, [lmin, lmax], encoding, phred))

                # Parse only every 2nd line of the file for the coverage
                # e.g.: GGATAATCTACCTTGACGATTTGTACTGGCGTTGGTTTCTTA (...)
                if (i + 3) % 4 == 0:
                    read_len = len(line.strip())
                    chars += read_len
                    nreads += 1

                    # Evaluate maximum read length for sample
                    if read_len > max_read_length:
                        logger.debug("Updating maximum read length at line "
                                     "{} to {}".format(i, read_len))
                        max_read_length = read_len

            # End of FastQ parsing
            logger.info("Finished FastQ file parsing")

            # The minimum expected coverage for a sample to pass
            exp_coverage = round(chars / (gsize * 1e6), 2)

            # Set json report
            if "-e" not in opts:

                json_dic = {
                    "tableRow": [{
                        "sample": sample_id,
                        "data": [
                            {"header": "Raw BP",
                             "value": chars,
                             "table": "qc",
                             "columnBar": True},
                            {"header": "Reads",
                             "value": nreads,
                             "table": "qc",
                             "columnBar": True},
                            {"header": "Coverage",
                             "value": exp_coverage,
                             "table": "qc",
                             "columnBar": True,
                             "failThreshold": minimum_coverage
                             }
                        ]
                    }],
                    "plotData": [{
                        "sample": sample_id,
                        "data": {
                            "sparkline": chars
                        }
                    }],
                }
            else:
                json_dic = {
                    "tableRow": [{
                        "sample": sample_id,
                        "data": [
                            {"header": "Coverage",
                             "value": exp_coverage,
                             "table": "qc",
                             "columnBar": True,
                             "failThreshold": minimum_coverage
                             }
                        ],
                    }],
                }

            # Get encoding
            if len(encoding) > 0:
                encoding = set(encoding)
                phred = set(phred)
                # Get encoding and phred as strings
                # e.g. enc: Sanger, Illumina-1.8
                # e.g. phred: 64
                enc = "{}".format(",".join([x for x in encoding]))
                phred = "{}".format(",".join(str(x) for x in phred))
                logger.info("Encoding set to {}".format(enc))
                logger.info("Phred set to {}".format(enc))

                enc_fh.write(enc)
                phred_fh.write(phred)
            # Encoding not found
            else:
                if not skip_encoding:
                    encoding_msg = "Could not guess encoding and phred from " \
                                   "FastQ"
                    logger.warning(encoding_msg)
                    json_dic["warnings"] = [{
                        "sample": sample_id,
                        "table": "qc",
                        "value": [encoding_msg]
                    }]
                    enc_fh.write("None")
                    phred_fh.write("None")

            # Estimate coverage
            logger.info("Estimating coverage based on a genome size of "
                        "{}".format(gsize))
            logger.info("Expected coverage is {}".format(exp_coverage))

            if exp_coverage >= minimum_coverage:
                cov_rep.write("{},{},{}\\n".format(
                    sample_id, str(exp_coverage), "PASS"))
                cov_fh.write(str(exp_coverage))
                status_fh.write("pass")
            # Estimated coverage does not pass minimum threshold
            else:
                fail_msg = "Sample with low coverage ({}), below the {} " \
                           "threshold".format(exp_coverage, minimum_coverage)
                logger.error(fail_msg)
                fail_fh.write(fail_msg)
                cov_fh.write("fail")
                status_fh.write("fail")
                cov_rep.write("{},{},{}\\n".format(
                    sample_id, str(exp_coverage), "FAIL"))
                json_dic["fail"] = [{
                    "sample": sample_id,
                    "table": "qc",
                    "value": [fail_msg]
                }]

            json_report.write(json.dumps(json_dic, separators=(",", ":")))
            # Maximum read length
            len_fh.write("{}".format(max_read_length))

        # This exception is raised when the input FastQ files are corrupted
        except EOFError:
            logger.error("The FastQ files could not be correctly "
                         "parsed. They may be corrupt")
            for fh in [enc_fh, phred_fh, cov_fh, cov_rep, len_fh]:
                fh.write("corrupt")
                status_fh.write("fail")
                fail_fh.write("Could not read/parse FastQ. "
                              "Possibly corrupt file")
Main executor of the integrity_coverage template.

    Parameters
    ----------
    sample_id : str
        Sample Identification string.
    fastq_pair : list
        Two element list containing the paired FastQ files.
    gsize : float or int
        Estimate of genome size in Mb.
    minimum_coverage : float or int
        Minimum coverage required for a sample to pass the coverage check
    opts : list
        List of arbitrary options. See `Expected input`_.