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")
def main ( sample_id , fastq_pair , gsize , minimum_coverage , opts ) : 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`_.
Main executor of the integrity_coverage template .