comparison mqppep_mrgfltr.py @ 1:b76c75521d91 draft

planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/mqppep commit 43e7a43b545c24b2dc33d039198551c032aa79be

0:8dfd5d2b5903

        "--phosphopeptides",
        "-p",
        nargs=1,
        required=True,
        dest="phosphopeptides",
        help="Phosphopeptide data for experimental results, including the intensities and the mapping to kinase domains, in tabular format",
    )
    #   UniProtKB/SwissProt DB input, SQLite
    parser.add_argument(
        "--ppep_mapping_db",
        "-d",

        "-x",
        nargs=1,
        required=False,
        default=[],
        dest="species",
        help="limit PhosphoSitePlus records to indicated species (field may be empty)",
    )

    # outputs:
    #   tabular output
    parser.add_argument(

        exit("Error copying ppep_mapping_db to mrgfltr_sqlite: %s" % str(e))

    # determine species to limit records from PSP_Regulatory_Sites
    if options.species is None:
        exit(
            'Argument "species" is required (and may be empty) but not supplied'
        )
    try:
        if len(options.species) > 0:
            species = options.species[0]
        else:

        DEPHOSPHOPEP = "DephosphoPep"
        DESCRIPTION = "Description"
        FUNCTION_PHOSPHORESIDUE = (
            "Function Phosphoresidue(PSP=PhosphoSitePlus.org)"
        )
        GENE_NAME = "Gene_Name"  # Gene Name from UniProtKB
        ON_FUNCTION = (
            "ON_FUNCTION"  # ON_FUNCTION column from PSP_Regulatory_Sites
        )
        ON_NOTES = "NOTES"  # NOTES column from PSP_Regulatory_Sites
        ON_OTHER_INTERACT = "ON_OTHER_INTERACT"  # ON_OTHER_INTERACT column from PSP_Regulatory_Sites
        ON_PROCESS = (
            "ON_PROCESS"  # ON_PROCESS column from PSP_Regulatory_Sites
        )
        ON_PROT_INTERACT = "ON_PROT_INTERACT"  # ON_PROT_INTERACT column from PSP_Regulatory_Sites
        PHOSPHOPEPTIDE = "Phosphopeptide"
        PHOSPHOPEPTIDE_MATCH = "Phosphopeptide_match"
        PHOSPHORESIDUE = "Phosphoresidue"
        PUTATIVE_UPSTREAM_DOMAINS = "Putative Upstream Kinases(PSP=PhosphoSitePlus.org)/Phosphatases/Binding Domains"
        SEQUENCE = "Sequence"
        SEQUENCE10 = "Sequence10"
        SEQUENCE7 = "Sequence7"
        SITE_PLUSMINUS_7AA_SQL = "SITE_PLUSMINUS_7AA"
        UNIPROT_ID = "UniProt_ID"

          INSERT INTO Citation (
            ObjectName,
            CitationData
          ) VALUES (?,?)
          """
        CITATION_INSERT_PSP = 'PhosphoSitePlus(R) (PSP) was created by Cell Signaling Technology Inc. It is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License. When using PSP data or analyses in printed publications or in online resources, the following acknowledgements must be included: (a) the words "PhosphoSitePlus(R), www.phosphosite.org" must be included at appropriate places in the text or webpage, and (b) the following citation must be included in the bibliography: "Hornbeck PV, Zhang B, Murray B, Kornhauser JM, Latham V, Skrzypek E PhosphoSitePlus, 2014: mutations, PTMs and recalibrations. Nucleic Acids Res. 2015 43:D512-20. PMID: 25514926."'
        CITATION_INSERT_PSP_REF = 'Hornbeck, 2014, "PhosphoSitePlus, 2014: mutations, PTMs and recalibrations.", https://pubmed.ncbi.nlm.nih.gov/22135298, https://doi.org/10.1093/nar/gkr1122'

        MRGFLTR_METADATA_COLUMNS = [
            "ppep_id",
            "Sequence10",
            "Sequence7",

                m = re_tab.match(line)
                upstream_map_p_peptide_list.append(m[0])
        file1.close()
        file1_encoded.close()

        # Get the list of phosphopeptides with the p's that represent the phosphorylation sites removed
        re_phos = re.compile("p")

        end_time = time.process_time()  # timer
        print(
            "%0.6f pre-read-SwissProt [0.1]" % (end_time - start_time,),
            file=sys.stderr,
        )

        # ----------- Get SwissProt data from SQLite database (start) -----------
        # build UniProt sequence LUT and list of unique SwissProt sequences

        # Open SwissProt SQLite database
        conn = sql.connect(uniprot_sqlite)
        cur = conn.cursor()

        Sequence = ""
        UniProt_ID = ""
        Description = ""
        Gene_Name = ""

        # ----------- Get SwissProt data from SQLite database (finish) -----------

        end_time = time.process_time()  # timer
        print(
            "%0.6f post-read-SwissProt [0.2]" % (end_time - start_time,),
            file=sys.stderr,
        )

        # ----------- Get SwissProt data from SQLite database (start) -----------
        # Open SwissProt SQLite database
        conn = sql.connect(uniprot_sqlite)
        cur = conn.cursor()

        # Set up dictionary to aggregate results for phosphopeptides correspounding to dephosphoeptide
        DephosphoPep_UniProtSeq_LUT = {}

        # Set up dictionary to accumulate results
        PhosphoPep_UniProtSeq_LUT = {}

                    #   ref: https://stackoverflow.com/a/4174955/15509512
                    r = sorted(r, key=operator.itemgetter(0))
                    # Get one tuple for each `phospho_pep`
                    #   in DephosphoPep_UniProtSeq_LUT[dephospho_pep]
                    for (upid, gn, desc) in r:
                        # Append pseudo-tuple per UniProt_ID but only when it is not present
                        if (
                            upid
                            not in DephosphoPep_UniProtSeq_LUT[
                                (dephospho_pep, UNIPROT_ID)
                            ]

        conn.close()
        # wipe local variables
        phospho_pep = dephospho_pep = sequence = 0
        upid = gn = desc = r = ""

        # ----------- Get SwissProt data from SQLite database (finish) -----------

        end_time = time.process_time()  # timer
        print(
            "%0.6f finished reading and decoding '%s' [0.4]"
            % (end_time - start_time, upstream_map_filename_tab),

            "%0.6f added index to Upstream Map [1g_2]"
            % (end_time - start_time,),
            file=sys.stderr,
        )  # timer

        # ########################################################################
        # # trim upstream_data to include only the upstream map columns
        # old_cols = upstream_data.columns.tolist()
        # i = 0
        # first_intensity = -1
        # last_intensity = -1
        # intensity_re = re.compile("Intensity.*")
        # for col_name in old_cols:
        #     m = intensity_re.match(col_name)
        #     if m:
        #         last_intensity = i
        #         if first_intensity == -1:
        #             first_intensity = i
        #     i += 1
        # # print('last intensity = %d' % last_intensity)
        # col_PKCalpha = last_intensity + 2
        #
        # data_in_cols = [old_cols[0]] + old_cols[
        #     first_intensity: last_intensity + 1
        # ]
        #
        # if upstream_data.empty:
        #     print("upstream_data is empty")
        #     exit(0)
        #
        # data_in = upstream_data.copy(deep=True)[data_in_cols]
        ########################################################################
        # trim upstream_data to include only the upstream map columns
        old_cols = upstream_data.columns.tolist()
        i = 0
        first_intensity = -1
        last_intensity = -1

        data_in = upstream_data.copy(deep=True)[data_in_cols]
        data_in.columns = data_col_names
        print("data_in")
        print(data_in)
        ########################################################################

        # Convert floating-point integers to int64 integers
        #   ref: https://stackoverflow.com/a/68497603/15509512
        data_in[list(data_in.columns[1:])] = (
            data_in[list(data_in.columns[1:])]

            % (end_time - start_time,),
            file=sys.stderr,
        )  # timer

        # Produce a dictionary of metadata for a single phosphopeptide.
        #   This is a replacement of `UniProtInfo_subdict` in the original code.
        def pseq_to_subdict(phospho_pep):
            # Strip "p" from phosphopeptide sequence
            dephospho_pep = re_phos.sub("", phospho_pep)

            # Determine number of phosphoresidues in phosphopeptide

            #             )
            #         )
            if phospho_pep not in PhosphoPep_UniProtSeq_LUT:
                raise PreconditionError(
                    phospho_pep,
                    "no matching phosphopeptide found in PhosphoPep_UniProtSeq_LUT",
                )
            if dephospho_pep not in DephosphoPep_UniProtSeq_LUT:
                raise PreconditionError(
                    dephospho_pep,
                    "dephosphorylated phosphopeptide not found in DephosphoPep_UniProtSeq_LUT",
                )
            if (
                dephospho_pep != PhosphoPep_UniProtSeq_LUT[(phospho_pep, DEPHOSPHOPEP)]
            ):
                my_err_msg = "dephosphorylated phosphopeptide does not match "
                my_err_msg += "PhosphoPep_UniProtSeq_LUT[(phospho_pep,DEPHOSPHOPEP)] = "
                my_err_msg += PhosphoPep_UniProtSeq_LUT[(phospho_pep, DEPHOSPHOPEP)]
                raise PreconditionError(dephospho_pep, my_err_msg)

            result[SEQUENCE] = [dephospho_pep]
            result[UNIPROT_ID] = DephosphoPep_UniProtSeq_LUT[
                (dephospho_pep, UNIPROT_ID)

                (dephospho_pep, GENE_NAME)
            ]
            if (dephospho_pep, SEQUENCE) not in DephosphoPep_UniProtSeq_LUT:
                raise PreconditionError(
                    dephospho_pep,
                    "no matching phosphopeptide found in DephosphoPep_UniProtSeq_LUT",
                )
            UniProtSeqList = DephosphoPep_UniProtSeq_LUT[
                (dephospho_pep, SEQUENCE)
            ]
            if len(UniProtSeqList) < 1:
                print(
                    "Skipping DephosphoPep_UniProtSeq_LUT[('%s',SEQUENCE)] because value has zero length"
                    % dephospho_pep
                )
                # raise PreconditionError(
                #     "DephosphoPep_UniProtSeq_LUT[('" + dephospho_pep + ",SEQUENCE)",
                #      'value has zero length'
                #      )
            for UniProtSeq in UniProtSeqList:
                i = 0
                phosphoresidues = []
                seq7s_set = set()
                seq7s = []

                        seq7s_set.add(seq7)

                    # add Sequence10
                    if psite < 10:  # phospho_pep at N terminus
                        seq10 = (
                            str(UniProtSeq)[:psite] + "p" + str(UniProtSeq)[psite: psite + 11]
                        )
                    elif (
                        len(UniProtSeq) - psite < 11
                    ):  # phospho_pep at C terminus
                        seq10 = (
                            str(UniProtSeq)[psite - 10: psite] + "p" + str(UniProtSeq)[psite:]
                        )
                    else:
                        seq10 = str(UniProtSeq)[psite - 10: psite + 11]
                        seq10 = seq10[:10] + "p" + seq10[10:]
                    if seq10 not in seq10s_set:

                    result[key] = joined_value

            newstring = "; ".join(
                [", ".join(prez) for prez in result[PHOSPHORESIDUE]]
            )
            # #separate the isoforms in PHOSPHORESIDUE column with ";"
            # oldstring = result[PHOSPHORESIDUE]
            # oldlist = list(oldstring)
            # newstring = ""
            # i = 0
            # for e in oldlist:
            #     if e == ";":
            #         if numps > 1:
            #             if i%numps:
            #                 newstring = newstring + ";"
            #             else:
            #                 newstring = newstring + ","
            #         else:
            #             newstring = newstring + ";"
            #         i +=1
            #     else:
            #         newstring = newstring + e
            result[PHOSPHORESIDUE] = newstring

            # separate sequence7's by |
            oldstring = result[SEQUENCE7]
            oldlist = oldstring

            % (end_time - start_time,),
            file=sys.stderr,
        )  # timer

        # Construct dictionary from list of lists
        #   ref: https://www.8bitavenue.com/how-to-convert-list-of-lists-to-dictionary-in-python/
        UniProt_Info = {
            result[0]: result[1]
            for result in result_list
            if result is not None
        }

        end_time = time.process_time()  # timer
        print(
            "%0.6f create dictionary mapping phosphopeptide to metadata dictionary [C]"
            % (end_time - start_time,),
            file=sys.stderr,
        )  # timer

        # cosmetic: add N_A to phosphopeptide rows with no hits
        p_peptide_list = []

                for _, row in uniprot_df.iterrows()
            ]
        ).reset_index()
        seq7_df.columns = [SEQUENCE7, PHOSPHOPEPTIDE]

        # --- -------------- begin read PSP_Regulatory_sites ---------------------------------
        # read in PhosphoSitePlus Regulatory Sites dataset
        # ----------- Get PhosphoSitePlus Regulatory Sites data from SQLite database (start) -----------
        conn = sql.connect(uniprot_sqlite)
        regsites_df = pandas.read_sql_query(PSP_REGSITE_SQL, conn)
        # Close SwissProt SQLite database
        conn.close()
        # ... -------------- end read PSP_Regulatory_sites ------------------------------------

        # keep only the human entries in dataframe
        if len(species) > 0:
            print(
                'Limit PhosphoSitesPlus records to species "' + species + '"'

                fp_series[i] = fp_series[i][2:]
            fp_series[i] = fp_series[i].replace("; ; ; ; ", "; ")
            fp_series[i] = fp_series[i].replace("; ; ; ", "; ")
            fp_series[i] = fp_series[i].replace("; ; ", "; ")

            # turn blanks into N_A to signify the info was searched for but cannot be found
            if fp_series[i] == "":
                fp_series[i] = N_A

            i += 1
        merge_df[FUNCTION_PHOSPHORESIDUE] = fp_series

            {"Protein description": PHOSPHOPEPTIDE},
            axis="columns",
            inplace=True,
        )

        # data_in.sort_values(PHOSPHOPEPTIDE_MATCH, inplace=True, kind='mergesort')
        res2 = sorted(
            data_in[PHOSPHOPEPTIDE_MATCH].tolist(), key=lambda s: s.casefold()
        )
        data_in = data_in.loc[res2]

        print(
            "%0.6f refactored columns for Upstream Map [1g]"
            % (end_time - start_time,),
            file=sys.stderr,
        )  # timer

        # #rename upstream columns in new list
        # new_cols = []
        # for name in cols:
        #     if "_NetworKIN" in name:
        #         name = name.split("_")[0]
        #     if " motif" in name:
        #         name = name.split(" motif")[0]
        #     if " sequence " in name:
        #         name = name.split(" sequence")[0]
        #     if "_Phosida" in name:
        #         name = name.split("_")[0]
        #     if "_PhosphoSite" in name:
        #         name = name.split("_")[0]
        #     new_cols.append(name)

        # rename upstream columns in new list
        def col_rename(name):
            if "_NetworKIN" in name:
                name = name.split("_")[0]

        new_cols_cast = new_cols
        new_cols_cast[0] = "p_peptide"
        upstream_data_cast.columns = new_cols_cast
        upstream_data_cast["p_peptide"] = upstream_data.index

        # --- -------------- begin read upstream_data_melt ------------------------------------
        # ----------- Get melted kinase mapping data from SQLite database (start) -----------
        conn = sql.connect(uniprot_sqlite)
        upstream_data_melt_df = pandas.read_sql_query(PPEP_MELT_SQL, conn)
        # Close SwissProt SQLite database
        conn.close()
        upstream_data_melt = upstream_data_melt_df.copy()

        print(
            "%0.6f upstream_data_melt_df pre-dedup has %d rows"
            % (end_time - start_time, len(upstream_data_melt.axes[0])),
            file=sys.stderr,
        )
        # ----------- Get melted kinase mapping data from SQLite database (finish) -----------
        # ... -------------- end read upstream_data_melt --------------------------------------

        end_time = time.process_time()  # timer
        print(
            "%0.6f melted and minimized Upstream Map dataframe [1h_2]"
            % (end_time - start_time,),

        for el in upstream_delta_melt_LoL:
            (p_peptide, characterization, X) = tuple(el)
            if p_peptide in melt_dict:
                melt_dict[p_peptide].append(characterization)
            else:
                exit(
                    'Phosphopeptide %s not found in ppep_mapping_db: "phopsphopeptides" and "ppep_mapping_db" must both originate from the same run of mqppep_kinase_mapping'
                    % (p_peptide)
                )

        end_time = time.process_time()  # timer
        print(
            "%0.6f appended peptide characterizations [1h_2b]"
            % (end_time - start_time,),

                FUNCTION_PHOSPHORESIDUE,
                PUTATIVE_UPSTREAM_DOMAINS,
            ]
        ]

        # cols_output_prelim = output_df.columns.tolist()
        #
        # print("cols_output_prelim")
        # print(cols_output_prelim)
        #
        # cols_output = cols_output_prelim[:8]+[cols_output_prelim[9]]+[cols_output_prelim[10]]
        #
        # print("cols_output with p-peptide")
        # print(cols_output)
        #
        # cols_output = [col for col in cols_output if not col == "p-peptide"]
        #
        # print("cols_output")
        # print(cols_output)
        #
        # output_df = output_df[cols_output]

        # join output_df back to quantitative columns in data_in df
        quant_cols = data_in.columns.tolist()
        quant_cols = quant_cols[1:]
        quant_data = data_in[quant_cols]

        # ----------- Write merge/filter metadata to SQLite database (start) -----------
        # Open SwissProt SQLite database
        conn = sql.connect(output_sqlite)
        cur = conn.cursor()

        cur.executescript(MRGFLTR_DDL)

            method="multi",
        )

        # Close SwissProt SQLite database
        conn.close()
        # ----------- Write merge/filter metadata to SQLite database (finish) -----------

        output_df = output_df.merge(
            quant_data,
            how="right",
            left_on=PHOSPHOPEPTIDE,

        output_cols = output_df.columns.tolist()
        output_cols = output_cols[:-1]
        output_df = output_df[output_cols]

        # cosmetic changes to Upstream column
        output_df[PUTATIVE_UPSTREAM_DOMAINS] = output_df[
            PUTATIVE_UPSTREAM_DOMAINS
        ].fillna(
            ""
        )  # fill the NaN with "" for those Phosphopeptides that got a "WARNING: Failed match for " in the upstream mapping
        us_series = pandas.Series(output_df[PUTATIVE_UPSTREAM_DOMAINS])
        i = 0
        while i < len(us_series):
            # turn blanks into N_A to signify the info was searched for but cannot be found
            if us_series[i] == "":
                us_series[i] = N_A
            i += 1
        output_df[PUTATIVE_UPSTREAM_DOMAINS] = us_series

        )  # timer

        # Rev. 7/1/2016
        # Rev. 7/3/2016 : fill NaN in Upstream column to replace to N/A's
        # Rev. 7/3/2016:  renamed Upstream column to PUTATIVE_UPSTREAM_DOMAINS
        # Rev. 12/2/2021: Converted to Python from ipynb; use fast Aho-Corasick searching; \
        #                read from SwissProt SQLite database
        # Rev. 12/9/2021: Transfer code to Galaxy tool wrapper

        #
        # copied from Excel Output Script.ipynb END #
        #

1:b76c75521d91

        "--phosphopeptides",
        "-p",
        nargs=1,
        required=True,
        dest="phosphopeptides",
        help=" ".join([
            "Phosphopeptide data for experimental results, including the",
            "intensities and the mapping to kinase domains, in tabular format"
        ]),
    )
    #   UniProtKB/SwissProt DB input, SQLite
    parser.add_argument(
        "--ppep_mapping_db",
        "-d",

        "-x",
        nargs=1,
        required=False,
        default=[],
        dest="species",
        help=" ".join([
            "limit PhosphoSitePlus records to indicated species",
            "(field may be empty)"
        ]),
    )

    # outputs:
    #   tabular output
    parser.add_argument(

        exit("Error copying ppep_mapping_db to mrgfltr_sqlite: %s" % str(e))

    # determine species to limit records from PSP_Regulatory_Sites
    if options.species is None:
        exit(
            'Argument "species" is required (& may be empty) but not supplied'
        )
    try:
        if len(options.species) > 0:
            species = options.species[0]
        else:

        DEPHOSPHOPEP = "DephosphoPep"
        DESCRIPTION = "Description"
        FUNCTION_PHOSPHORESIDUE = (
            "Function Phosphoresidue(PSP=PhosphoSitePlus.org)"
        )
        # Gene Name from UniProtKB
        GENE_NAME = "Gene_Name"
        # ON_FUNCTION column from PSP_Regulatory_Sites
        ON_FUNCTION = ("ON_FUNCTION")
        # NOTES column from PSP_Regulatory_Sites
        ON_NOTES = "NOTES"
        # ON_OTHER_INTERACT column from PSP_Regulatory_Sites
        ON_OTHER_INTERACT = "ON_OTHER_INTERACT"
        # ON_PROCESS column from PSP_Regulatory_Sites
        ON_PROCESS = ("ON_PROCESS")
        # ON_PROT_INTERACT column from PSP_Regulatory_Sites
        ON_PROT_INTERACT = "ON_PROT_INTERACT"
        PHOSPHOPEPTIDE = "Phosphopeptide"
        PHOSPHOPEPTIDE_MATCH = "Phosphopeptide_match"
        PHOSPHORESIDUE = "Phosphoresidue"
        PUTATIVE_UPSTREAM_DOMAINS = " ".join([
            "Putative Upstream Kinases(PSP=PhosphoSitePlus.org)/",
            "Phosphatases/Binding Domains"
        ])
        SEQUENCE = "Sequence"
        SEQUENCE10 = "Sequence10"
        SEQUENCE7 = "Sequence7"
        SITE_PLUSMINUS_7AA_SQL = "SITE_PLUSMINUS_7AA"
        UNIPROT_ID = "UniProt_ID"

          INSERT INTO Citation (
            ObjectName,
            CitationData
          ) VALUES (?,?)
          """
        CITATION_INSERT_PSP = " ".join([
            "PhosphoSitePlus(R) (PSP) was created by Cell Signaling",
            "Technology Inc. It is licensed under a Creative Commons",
            "Attribution-NonCommercial-ShareAlike 3.0 Unported License.",
            "When using PSP data or analyses in printed publications or",
            "in online resources, the following acknowledgements must be",
            "included: (a) the words",
            '"PhosphoSitePlus(R), www.phosphosite.org" must',
            "be included at appropriate places in the text or webpage,",
            "and (b) the following citation must be included in the",
            'bibliography: "Hornbeck PV, Zhang B, Murray B, Kornhauser',
            "JM, Latham V, Skrzypek E PhosphoSitePlus, 2014: mutations,",
            "PTMs and recalibrations. Nucleic Acids Res. 2015",
            '43:D512-20. PMID: 25514926."'
        ])
        CITATION_INSERT_PSP_REF = " ".join([
            'Hornbeck, 2014, "PhosphoSitePlus, 2014: mutations, PTMs and',
            'recalibrations.", https://pubmed.ncbi.nlm.nih.gov/22135298,',
            "https://doi.org/10.1093/nar/gkr1122"
        ])

        MRGFLTR_METADATA_COLUMNS = [
            "ppep_id",
            "Sequence10",
            "Sequence7",

                m = re_tab.match(line)
                upstream_map_p_peptide_list.append(m[0])
        file1.close()
        file1_encoded.close()

        # Get the list of phosphopeptides with the p's that represent
        #     the phosphorylation sites removed
        re_phos = re.compile("p")

        end_time = time.process_time()  # timer
        print(
            "%0.6f pre-read-SwissProt [0.1]" % (end_time - start_time,),
            file=sys.stderr,
        )

        # -------- Get SwissProt data from SQLite database (start) -----------
        # build UniProt sequence LUT and list of unique SwissProt sequences

        # Open SwissProt SQLite database
        conn = sql.connect(uniprot_sqlite)
        cur = conn.cursor()

        Sequence = ""
        UniProt_ID = ""
        Description = ""
        Gene_Name = ""

        # -------- Get SwissProt data from SQLite database (finish) -----------

        end_time = time.process_time()  # timer
        print(
            "%0.6f post-read-SwissProt [0.2]" % (end_time - start_time,),
            file=sys.stderr,
        )

        # -------- Get SwissProt data from SQLite database (start) -----------
        # Open SwissProt SQLite database
        conn = sql.connect(uniprot_sqlite)
        cur = conn.cursor()

        # Set up dictionary to aggregate results for phosphopeptides
        #     corresponding to dephosphoeptide
        DephosphoPep_UniProtSeq_LUT = {}

        # Set up dictionary to accumulate results
        PhosphoPep_UniProtSeq_LUT = {}

                    #   ref: https://stackoverflow.com/a/4174955/15509512
                    r = sorted(r, key=operator.itemgetter(0))
                    # Get one tuple for each `phospho_pep`
                    #   in DephosphoPep_UniProtSeq_LUT[dephospho_pep]
                    for (upid, gn, desc) in r:
                        # Append pseudo-tuple per UniProt_ID iff not present
                        if (
                            upid
                            not in DephosphoPep_UniProtSeq_LUT[
                                (dephospho_pep, UNIPROT_ID)
                            ]

        conn.close()
        # wipe local variables
        phospho_pep = dephospho_pep = sequence = 0
        upid = gn = desc = r = ""

        # -------- Get SwissProt data from SQLite database (finish) -----------

        end_time = time.process_time()  # timer
        print(
            "%0.6f finished reading and decoding '%s' [0.4]"
            % (end_time - start_time, upstream_map_filename_tab),

            "%0.6f added index to Upstream Map [1g_2]"
            % (end_time - start_time,),
            file=sys.stderr,
        )  # timer

        # trim upstream_data to include only the upstream map columns
        old_cols = upstream_data.columns.tolist()
        i = 0
        first_intensity = -1
        last_intensity = -1

        data_in = upstream_data.copy(deep=True)[data_in_cols]
        data_in.columns = data_col_names
        print("data_in")
        print(data_in)
        ######################################################################

        # Convert floating-point integers to int64 integers
        #   ref: https://stackoverflow.com/a/68497603/15509512
        data_in[list(data_in.columns[1:])] = (
            data_in[list(data_in.columns[1:])]

            % (end_time - start_time,),
            file=sys.stderr,
        )  # timer

        # Produce a dictionary of metadata for a single phosphopeptide.
        #   This is replaces `UniProtInfo_subdict` from the original code.
        def pseq_to_subdict(phospho_pep):
            # Strip "p" from phosphopeptide sequence
            dephospho_pep = re_phos.sub("", phospho_pep)

            # Determine number of phosphoresidues in phosphopeptide

            #             )
            #         )
            if phospho_pep not in PhosphoPep_UniProtSeq_LUT:
                raise PreconditionError(
                    phospho_pep,
                    " ".join([
                        "no matching phosphopeptide found in",
                        "PhosphoPep_UniProtSeq_LUT"
                    ])
                )
            if dephospho_pep not in DephosphoPep_UniProtSeq_LUT:
                raise PreconditionError(
                    dephospho_pep,
                    " ".join([
                        "dephosphorylated phosphopeptide not found",
                        "in DephosphoPep_UniProtSeq_LUT"
                    ])
                )
            if (
                dephospho_pep != PhosphoPep_UniProtSeq_LUT[
                    (phospho_pep, DEPHOSPHOPEP)]
            ):
                my_err_msg = "dephosphorylated phosphopeptide does not match "
                my_err_msg += "PhosphoPep_UniProtSeq_LUT"
                my_err_msg += "[(phospho_pep,DEPHOSPHOPEP)] = "
                my_err_msg += PhosphoPep_UniProtSeq_LUT[
                    (phospho_pep, DEPHOSPHOPEP)]
                raise PreconditionError(dephospho_pep, my_err_msg)

            result[SEQUENCE] = [dephospho_pep]
            result[UNIPROT_ID] = DephosphoPep_UniProtSeq_LUT[
                (dephospho_pep, UNIPROT_ID)

                (dephospho_pep, GENE_NAME)
            ]
            if (dephospho_pep, SEQUENCE) not in DephosphoPep_UniProtSeq_LUT:
                raise PreconditionError(
                    dephospho_pep,
                    "".join(
                        "no matching phosphopeptide found",
                        "in DephosphoPep_UniProtSeq_LUT")
                )
            UniProtSeqList = DephosphoPep_UniProtSeq_LUT[
                (dephospho_pep, SEQUENCE)
            ]
            if len(UniProtSeqList) < 1:
                print(
                    "Skipping DephosphoPep_UniProtSeq_LUT[('%s',SEQUENCE)] %s"
                    % (dephospho_pep, "because value has zero length")
                )
            for UniProtSeq in UniProtSeqList:
                i = 0
                phosphoresidues = []
                seq7s_set = set()
                seq7s = []

                        seq7s_set.add(seq7)

                    # add Sequence10
                    if psite < 10:  # phospho_pep at N terminus
                        seq10 = (
                            str(UniProtSeq)[:psite] + "p" + str(
                                UniProtSeq)[psite: psite + 11]
                        )
                    elif (
                        len(UniProtSeq) - psite < 11
                    ):  # phospho_pep at C terminus
                        seq10 = (
                            str(UniProtSeq)[psite - 10: psite] + "p" + str(
                                UniProtSeq)[psite:]
                        )
                    else:
                        seq10 = str(UniProtSeq)[psite - 10: psite + 11]
                        seq10 = seq10[:10] + "p" + seq10[10:]
                    if seq10 not in seq10s_set:

                    result[key] = joined_value

            newstring = "; ".join(
                [", ".join(prez) for prez in result[PHOSPHORESIDUE]]
            )
            result[PHOSPHORESIDUE] = newstring

            # separate sequence7's by |
            oldstring = result[SEQUENCE7]
            oldlist = oldstring

            % (end_time - start_time,),
            file=sys.stderr,
        )  # timer

        # Construct dictionary from list of lists
        #   ref: https://www.8bitavenue.com/\
        #        how-to-convert-list-of-lists-to-dictionary-in-python/
        UniProt_Info = {
            result[0]: result[1]
            for result in result_list
            if result is not None
        }

        end_time = time.process_time()  # timer
        print(
            "%0.6f create dictionary mapping phosphopeptide %s"
            % (end_time - start_time, "to metadata dictionary [C]"),
            file=sys.stderr,
        )  # timer

        # cosmetic: add N_A to phosphopeptide rows with no hits
        p_peptide_list = []

                for _, row in uniprot_df.iterrows()
            ]
        ).reset_index()
        seq7_df.columns = [SEQUENCE7, PHOSPHOPEPTIDE]

        # read in PhosphoSitePlus Regulatory Sites dataset from SQLite
        # --- -------------- begin read PSP_Regulatory_sites -----
        conn = sql.connect(uniprot_sqlite)
        regsites_df = pandas.read_sql_query(PSP_REGSITE_SQL, conn)
        # Close SwissProt SQLite database
        conn.close()
        # ... -------------- end read PSP_Regulatory_sites -------

        # keep only the human entries in dataframe
        if len(species) > 0:
            print(
                'Limit PhosphoSitesPlus records to species "' + species + '"'

                fp_series[i] = fp_series[i][2:]
            fp_series[i] = fp_series[i].replace("; ; ; ; ", "; ")
            fp_series[i] = fp_series[i].replace("; ; ; ", "; ")
            fp_series[i] = fp_series[i].replace("; ; ", "; ")

            # turn blanks into N_A to signify the info
            #     that was searched for but cannot be found
            if fp_series[i] == "":
                fp_series[i] = N_A

            i += 1
        merge_df[FUNCTION_PHOSPHORESIDUE] = fp_series

            {"Protein description": PHOSPHOPEPTIDE},
            axis="columns",
            inplace=True,
        )

        res2 = sorted(
            data_in[PHOSPHOPEPTIDE_MATCH].tolist(), key=lambda s: s.casefold()
        )
        data_in = data_in.loc[res2]

        print(
            "%0.6f refactored columns for Upstream Map [1g]"
            % (end_time - start_time,),
            file=sys.stderr,
        )  # timer

        # rename upstream columns in new list
        def col_rename(name):
            if "_NetworKIN" in name:
                name = name.split("_")[0]

        new_cols_cast = new_cols
        new_cols_cast[0] = "p_peptide"
        upstream_data_cast.columns = new_cols_cast
        upstream_data_cast["p_peptide"] = upstream_data.index

        # Get melted kinase mapping data from SQLite database
        # --- begin read upstream_data_melt ---------------------------------
        conn = sql.connect(uniprot_sqlite)
        upstream_data_melt_df = pandas.read_sql_query(PPEP_MELT_SQL, conn)
        # Close SwissProt SQLite database
        conn.close()
        upstream_data_melt = upstream_data_melt_df.copy()

        print(
            "%0.6f upstream_data_melt_df pre-dedup has %d rows"
            % (end_time - start_time, len(upstream_data_melt.axes[0])),
            file=sys.stderr,
        )
        # ... end read upstream_data_melt ---------------------------------

        end_time = time.process_time()  # timer
        print(
            "%0.6f melted and minimized Upstream Map dataframe [1h_2]"
            % (end_time - start_time,),

        for el in upstream_delta_melt_LoL:
            (p_peptide, characterization, X) = tuple(el)
            if p_peptide in melt_dict:
                melt_dict[p_peptide].append(characterization)
            else:
                los = [
                    "Phosphopeptide %s" % p_peptide,
                    "not found in ppep_mapping_db:",
                    '"phopsphopeptides" and "ppep_mapping_db" must both',
                    "originate from the same run of mqppep_kinase_mapping"
                ]
                exit(" ".join(los))

        end_time = time.process_time()  # timer
        print(
            "%0.6f appended peptide characterizations [1h_2b]"
            % (end_time - start_time,),

                FUNCTION_PHOSPHORESIDUE,
                PUTATIVE_UPSTREAM_DOMAINS,
            ]
        ]

        # join output_df back to quantitative columns in data_in df
        quant_cols = data_in.columns.tolist()
        quant_cols = quant_cols[1:]
        quant_data = data_in[quant_cols]

        # ---- Write merge/filter metadata to SQLite database (start) ----
        # Open SwissProt SQLite database
        conn = sql.connect(output_sqlite)
        cur = conn.cursor()

        cur.executescript(MRGFLTR_DDL)

            method="multi",
        )

        # Close SwissProt SQLite database
        conn.close()
        # ---- Write merge/filter metadata to SQLite database (finish) ----

        output_df = output_df.merge(
            quant_data,
            how="right",
            left_on=PHOSPHOPEPTIDE,

        output_cols = output_df.columns.tolist()
        output_cols = output_cols[:-1]
        output_df = output_df[output_cols]

        # cosmetic changes to Upstream column
        # fill the NaN with "" for those Phosphopeptides that got a
        #   "WARNING: Failed match for " in the upstream mapping
        output_df[PUTATIVE_UPSTREAM_DOMAINS] = output_df[
            PUTATIVE_UPSTREAM_DOMAINS
        ].fillna(
            ""
        )
        us_series = pandas.Series(output_df[PUTATIVE_UPSTREAM_DOMAINS])
        i = 0
        while i < len(us_series):
            # turn blanks into N_A to signify the info
            #   that was searched for but cannot be found
            if us_series[i] == "":
                us_series[i] = N_A
            i += 1
        output_df[PUTATIVE_UPSTREAM_DOMAINS] = us_series

        )  # timer

        # Rev. 7/1/2016
        # Rev. 7/3/2016 : fill NaN in Upstream column to replace to N/A's
        # Rev. 7/3/2016:  renamed Upstream column to PUTATIVE_UPSTREAM_DOMAINS
        # Rev. 12/2/2021: Converted to Python from ipynb; use fast \
        #                 Aho-Corasick searching; \
        #                 read from SwissProt SQLite database
        # Rev. 12/9/2021: Transfer code to Galaxy tool wrapper

        #
        # copied from Excel Output Script.ipynb END #
        #