comparison graphprot_predict_wrapper.py @ 5:ddcf35a868b8 draft

planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/tools/rna_tools/graphprot commit ad60258f5759eaa205fec4af6143c728ea131419

4:4ad83aed5c3c

import os
import subprocess
import sys

import gplib


"""

TOOL DEPENDENCIES
=================

--sc-thr 0.0 --max-merge-dist 0 --conf-out  --ap-extlr 5

"""


###############################################################################

def setup_argument_parser():
    """Setup argparse parser."""
    help_description = """
    Galaxy wrapper script for GraphProt (-action predict and -action

    If --gen-site-bed .bed file is provided, peak regions will be output
    with genomic coordinates too.

    """
    # Define argument parser.
    p = ap.ArgumentParser(add_help=False,
                          prog="graphprot_predict_wrapper.py",
                          description=help_description,
                          formatter_class=ap.MetavarTypeHelpFormatter)

    # Argument groups.
    p_man = p.add_argument_group("REQUIRED ARGUMENTS")
    p_opt = p.add_argument_group("OPTIONAL ARGUMENTS")

    # Required arguments.
    p_opt.add_argument("-h", "--help",
                       action="help",
                       help="Print help message")
    p_man.add_argument("--fasta",
                       dest="in_fa",
                       type=str,
                       required=True,
                       help="Sequences .fa file to predict"
                            " on (option -fasta)")
    p_man.add_argument("--model",
                       dest="in_model",
                       type=str,
                       required=True,
                       help="GraphProt model file to use for predictions"
                            " (option -model)")
    p_man.add_argument("--params",
                       dest="in_params",
                       type=str,
                       required=True,
                       help="Parameter file for given model")
    p_man.add_argument("--data-id",
                       dest="data_id",
                       type=str,
                       required=True,
                       help="Data ID (option -prefix)")
    # ---> I'm  a conditional argument <---
    p_opt.add_argument("--ws-pred",
                       dest="ws_pred",
                       default=False,
                       action="store_true",
                       help="Run a whole site prediction instead "
                            "of calculating profiles (default: false)")
    # Additional arguments.
    p_opt.add_argument("--sc-thr",
                       dest="score_thr",
                       type=float,
                       default=0,
                       help="Score threshold for extracting "
                            "average profile peak regions (default: 0)")
    p_opt.add_argument("--max-merge-dist",
                       dest="max_merge_dist",
                       type=int,
                       default=0,
                       choices=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
                       help="Maximum merge distance for nearby peak regions"
                            " (default: report all non-overlapping regions)")
    p_opt.add_argument("--gen-site-bed",
                       dest="genomic_sites_bed",
                       type=str,
                       help=".bed file specifying the genomic regions of "
                            "the input .fa sequences. Corrupt .bed "
                            "information will be punished (default: false)")
    p_opt.add_argument("--conf-out",
                       dest="conf_out",
                       default=False,
                       action="store_true",
                       help="Output filtered peak regions BED file or "
                            "predictions file (if --ws-pred) using the median "
                            "positive training site score for filtering "
                            "(default: false)")
    p_opt.add_argument("--gp-output",
                       dest="gp_output",
                       default=False,
                       action="store_true",
                       help="Print output produced by GraphProt "
                            "(default: false)")
    p_opt.add_argument("--ap-extlr",
                       dest="ap_extlr",
                       type=int,
                       default=5,
                       choices=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
                       help="Define average profile up- and downstream "
                            "extension to produce the average profile. "
                            "The mean over small sequence windows "
                            "(window length = --ap-extlr*2 + 1) is used to "
                            "get position scores, thus the average profile "
                            "is more smooth than the initial profile output "
                            "by GraphProt (default: 5)")
    return p


###############################################################################

if __name__ == '__main__':

    # Setup argparse.
    parser = setup_argument_parser()
    # Read in command line arguments.
    args = parser.parse_args()

    # Check for Linux.
    assert "linux" in sys.platform, "please use Linux"
    # Check tool availability.
    assert gplib.is_tool("GraphProt.pl"), "GraphProt.pl not in PATH"
    # Check file inputs.
    assert os.path.exists(args.in_fa), \
        "input .fa file \"%s\" not found" % (args.in_fa)
    assert os.path.exists(args.in_model), \
        "input .model file \"%s\" not found" % (args.in_model)
    assert os.path.exists(args.in_params), \
        "input .params file \"%s\" not found" % (args.in_params)
    # Count .fa entries.
    c_in_fa = gplib.count_fasta_headers(args.in_fa)
    assert c_in_fa, "input .fa file \"%s\" no headers found" % (args.in_fa)
    print("# input .fa sequences:   %i" % (c_in_fa))
    # Read in FASTA sequences to check for uppercase sequences.
    seqs_dic = gplib.read_fasta_into_dic(args.in_fa)
    # Check for lowercase only sequences, which cause GP to crash.
    error_mess = "input sequences encountered containing "\
        "only lowercase characters or lowercase characters in between "\
        "uppercase characters. Please provide either all uppercase "\
        "sequences or sequences containing uppercase regions surrounded "\
        "by lowercase context regions for structure calculation (see "\
        "viewpoint concept in original GraphProt publication "\
        "for more details)"
    if args.ws_pred:
        bad_ids = gplib.check_seqs_dic_format(seqs_dic)
        assert not bad_ids, "%s" % (error_mess)

    c_uc_nt = gplib.seqs_dic_count_uc_nts(seqs_dic)
    assert c_uc_nt, "no uppercase nucleotides in input .fa sequences. "\
                    "Please change sequences to uppercase (keep in mind "\
                    "GraphProt only scores uppercase regions (according "\
                    "to its viewpoint concept)"
    if not args.ws_pred:
        # Check for lowercase sequences.
        c_lc_nt = gplib.seqs_dic_count_lc_nts(seqs_dic)
        assert not c_lc_nt, "lowercase nucleotides in input .fa not allowed"\
            " in profile predictions, since GraphProt only scores "\
            "uppercase regions (according to its viewpoint concept)"
    # Check .bed.
    if args.genomic_sites_bed:
        # An array of checks, marvelous.
        assert \
            os.path.exists(args.genomic_sites_bed), \
            "genomic .bed file \"%s\" not found" % (args.genomic_sites_bed)
        # Check .bed for content.
        assert gplib.count_file_rows(args.genomic_sites_bed), \
            "genomic .bed file \"%s\" is empty" % (args.genomic_sites_bed)
        # Check .bed for 6-column format.
        assert gplib.bed_check_six_col_format(args.genomic_sites_bed), \
            "genomic .bed file \"%s\" appears to not "\
            "be in 6-column .bed format" % (args.genomic_sites_bed)
        # Check for unique column 4 IDs.
        assert gplib.bed_check_unique_ids(args.genomic_sites_bed), \
            "genomic .bed file \"%s\" column 4 IDs not unique" % \
            (args.genomic_sites_bed)
        # Read in .bed regions, compare to FASTA sequences.
        seq_len_dic = gplib.get_seq_lengths_from_seqs_dic(seqs_dic)
        reg_len_dic = gplib.bed_get_region_lengths(args.genomic_sites_bed)
        for seq_id in seq_len_dic:
            seq_l = seq_len_dic[seq_id]
            assert seq_id in reg_len_dic, \
                "sequence ID \"%s\" missing in input .bed \"%s\"" % \
                (seq_id, args.genomic_sites_bed)
            reg_l = reg_len_dic[seq_id]
            assert seq_l == reg_l, \
                "sequence length differs from .bed region length (%i != %i)" \
                % (seq_l, reg_l)
    # Read in model parameters.
    param_dic = gplib.graphprot_get_param_dic(args.in_params)
    # Create GraphProt parameter string.
    param_string = gplib.graphprot_get_param_string(args.in_params)

    Run predictions.

    """
    if args.ws_pred:
        # Do whole site prediction.
        print("Starting whole site predictions on input .fa file"
              " (-action predict) ... ")
        check_cmd = "GraphProt.pl -action predict -prefix " \
            + args.data_id + " -fasta " + args.in_fa + " " \
            + param_string + " -model " + args.in_model
        output = subprocess.getoutput(check_cmd)
        assert output, "the following call of GraphProt.pl "\
            "produced no output:\n%s" % (check_cmd)
        if args.gp_output:
            print(output)
        ws_predictions_file = args.data_id + ".predictions"
        assert os.path.exists(ws_predictions_file), \
            "Whole site prediction output .predictions file \"%s\" not found" \
            % (ws_predictions_file)
        if args.conf_out:
            # Filter by pos_train_ws_pred_median median.
            assert "pos_train_ws_pred_median" in param_dic, \
                "whole site top scores median information "\
                "missing in .params file"
            pos_tr_ws_pred_med = \
                float(param_dic["pos_train_ws_pred_median"])
            # Filtered file.
            filt_ws_predictions_file = args.data_id + ".p50.predictions"
            print("Extracting p50 sites from whole site predictions"
                  " (score threshold = %f) ... " % (pos_tr_ws_pred_med))
            gplib.graphprot_filter_predictions_file(ws_predictions_file,
                                                    filt_ws_predictions_file,
                                                    sc_thr=pos_tr_ws_pred_med)
    else:
        # Do profile prediction.
        print("Starting profile predictions on on input .fa file "
              "(-action predict_profile) ... ")
        check_cmd = "GraphProt.pl -action predict_profile -prefix " \
                    + args.data_id + " -fasta " + args.in_fa + " " + \
                    param_string + " -model " + args.in_model
        output = subprocess.getoutput(check_cmd)
        assert output, "the following call of GraphProt.pl produced "\
                       "no output:\n%s" % (check_cmd)
        if args.gp_output:
            print(output)
        profile_predictions_file = args.data_id + ".profile"
        assert os.path.exists(profile_predictions_file), \
            "Profile prediction output .profile file \"%s\" not found" % \
            (profile_predictions_file)

        # Profile prediction output files.
        avg_prof_file = args.data_id + ".avg_profile"
        avg_prof_peaks_file = args.data_id + ".avg_profile.peaks.bed"
        avg_prof_gen_peaks_file = args.data_id + \
            ".avg_profile.genomic_peaks.bed"
        avg_prof_peaks_p50_file = args.data_id + \
            ".avg_profile.p50.peaks.bed"
        avg_prof_gen_peaks_p50_file = args.data_id + \
            ".avg_profile.p50.genomic_peaks.bed"

        # Get sequence IDs in order from input .fa file.
        seq_ids_list = gplib.fasta_read_in_ids(args.in_fa)
        # Calculate average profiles.
        print("Getting average profile from profile "
              "(extlr for smoothing: %i) ... " % (args.ap_extlr))
        gplib.graphprot_profile_calc_avg_profile(profile_predictions_file,
                                                 avg_prof_file,
                                                 ap_extlr=args.ap_extlr,
                                                 seq_ids_list=seq_ids_list,
                                                 method=2)
        # Extract peak regions on sequences with threshold score 0.
        print("Extracting peak regions from average profile "
              "(score threshold = 0) ... ")
        killpep8 = args.max_merge_dist
        gplib.graphprot_profile_extract_peak_regions(avg_prof_file,
                                                     avg_prof_peaks_file,
                                                     max_merge_dist=killpep8,
                                                     sc_thr=args.score_thr)
        # Convert peaks to genomic coordinates.
        if args.genomic_sites_bed:
            print("Converting peak regions to genomic coordinates ... ")
            killit = args.genomic_sites_bed
            gplib.bed_peaks_to_genomic_peaks(avg_prof_peaks_file,
                                             avg_prof_gen_peaks_file,
                                             print_rows=False,
                                             genomic_sites_bed=killit)
        # Extract peak regions with threshold score p50.
        if args.conf_out:
            sc_id = "pos_train_avg_profile_median_%i" % (args.ap_extlr)
            # Filter by pos_tr_ws_pred_med median.
            assert sc_id in param_dic, "average profile extlr %i median "\
                "information missing in .params file" % (args.ap_extlr)
            p50_sc_thr = float(param_dic[sc_id])
            print("Extracting p50 peak regions from average profile "
                  "(score threshold = %f) ... " % (p50_sc_thr))
            despair = avg_prof_peaks_p50_file
            pain = args.max_merge_dist
            gplib.graphprot_profile_extract_peak_regions(avg_prof_file,
                                                         despair,
                                                         max_merge_dist=pain,
                                                         sc_thr=p50_sc_thr)
            # Convert peaks to genomic coordinates.
            if args.genomic_sites_bed:
                print("Converting p50 peak regions to "
                      "genomic coordinates ... ")
                madness = args.genomic_sites_bed
                gplib.bed_peaks_to_genomic_peaks(avg_prof_peaks_p50_file,
                                                 avg_prof_gen_peaks_p50_file,
                                                 genomic_sites_bed=madness)
    # Done.
    print("Script: I'm done.")
    print("Author: ... ")

5:ddcf35a868b8

import os
import subprocess
import sys

import gplib

"""

TOOL DEPENDENCIES
=================

--sc-thr 0.0 --max-merge-dist 0 --conf-out  --ap-extlr 5

"""


#######################################################################


def setup_argument_parser():
    """Setup argparse parser."""
    help_description = """
    Galaxy wrapper script for GraphProt (-action predict and -action

    If --gen-site-bed .bed file is provided, peak regions will be output
    with genomic coordinates too.

    """
    # Define argument parser.
    p = ap.ArgumentParser(
        add_help=False,
        prog="graphprot_predict_wrapper.py",
        description=help_description,
        formatter_class=ap.MetavarTypeHelpFormatter,
    )

    # Argument groups.
    p_man = p.add_argument_group("REQUIRED ARGUMENTS")
    p_opt = p.add_argument_group("OPTIONAL ARGUMENTS")

    # Required arguments.
    p_opt.add_argument("-h", "--help", action="help", help="Print help message")
    p_man.add_argument(
        "--fasta",
        dest="in_fa",
        type=str,
        required=True,
        help="Sequences .fa file to predict" " on (option -fasta)",
    )
    p_man.add_argument(
        "--model",
        dest="in_model",
        type=str,
        required=True,
        help="GraphProt model file to use for predictions" " (option -model)",
    )
    p_man.add_argument(
        "--params",
        dest="in_params",
        type=str,
        required=True,
        help="Parameter file for given model",
    )
    p_man.add_argument(
        "--data-id",
        dest="data_id",
        type=str,
        required=True,
        help="Data ID (option -prefix)",
    )
    # ---> I'm  a conditional argument <---
    p_opt.add_argument(
        "--ws-pred",
        dest="ws_pred",
        default=False,
        action="store_true",
        help="Run a whole site prediction instead "
        "of calculating profiles (default: false)",
    )
    # Additional arguments.
    p_opt.add_argument(
        "--sc-thr",
        dest="score_thr",
        type=float,
        default=0,
        help="Score threshold for extracting "
        "average profile peak regions (default: 0)",
    )
    p_opt.add_argument(
        "--max-merge-dist",
        dest="max_merge_dist",
        type=int,
        default=0,
        choices=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
        help="Maximum merge distance for nearby peak regions"
        " (default: report all non-overlapping regions)",
    )
    p_opt.add_argument(
        "--gen-site-bed",
        dest="genomic_sites_bed",
        type=str,
        help=".bed file specifying the genomic regions of "
        "the input .fa sequences. Corrupt .bed "
        "information will be punished (default: false)",
    )
    p_opt.add_argument(
        "--conf-out",
        dest="conf_out",
        default=False,
        action="store_true",
        help="Output filtered peak regions BED file or "
        "predictions file (if --ws-pred) using the median "
        "positive training site score for filtering "
        "(default: false)",
    )
    p_opt.add_argument(
        "--gp-output",
        dest="gp_output",
        default=False,
        action="store_true",
        help="Print output produced by GraphProt " "(default: false)",
    )
    p_opt.add_argument(
        "--ap-extlr",
        dest="ap_extlr",
        type=int,
        default=5,
        choices=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
        help="Define average profile up- and downstream "
        "extension to produce the average profile. "
        "The mean over small sequence windows "
        "(window length = --ap-extlr*2 + 1) is used to "
        "get position scores, thus the average profile "
        "is more smooth than the initial profile output "
        "by GraphProt (default: 5)",
    )
    return p


#######################################################################

if __name__ == "__main__":

    # Setup argparse.
    parser = setup_argument_parser()
    # Read in command line arguments.
    args = parser.parse_args()

    # Check for Linux.
    assert "linux" in sys.platform, "please use Linux"
    # Check tool availability.
    assert gplib.is_tool("GraphProt.pl"), "GraphProt.pl not in PATH"
    # Check file inputs.
    assert os.path.exists(args.in_fa), 'input .fa file "%s" not found' % (args.in_fa)
    assert os.path.exists(args.in_model), 'input .model file "%s" not found' % (
        args.in_model
    )
    assert os.path.exists(args.in_params), 'input .params file "%s" not found' % (
        args.in_params
    )
    # Count .fa entries.
    c_in_fa = gplib.count_fasta_headers(args.in_fa)
    assert c_in_fa, 'input .fa file "%s" no headers found' % (args.in_fa)
    print("# input .fa sequences:   %i" % (c_in_fa))
    # Read in FASTA sequences to check for uppercase sequences.
    seqs_dic = gplib.read_fasta_into_dic(args.in_fa)
    # Check for lowercase only sequences, which cause GP to crash.
    error_mess = (
        "input sequences encountered containing "
        "only lowercase characters or lowercase characters in between "
        "uppercase characters. Please provide either all uppercase "
        "sequences or sequences containing uppercase regions surrounded "
        "by lowercase context regions for structure calculation (see "
        "viewpoint concept in original GraphProt publication "
        "for more details)"
    )
    if args.ws_pred:
        bad_ids = gplib.check_seqs_dic_format(seqs_dic)
        assert not bad_ids, "%s" % (error_mess)

    c_uc_nt = gplib.seqs_dic_count_uc_nts(seqs_dic)
    assert c_uc_nt, (
        "no uppercase nucleotides in input .fa sequences. "
        "Please change sequences to uppercase (keep in mind "
        "GraphProt only scores uppercase regions (according "
        "to its viewpoint concept)"
    )
    if not args.ws_pred:
        # Check for lowercase sequences.
        c_lc_nt = gplib.seqs_dic_count_lc_nts(seqs_dic)
        assert not c_lc_nt, (
            "lowercase nucleotides in input .fa not allowed"
            " in profile predictions, since GraphProt only scores "
            "uppercase regions (according to its viewpoint concept)"
        )
    # Check .bed.
    if args.genomic_sites_bed:
        # An array of checks, marvelous.
        assert os.path.exists(
            args.genomic_sites_bed
        ), 'genomic .bed file "%s" not found' % (args.genomic_sites_bed)
        # Check .bed for content.
        assert gplib.count_file_rows(
            args.genomic_sites_bed
        ), 'genomic .bed file "%s" is empty' % (args.genomic_sites_bed)
        # Check .bed for 6-column format.
        assert gplib.bed_check_six_col_format(
            args.genomic_sites_bed
        ), 'genomic .bed file "%s" appears to not ' "be in 6-column .bed format" % (
            args.genomic_sites_bed
        )
        # Check for unique column 4 IDs.
        assert gplib.bed_check_unique_ids(
            args.genomic_sites_bed
        ), 'genomic .bed file "%s" column 4 IDs not unique' % (args.genomic_sites_bed)
        # Read in .bed regions, compare to FASTA sequences.
        seq_len_dic = gplib.get_seq_lengths_from_seqs_dic(seqs_dic)
        reg_len_dic = gplib.bed_get_region_lengths(args.genomic_sites_bed)
        for seq_id in seq_len_dic:
            seq_l = seq_len_dic[seq_id]
            assert (
                seq_id in reg_len_dic
            ), 'sequence ID "%s" missing in input .bed "%s"' % (
                seq_id,
                args.genomic_sites_bed,
            )
            reg_l = reg_len_dic[seq_id]
            assert (
                seq_l == reg_l
            ), "sequence length differs from .bed region length (%i != %i)" % (
                seq_l,
                reg_l,
            )
    # Read in model parameters.
    param_dic = gplib.graphprot_get_param_dic(args.in_params)
    # Create GraphProt parameter string.
    param_string = gplib.graphprot_get_param_string(args.in_params)

    Run predictions.

    """
    if args.ws_pred:
        # Do whole site prediction.
        print(
            "Starting whole site predictions on input .fa file"
            " (-action predict) ... "
        )
        check_cmd = (
            "GraphProt.pl -action predict -prefix "
            + args.data_id
            + " -fasta "
            + args.in_fa
            + " "
            + param_string
            + " -model "
            + args.in_model
        )
        output = subprocess.getoutput(check_cmd)
        assert (
            output
        ), "the following call of GraphProt.pl " "produced no output:\n%s" % (check_cmd)
        if args.gp_output:
            print(output)
        ws_predictions_file = args.data_id + ".predictions"
        assert os.path.exists(
            ws_predictions_file
        ), 'Whole site prediction output .predictions file "%s" not found' % (
            ws_predictions_file
        )
        if args.conf_out:
            # Filter by pos_train_ws_pred_median median.
            assert "pos_train_ws_pred_median" in param_dic, (
                "whole site top scores median information " "missing in .params file"
            )
            pos_tr_ws_pred_med = float(param_dic["pos_train_ws_pred_median"])
            # Filtered file.
            filt_ws_predictions_file = args.data_id + ".p50.predictions"
            print(
                "Extracting p50 sites from whole site predictions"
                " (score threshold = %f) ... " % (pos_tr_ws_pred_med)
            )
            gplib.graphprot_filter_predictions_file(
                ws_predictions_file, filt_ws_predictions_file, sc_thr=pos_tr_ws_pred_med
            )
    else:
        # Do profile prediction.
        print(
            "Starting profile predictions on on input .fa file "
            "(-action predict_profile) ... "
        )
        check_cmd = (
            "GraphProt.pl -action predict_profile -prefix "
            + args.data_id
            + " -fasta "
            + args.in_fa
            + " "
            + param_string
            + " -model "
            + args.in_model
        )
        output = subprocess.getoutput(check_cmd)
        assert (
            output
        ), "the following call of GraphProt.pl produced " "no output:\n%s" % (check_cmd)
        if args.gp_output:
            print(output)
        profile_predictions_file = args.data_id + ".profile"
        assert os.path.exists(
            profile_predictions_file
        ), 'Profile prediction output .profile file "%s" not found' % (
            profile_predictions_file
        )

        # Profile prediction output files.
        avg_prof_file = args.data_id + ".avg_profile"
        avg_prof_peaks_file = args.data_id + ".avg_profile.peaks.bed"
        avg_prof_gen_peaks_file = args.data_id + ".avg_profile.genomic_peaks.bed"
        avg_prof_peaks_p50_file = args.data_id + ".avg_profile.p50.peaks.bed"
        avg_prof_gen_peaks_p50_file = (
            args.data_id + ".avg_profile.p50.genomic_peaks.bed"
        )

        # Get sequence IDs in order from input .fa file.
        seq_ids_list = gplib.fasta_read_in_ids(args.in_fa)
        # Calculate average profiles.
        print(
            "Getting average profile from profile "
            "(extlr for smoothing: %i) ... " % (args.ap_extlr)
        )
        gplib.graphprot_profile_calc_avg_profile(
            profile_predictions_file,
            avg_prof_file,
            ap_extlr=args.ap_extlr,
            seq_ids_list=seq_ids_list,
            method=2,
        )
        # Extract peak regions on sequences with threshold score 0.
        print(
            "Extracting peak regions from average profile " "(score threshold = 0) ... "
        )
        killpep8 = args.max_merge_dist
        gplib.graphprot_profile_extract_peak_regions(
            avg_prof_file,
            avg_prof_peaks_file,
            max_merge_dist=killpep8,
            sc_thr=args.score_thr,
        )
        # Convert peaks to genomic coordinates.
        if args.genomic_sites_bed:
            print("Converting peak regions to genomic coordinates ... ")
            killit = args.genomic_sites_bed
            gplib.bed_peaks_to_genomic_peaks(
                avg_prof_peaks_file,
                avg_prof_gen_peaks_file,
                print_rows=False,
                genomic_sites_bed=killit,
            )
        # Extract peak regions with threshold score p50.
        if args.conf_out:
            sc_id = "pos_train_avg_profile_median_%i" % (args.ap_extlr)
            # Filter by pos_tr_ws_pred_med median.
            assert sc_id in param_dic, (
                "average profile extlr %i median "
                "information missing in .params file" % (args.ap_extlr)
            )
            p50_sc_thr = float(param_dic[sc_id])
            print(
                "Extracting p50 peak regions from average profile "
                "(score threshold = %f) ... " % (p50_sc_thr)
            )
            despair = avg_prof_peaks_p50_file
            pain = args.max_merge_dist
            gplib.graphprot_profile_extract_peak_regions(
                avg_prof_file, despair, max_merge_dist=pain, sc_thr=p50_sc_thr
            )
            # Convert peaks to genomic coordinates.
            if args.genomic_sites_bed:
                print("Converting p50 peak regions to " "genomic coordinates ... ")
                madness = args.genomic_sites_bed
                gplib.bed_peaks_to_genomic_peaks(
                    avg_prof_peaks_p50_file,
                    avg_prof_gen_peaks_p50_file,
                    genomic_sites_bed=madness,
                )
    # Done.
    print("Script: I'm done.")
    print("Author: ... ")