flowcraft/templates/process_assembly_mapping.py
def check_filtered_assembly(coverage_info, coverage_bp, minimum_coverage,
                            genome_size, contig_size, max_contigs,
                            sample_id):
    """Checks whether a filtered assembly passes a size threshold

    Given a minimum coverage threshold, this function evaluates whether an
    assembly will pass the minimum threshold of ``genome_size * 1e6 * 0.8``,
    which means 80% of the expected genome size or the maximum threshold
    of ``genome_size * 1e6 * 1.5``, which means 150% of the expected genome
    size. It will issue a warning if any of these thresholds is crossed.
    In the case of an expected genome size below 80% it will return False.

    Parameters
    ----------
    coverage_info : OrderedDict or dict
        Dictionary containing the coverage information for each contig.
    coverage_bp : dict
        Dictionary containing the per base coverage information for each
        contig. Used to determine the total number of base pairs in the
        final assembly.
    minimum_coverage : int
        Minimum coverage required for a contig to pass the filter.
    genome_size : int
        Expected genome size.
    contig_size : dict
        Dictionary with the len of each contig. Contig headers as keys and
        the corresponding lenght as values.
    max_contigs : int
        Maximum threshold for contig number. A warning is issued if this
        threshold is crossed.
    sample_id : str
        Id or name of the current sample

    Returns
    -------
    x : bool
        True if the filtered assembly size is higher than 80% of the
        expected genome size.

    """

    # Get size of assembly after filtering contigs below minimum_coverage
    assembly_len = sum([v for k, v in contig_size.items()
                        if coverage_info[k]["cov"] >= minimum_coverage])
    logger.debug("Assembly length after filtering for minimum coverage of"
                 " {}: {}".format(minimum_coverage, assembly_len))
    # Get number of contigs after filtering
    ncontigs = len([x for x in coverage_info.values()
                    if x["cov"] >= minimum_coverage])
    logger.debug("Number of contigs: {}".format(ncontigs))
    # Get number of bp after filtering
    filtered_contigs = [k for k, v in coverage_info.items()
                        if v["cov"] >= minimum_coverage]
    logger.debug("Filtered contigs for minimum coverage of "
                 "{}: {}".format(minimum_coverage, filtered_contigs))
    total_assembled_bp = sum([sum(coverage_bp[x]) for x in filtered_contigs
                              if x in coverage_bp])
    logger.debug("Total number of assembled base pairs:"
                 "{}".format(total_assembled_bp))

    warnings = []
    fails = []
    health = True

    with open(".warnings", "w") as warn_fh, \
            open(".report.json", "w") as json_report:

        logger.debug("Checking assembly size after filtering : {}".format(
            assembly_len))

        # If the filtered assembly size is above the 150% genome size
        # threshold, issue a warning
        if assembly_len > genome_size * 1e6 * 1.5:
            warn_msg = "Assembly size ({}) smaller than the maximum" \
                       " threshold of 150% of expected genome size.".format(
                            assembly_len)
            logger.warning(warn_msg)
            warn_fh.write(warn_msg)
            fails.append("Large_genome_size_({})".format(assembly_len))

        # If the number of contigs in the filtered assembly size crosses the
        # max_contigs threshold, issue a warning
        logger.debug("Checking number of contigs: {}".format(
                len(coverage_info)))
        contig_threshold = max_contigs * genome_size / 1.5
        if ncontigs > contig_threshold:
            warn_msg = "The number of contigs ({}) exceeds the threshold of " \
                       "100 contigs per 1.5Mb ({})".format(
                            ncontigs, round(contig_threshold, 1))
            logger.warning(warn_msg)
            warn_fh.write(warn_msg)
            warnings.append(warn_msg)

        # If the filtered assembly size falls below the 80% genome size
        # threshold, fail this check and return False
        if assembly_len < genome_size * 1e6 * 0.8:
            warn_msg = "Assembly size smaller than the minimum" \
                       " threshold of 80% of expected genome size: {}".format(
                            assembly_len)
            logger.warning(warn_msg)
            warn_fh.write(warn_msg)
            fails.append("Small_genome_size_({})".format(assembly_len))
            assembly_len = sum([v for v in contig_size.values()])
            total_assembled_bp = sum(
                [sum(coverage_bp[x]) for x in coverage_info if x in
                 coverage_bp])
            logger.debug("Assembly length without coverage filtering: "
                         "{}".format(assembly_len))
            logger.debug("Total number of assembled base pairs without"
                         " filtering: {}".format(total_assembled_bp))

            health = False

        json_dic = {
            "plotData": [{
                "sample": sample_id,
                "data": {
                    "sparkline": total_assembled_bp
                }
            }]
        }

        if warnings:
            json_dic["warnings"] = [{
                "sample": sample_id,
                "table": "assembly",
                "value": warnings
            }]
        if fails:
            json_dic["fail"] = [{
                "sample": sample_id,
                "table": "assembly",
                "value": [fails]
            }]

        json_report.write(json.dumps(json_dic, separators=(",", ":")))

    return health
def check_filtered_assembly(coverage_info, coverage_bp, minimum_coverage,
                            genome_size, contig_size, max_contigs,
                            sample_id):
    """Checks whether a filtered assembly passes a size threshold

    Given a minimum coverage threshold, this function evaluates whether an
    assembly will pass the minimum threshold of ``genome_size * 1e6 * 0.8``,
    which means 80% of the expected genome size or the maximum threshold
    of ``genome_size * 1e6 * 1.5``, which means 150% of the expected genome
    size. It will issue a warning if any of these thresholds is crossed.
    In the case of an expected genome size below 80% it will return False.

    Parameters
    ----------
    coverage_info : OrderedDict or dict
        Dictionary containing the coverage information for each contig.
    coverage_bp : dict
        Dictionary containing the per base coverage information for each
        contig. Used to determine the total number of base pairs in the
        final assembly.
    minimum_coverage : int
        Minimum coverage required for a contig to pass the filter.
    genome_size : int
        Expected genome size.
    contig_size : dict
        Dictionary with the len of each contig. Contig headers as keys and
        the corresponding lenght as values.
    max_contigs : int
        Maximum threshold for contig number. A warning is issued if this
        threshold is crossed.
    sample_id : str
        Id or name of the current sample

    Returns
    -------
    x : bool
        True if the filtered assembly size is higher than 80% of the
        expected genome size.

    """

    # Get size of assembly after filtering contigs below minimum_coverage
    assembly_len = sum([v for k, v in contig_size.items()
                        if coverage_info[k]["cov"] >= minimum_coverage])
    logger.debug("Assembly length after filtering for minimum coverage of"
                 " {}: {}".format(minimum_coverage, assembly_len))
    # Get number of contigs after filtering
    ncontigs = len([x for x in coverage_info.values()
                    if x["cov"] >= minimum_coverage])
    logger.debug("Number of contigs: {}".format(ncontigs))
    # Get number of bp after filtering
    filtered_contigs = [k for k, v in coverage_info.items()
                        if v["cov"] >= minimum_coverage]
    logger.debug("Filtered contigs for minimum coverage of "
                 "{}: {}".format(minimum_coverage, filtered_contigs))
    total_assembled_bp = sum([sum(coverage_bp[x]) for x in filtered_contigs
                              if x in coverage_bp])
    logger.debug("Total number of assembled base pairs:"
                 "{}".format(total_assembled_bp))

    warnings = []
    fails = []
    health = True

    with open(".warnings", "w") as warn_fh, \
            open(".report.json", "w") as json_report:

        logger.debug("Checking assembly size after filtering : {}".format(
            assembly_len))

        # If the filtered assembly size is above the 150% genome size
        # threshold, issue a warning
        if assembly_len > genome_size * 1e6 * 1.5:
            warn_msg = "Assembly size ({}) smaller than the maximum" \
                       " threshold of 150% of expected genome size.".format(
                            assembly_len)
            logger.warning(warn_msg)
            warn_fh.write(warn_msg)
            fails.append("Large_genome_size_({})".format(assembly_len))

        # If the number of contigs in the filtered assembly size crosses the
        # max_contigs threshold, issue a warning
        logger.debug("Checking number of contigs: {}".format(
                len(coverage_info)))
        contig_threshold = max_contigs * genome_size / 1.5
        if ncontigs > contig_threshold:
            warn_msg = "The number of contigs ({}) exceeds the threshold of " \
                       "100 contigs per 1.5Mb ({})".format(
                            ncontigs, round(contig_threshold, 1))
            logger.warning(warn_msg)
            warn_fh.write(warn_msg)
            warnings.append(warn_msg)

        # If the filtered assembly size falls below the 80% genome size
        # threshold, fail this check and return False
        if assembly_len < genome_size * 1e6 * 0.8:
            warn_msg = "Assembly size smaller than the minimum" \
                       " threshold of 80% of expected genome size: {}".format(
                            assembly_len)
            logger.warning(warn_msg)
            warn_fh.write(warn_msg)
            fails.append("Small_genome_size_({})".format(assembly_len))
            assembly_len = sum([v for v in contig_size.values()])
            total_assembled_bp = sum(
                [sum(coverage_bp[x]) for x in coverage_info if x in
                 coverage_bp])
            logger.debug("Assembly length without coverage filtering: "
                         "{}".format(assembly_len))
            logger.debug("Total number of assembled base pairs without"
                         " filtering: {}".format(total_assembled_bp))

            health = False

        json_dic = {
            "plotData": [{
                "sample": sample_id,
                "data": {
                    "sparkline": total_assembled_bp
                }
            }]
        }

        if warnings:
            json_dic["warnings"] = [{
                "sample": sample_id,
                "table": "assembly",
                "value": warnings
            }]
        if fails:
            json_dic["fail"] = [{
                "sample": sample_id,
                "table": "assembly",
                "value": [fails]
            }]

        json_report.write(json.dumps(json_dic, separators=(",", ":")))

    return health
def check_filtered_assembly ( coverage_info , coverage_bp , minimum_coverage , genome_size , contig_size , max_contigs , sample_id ) : # Get size of assembly after filtering contigs below minimum_coverage assembly_len = sum ( [ v for k , v in contig_size . items ( ) if coverage_info [ k ] [ "cov" ] >= minimum_coverage ] ) logger . debug ( "Assembly length after filtering for minimum coverage of" " {}: {}" . format ( minimum_coverage , assembly_len ) ) # Get number of contigs after filtering ncontigs = len ( [ x for x in coverage_info . values ( ) if x [ "cov" ] >= minimum_coverage ] ) logger . debug ( "Number of contigs: {}" . format ( ncontigs ) ) # Get number of bp after filtering filtered_contigs = [ k for k , v in coverage_info . items ( ) if v [ "cov" ] >= minimum_coverage ] logger . debug ( "Filtered contigs for minimum coverage of " "{}: {}" . format ( minimum_coverage , filtered_contigs ) ) total_assembled_bp = sum ( [ sum ( coverage_bp [ x ] ) for x in filtered_contigs if x in coverage_bp ] ) logger . debug ( "Total number of assembled base pairs:" "{}" . format ( total_assembled_bp ) ) warnings = [ ] fails = [ ] health = True with open ( ".warnings" , "w" ) as warn_fh , open ( ".report.json" , "w" ) as json_report : logger . debug ( "Checking assembly size after filtering : {}" . format ( assembly_len ) ) # If the filtered assembly size is above the 150% genome size # threshold, issue a warning if assembly_len > genome_size * 1e6 * 1.5 : warn_msg = "Assembly size ({}) smaller than the maximum" " threshold of 150% of expected genome size." . format ( assembly_len ) logger . warning ( warn_msg ) warn_fh . write ( warn_msg ) fails . append ( "Large_genome_size_({})" . format ( assembly_len ) ) # If the number of contigs in the filtered assembly size crosses the # max_contigs threshold, issue a warning logger . debug ( "Checking number of contigs: {}" . format ( len ( coverage_info ) ) ) contig_threshold = max_contigs * genome_size / 1.5 if ncontigs > contig_threshold : warn_msg = "The number of contigs ({}) exceeds the threshold of " "100 contigs per 1.5Mb ({})" . format ( ncontigs , round ( contig_threshold , 1 ) ) logger . warning ( warn_msg ) warn_fh . write ( warn_msg ) warnings . append ( warn_msg ) # If the filtered assembly size falls below the 80% genome size # threshold, fail this check and return False if assembly_len < genome_size * 1e6 * 0.8 : warn_msg = "Assembly size smaller than the minimum" " threshold of 80% of expected genome size: {}" . format ( assembly_len ) logger . warning ( warn_msg ) warn_fh . write ( warn_msg ) fails . append ( "Small_genome_size_({})" . format ( assembly_len ) ) assembly_len = sum ( [ v for v in contig_size . values ( ) ] ) total_assembled_bp = sum ( [ sum ( coverage_bp [ x ] ) for x in coverage_info if x in coverage_bp ] ) logger . debug ( "Assembly length without coverage filtering: " "{}" . format ( assembly_len ) ) logger . debug ( "Total number of assembled base pairs without" " filtering: {}" . format ( total_assembled_bp ) ) health = False json_dic = { "plotData" : [ { "sample" : sample_id , "data" : { "sparkline" : total_assembled_bp } } ] } if warnings : json_dic [ "warnings" ] = [ { "sample" : sample_id , "table" : "assembly" , "value" : warnings } ] if fails : json_dic [ "fail" ] = [ { "sample" : sample_id , "table" : "assembly" , "value" : [ fails ] } ] json_report . write ( json . dumps ( json_dic , separators = ( "," , ":" ) ) ) return health
Checks whether a filtered assembly passes a size threshold

    Given a minimum coverage threshold, this function evaluates whether an
    assembly will pass the minimum threshold of ``genome_size * 1e6 * 0.8``,
    which means 80% of the expected genome size or the maximum threshold
    of ``genome_size * 1e6 * 1.5``, which means 150% of the expected genome
    size. It will issue a warning if any of these thresholds is crossed.
    In the case of an expected genome size below 80% it will return False.

    Parameters
    ----------
    coverage_info : OrderedDict or dict
        Dictionary containing the coverage information for each contig.
    coverage_bp : dict
        Dictionary containing the per base coverage information for each
        contig. Used to determine the total number of base pairs in the
        final assembly.
    minimum_coverage : int
        Minimum coverage required for a contig to pass the filter.
    genome_size : int
        Expected genome size.
    contig_size : dict
        Dictionary with the len of each contig. Contig headers as keys and
        the corresponding lenght as values.
    max_contigs : int
        Maximum threshold for contig number. A warning is issued if this
        threshold is crossed.
    sample_id : str
        Id or name of the current sample

    Returns
    -------
    x : bool
        True if the filtered assembly size is higher than 80% of the
        expected genome size.
Checks whether a filtered assembly passes a size threshold