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")
def main ( sample_id , assembly_file , gsize , opts , assembler ) : 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
Main executor of the process_spades template .