flowcraft/templates/process_assembly.py
def main(sample_id, assembly_file, gsize, opts, assembler):
    """Main executor of the process_spades template.

    Parameters
    ----------
    sample_id : str
        Sample Identification string.
    assembly_file : str
        Path to the assembly file generated by Spades.
    gsize : int
        Estimate of genome size.
    opts : list
        List of options for processing spades assembly.
    assembler : str
        Name of the assembler, for logging purposes

    """

    logger.info("Starting assembly file processing")
    warnings = []
    fails = ""

    min_contig_len, min_kmer_cov, max_contigs = [int(x) for x in opts]
    logger.debug("Setting minimum conting length to: {}".format(
        min_contig_len))
    logger.debug("Setting minimum kmer coverage: {}".format(min_kmer_cov))

    # Parse the spades assembly file and perform the first filtering.
    logger.info("Starting assembly parsing")
    assembly_obj = Assembly(assembly_file, min_contig_len, min_kmer_cov,
                               sample_id)

    with open(".warnings", "w") as warn_fh:
        t_80 = gsize * 1000000 * 0.8
        t_150 = gsize * 1000000 * 1.5
        # Check if assembly size of the first assembly is lower than 80% of the
        # estimated genome size. If True, redo the filtering without the
        # k-mer coverage filter
        assembly_len = assembly_obj.get_assembly_length()
        logger.debug("Checking assembly length: {}".format(assembly_len))

        if assembly_len < t_80:

            logger.warning("Assembly size ({}) smaller than the minimum "
                           "threshold of 80% of expected genome size. "
                           "Applying contig filters without the k-mer "
                           "coverage filter".format(assembly_len))
            assembly_obj.filter_contigs(*[
                ["length", ">=", min_contig_len]
            ])

            assembly_len = assembly_obj.get_assembly_length()
            logger.debug("Checking updated assembly length: "
                         "{}".format(assembly_len))
            if assembly_len < t_80:

                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 = warn_msg

        if assembly_len > t_150:

            warn_msg = "Assembly size ({}) larger than the maximum" \
                       " threshold of 150% of expected genome size.".format(
                            assembly_len)
            logger.warning(warn_msg)
            warn_fh.write(warn_msg)
            fails = warn_msg

        logger.debug("Checking number of contigs: {}".format(
            len(assembly_obj.contigs)))
        contig_threshold = (max_contigs * gsize) / 1.5
        if len(assembly_obj.contigs) > contig_threshold:

            warn_msg = "The number of contigs ({}) exceeds the threshold of " \
                       "{} contigs per 1.5Mb ({})".format(
                            len(assembly_obj.contigs),
                            max_contigs,
                            round(contig_threshold, 1))

            logger.warning(warn_msg)
            warn_fh.write(warn_msg)
            warnings.append(warn_msg)

    # Write filtered assembly
    logger.debug("Renaming old assembly file to: {}".format(
        "{}.old".format(assembly_file)))
    assembly_obj.write_assembly("{}_proc.fasta".format(
        os.path.splitext(assembly_file)[0]))
    # Write report
    output_report = "{}.report.csv".format(sample_id)
    assembly_obj.write_report(output_report)
    # Write json report
    with open(".report.json", "w") as json_report:
        json_dic = {
            "tableRow": [{
                "sample": sample_id,
                "data": [
                    {"header": "Contigs ({})".format(assembler),
                     "value": len(assembly_obj.contigs),
                     "table": "assembly",
                     "columnBar": True},
                    {"header": "Assembled BP ({})".format(assembler),
                     "value": assembly_len,
                     "table": "assembly",
                     "columnBar": True}
                ]
            }],
        }

        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=(",", ":")))

    with open(".status", "w") as status_fh:
        status_fh.write("pass")
def main(sample_id, assembly_file, gsize, opts, assembler):
    """Main executor of the process_spades template.

    Parameters
    ----------
    sample_id : str
        Sample Identification string.
    assembly_file : str
        Path to the assembly file generated by Spades.
    gsize : int
        Estimate of genome size.
    opts : list
        List of options for processing spades assembly.
    assembler : str
        Name of the assembler, for logging purposes

    """

    logger.info("Starting assembly file processing")
    warnings = []
    fails = ""

    min_contig_len, min_kmer_cov, max_contigs = [int(x) for x in opts]
    logger.debug("Setting minimum conting length to: {}".format(
        min_contig_len))
    logger.debug("Setting minimum kmer coverage: {}".format(min_kmer_cov))

    # Parse the spades assembly file and perform the first filtering.
    logger.info("Starting assembly parsing")
    assembly_obj = Assembly(assembly_file, min_contig_len, min_kmer_cov,
                               sample_id)

    with open(".warnings", "w") as warn_fh:
        t_80 = gsize * 1000000 * 0.8
        t_150 = gsize * 1000000 * 1.5
        # Check if assembly size of the first assembly is lower than 80% of the
        # estimated genome size. If True, redo the filtering without the
        # k-mer coverage filter
        assembly_len = assembly_obj.get_assembly_length()
        logger.debug("Checking assembly length: {}".format(assembly_len))

        if assembly_len < t_80:

            logger.warning("Assembly size ({}) smaller than the minimum "
                           "threshold of 80% of expected genome size. "
                           "Applying contig filters without the k-mer "
                           "coverage filter".format(assembly_len))
            assembly_obj.filter_contigs(*[
                ["length", ">=", min_contig_len]
            ])

            assembly_len = assembly_obj.get_assembly_length()
            logger.debug("Checking updated assembly length: "
                         "{}".format(assembly_len))
            if assembly_len < t_80:

                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 = warn_msg

        if assembly_len > t_150:

            warn_msg = "Assembly size ({}) larger than the maximum" \
                       " threshold of 150% of expected genome size.".format(
                            assembly_len)
            logger.warning(warn_msg)
            warn_fh.write(warn_msg)
            fails = warn_msg

        logger.debug("Checking number of contigs: {}".format(
            len(assembly_obj.contigs)))
        contig_threshold = (max_contigs * gsize) / 1.5
        if len(assembly_obj.contigs) > contig_threshold:

            warn_msg = "The number of contigs ({}) exceeds the threshold of " \
                       "{} contigs per 1.5Mb ({})".format(
                            len(assembly_obj.contigs),
                            max_contigs,
                            round(contig_threshold, 1))

            logger.warning(warn_msg)
            warn_fh.write(warn_msg)
            warnings.append(warn_msg)

    # Write filtered assembly
    logger.debug("Renaming old assembly file to: {}".format(
        "{}.old".format(assembly_file)))
    assembly_obj.write_assembly("{}_proc.fasta".format(
        os.path.splitext(assembly_file)[0]))
    # Write report
    output_report = "{}.report.csv".format(sample_id)
    assembly_obj.write_report(output_report)
    # Write json report
    with open(".report.json", "w") as json_report:
        json_dic = {
            "tableRow": [{
                "sample": sample_id,
                "data": [
                    {"header": "Contigs ({})".format(assembler),
                     "value": len(assembly_obj.contigs),
                     "table": "assembly",
                     "columnBar": True},
                    {"header": "Assembled BP ({})".format(assembler),
                     "value": assembly_len,
                     "table": "assembly",
                     "columnBar": True}
                ]
            }],
        }

        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=(",", ":")))

    with open(".status", "w") as status_fh:
        status_fh.write("pass")
Main executor of the process_spades template.

    Parameters
    ----------
    sample_id : str
        Sample Identification string.
    assembly_file : str
        Path to the assembly file generated by Spades.
    gsize : int
        Estimate of genome size.
    opts : list
        List of options for processing spades assembly.
    assembler : str
        Name of the assembler, for logging purposes