test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys, re, numpy as np, pandas as pd\n",
    "from pyteomics import mzxml, mzml\n",
    "\n",
    "\n",
    "class progressBar:\n",
    "    def __init__(self, total):\n",
    "        self.total = total\n",
    "        self.barLength = 20\n",
    "        self.count = 0\n",
    "        self.progress = 0\n",
    "        self.block = 0\n",
    "        self.status = \"\"\n",
    "\n",
    "    def increment(self):\n",
    "        self.count += 1\n",
    "        self.progress = self.count / self.total\n",
    "        self.block = int(round(self.barLength * self.progress))\n",
    "        if self.progress == 1:\n",
    "            self.status = \"Done...\\r\\n\"\n",
    "        else:\n",
    "            self.status = \"\"\n",
    "        #         self.status = str(self.count) + \"/\" + str(self.total)\n",
    "        text = \"\\r  Progress: [{0}] {1}% {2}\".format(\"#\" * self.block + \"-\" * (self.barLength - self.block),\n",
    "                                                     int(self.progress * 100), self.status)\n",
    "        sys.stdout.write(text)\n",
    "        sys.stdout.flush()\n",
    "\n",
    "\n",
    "def getParams(paramFile):\n",
    "    parameters = dict()\n",
    "    with open(paramFile, 'r') as file:\n",
    "        for line in file:\n",
    "            if re.search(r'^#', line) or re.search(r'^\\s', line):\n",
    "                continue\n",
    "            line = re.sub(r'#.*', '', line)  # Remove comments (start from '#')\n",
    "            line = re.sub(r'\\s*', '', line)  # Remove all whitespaces\n",
    "\n",
    "            # Exception for \"feature_files\" parameter\n",
    "            if \"feature_files\" in parameters and line.endswith(\"feature\"):\n",
    "                parameters[\"feature_files\"].append(line)\n",
    "            else:\n",
    "                key = line.split('=')[0]\n",
    "                val = line.split('=')[1]\n",
    "                if key == \"feature_files\":\n",
    "                    parameters[key] = [val]\n",
    "                else:\n",
    "                    parameters[key] = val\n",
    "    return parameters\n",
    "\n",
    "\n",
    "def parseIdtxt(idTxt, params):\n",
    "    df = pd.read_table(idTxt, sep=\";\", skiprows=0, header=1)\n",
    "    df = df.groupby([\"Peptide\", \"Protein\", \"Outfile\"]).size().reset_index()\n",
    "#     outfiles = sorted(list(set(df[\"Outfile\"]))) # Unique PSMs (i.e. MS2 scan numbers)\n",
    "    psms = {}\n",
    "    pep2psm = {}\n",
    "    prot2psm = {}\n",
    "    for i in range(df.shape[0]):\n",
    "        pep = df[\"Peptide\"][i]\n",
    "        prot = df[\"Protein\"][i]\n",
    "        psm = df[\"Outfile\"][i]\n",
    "        fraction = re.sub(r\"(.\\d+)$\", \".mzXML\", os.path.dirname(psm))\n",
    "        \n",
    "        if fraction in psms:\n",
    "            psms[fraction].append(psm)\n",
    "        else:\n",
    "            psms[fraction] = [psm]\n",
    "        \n",
    "        if fraction in pep2psm:\n",
    "            if pep in pep2psm[fraction]:\n",
    "                pep2psm[fraction][pep].append(psm)\n",
    "                pep2psm[fraction][pep] = list(set(pep2psm[fraction][pep]))\n",
    "            else:\n",
    "                pep2psm[fraction][pep] = [psm]\n",
    "        else:\n",
    "            pep2psm[fraction] = {}\n",
    "            pep2psm[fraction][pep] = [psm]\n",
    "\n",
    "        if fraction in prot2psm:\n",
    "            if prot in prot2psm[fraction]:\n",
    "                prot2psm[fraction][prot].append(psm)\n",
    "                prot2psm[fraction][prot] = list(set(prot2psm[fraction][prot]))\n",
    "            else:\n",
    "                prot2psm[fraction][prot] = [psm]\n",
    "        else:\n",
    "            prot2psm[fraction] = {}\n",
    "            prot2psm[fraction][prot] = [psm]        \n",
    "\n",
    "#         if pep in pep2psm:\n",
    "#             pep2psm[pep].append(psm)\n",
    "#             pep2psm[pep] = list(set(pep2psm[pep]))\n",
    "#         else:\n",
    "#             pep2psm[pep] = [psm]\n",
    "#         if prot in prot2psm:\n",
    "#             prot2psm[prot].append(psm)\n",
    "#             prot2psm[prot] = list(set(prot2psm[prot]))\n",
    "#         else:\n",
    "#             prot2psm[prot] = [psm]\n",
    "\n",
    "    return psms, pep2psm, prot2psm\n",
    "\n",
    "\n",
    "def getReporterMz(name):\n",
    "    if name == \"sig126\":\n",
    "        return 126.127726\n",
    "    elif name == \"sig127\" or name == \"sig127N\":\n",
    "        return 127.124761\n",
    "    elif name == \"sig127C\":\n",
    "        return 127.131081\n",
    "    elif name == \"sig128N\":\n",
    "        return 128.128116\n",
    "    elif name == \"sig128\" or name == \"sig128C\":\n",
    "        return 128.134436\n",
    "    elif name == \"sig129\" or name == \"sig129N\":\n",
    "        return 129.131471\n",
    "    elif name == \"sig129C\":\n",
    "        return 129.137790\n",
    "    elif name == \"sig130N\":\n",
    "        return 130.134825\n",
    "    elif name == \"sig130\" or name == \"sig130C\":\n",
    "        return 130.141145\n",
    "    elif name == \"sig131\" or name == \"sig131N\":\n",
    "        return 131.138180\n",
    "    elif name == \"sig131C\":\n",
    "        return 131.144500\n",
    "    elif name == \"sig132N\":\n",
    "        return 132.141535\n",
    "    elif name == \"sig132C\":\n",
    "        return 132.147855\n",
    "    elif name == \"sig133N\":\n",
    "        return 133.144890\n",
    "    elif name == \"sig133C\":\n",
    "        return 133.151210\n",
    "    elif name == \"sig134N\":\n",
    "        return 134.148245\n",
    "\n",
    "\n",
    "def getReporterIntensity(scanNum, params, **kwargs):\n",
    "    tol = 10\n",
    "    spec = reader[scanNum]\n",
    "    reporterNames = params[\"tmt_reporters_used\"].split(\";\")\n",
    "\n",
    "    # Extract reporter ion m/z and intensity values\n",
    "    mzArray = []\n",
    "    intensityArray = []\n",
    "    for reporter in reporterNames:\n",
    "        if reporter in kwargs:\n",
    "            mz = getReporterMz(reporter) * (1 + kwargs[reporter] / 1e6)\n",
    "        else:\n",
    "            mz = getReporterMz(reporter)\n",
    "        lL = mz - mz * tol / 1e6\n",
    "        uL = mz + mz * tol / 1e6\n",
    "        ind = np.where((spec[\"m/z array\"] >= lL) & (spec[\"m/z array\"] <= uL))[0]\n",
    "        if len(ind) == 0:\n",
    "            reporterMz = 0\n",
    "        elif len(ind) == 1:\n",
    "            ind = ind[0]\n",
    "            reporterMz = spec[\"m/z array\"][ind]\n",
    "        elif len(ind) > 1:\n",
    "            ind2 = np.argmax(spec[\"intensity array\"][ind])\n",
    "            ind = ind[ind2]\n",
    "            reporterMz = spec[\"m/z array\"][ind]\n",
    "        if lL <= reporterMz < uL:\n",
    "            reporterIntensity = spec[\"intensity array\"][ind]\n",
    "        else:\n",
    "            reporterIntensity = 0\n",
    "        mzArray.append(reporterMz)\n",
    "        intensityArray.append(reporterIntensity)\n",
    "\n",
    "    # Extract K-TMT and R-y1 ion intensities\n",
    "    tmtNo = re.search(\"\\d+\", os.path.basename(params[\"impurity_matrix\"])).group(0)\n",
    "    if tmtNo == \"16\":\n",
    "        Kmz = 451.3199494907  # m/z of K-TMT-y1 ion for TMT16 (= 147.1128041645 + 304.2071453 (TMTpro-balancer))\n",
    "    else:\n",
    "        Kmz = 376.2757362992;  # m/z of K-TMT-y1 ion for TMT10, 11 (= 147.1128... + 229.1629321 (TMT-balancer))\n",
    "\n",
    "    Kintensity = 0\n",
    "    KlL = Kmz - Kmz * tol / 1e6\n",
    "    KuL = Kmz + Kmz * tol / 1e6\n",
    "    ind = np.where((spec[\"m/z array\"] >= KlL) & (spec[\"m/z array\"] <= KuL))[0]\n",
    "    if len(ind) > 1:\n",
    "        ind2 = np.argmax(spec[\"intensity array\"][ind])\n",
    "        ind = ind[ind2]\n",
    "        ind = ind[0]\n",
    "        Kintensity = spec[\"intensity array\"][ind]\n",
    "\n",
    "    Rintensity = 0\n",
    "    Rmz = 175.1189521741\n",
    "    RlL = Rmz - Rmz * tol / 1e6\n",
    "    RuL = Rmz + Rmz * tol / 1e6\n",
    "    ind = np.where((spec[\"m/z array\"] >= RlL) & (spec[\"m/z array\"] <= RuL))[0]\n",
    "    if len(ind) > 1:\n",
    "        ind2 = np.argmax(spec[\"intensity array\"][ind])\n",
    "        ind = ind[ind2]\n",
    "        ind = ind[0]\n",
    "        Rintensity = spec[\"intensity array\"][ind]\n",
    "\n",
    "    outArray = mzArray + intensityArray + [Kintensity, Rintensity]\n",
    "    return outArray\n",
    "\n",
    "\n",
    "def matchMs2ToMs3(psms, reader):\n",
    "    output = {} # Dictionary whose key = MS2 scan number (string), value = corresponding MS3 scan number (string)\n",
    "    for psm in psms:\n",
    "        ms2Num = int(os.path.basename(psm).split(\".\")[-3])  # It should be changed to \"int\" here to deal with numbers such as \"09999\"\n",
    "        ms2PrecMz = reader[str(ms2Num)][\"precursorMz\"][0][\"precursorMz\"]\n",
    "\n",
    "        # Find corresponding MS3 scan number\n",
    "        i = ms2Num + 1\n",
    "        spec = reader[str(i)]\n",
    "        while spec[\"msLevel\"] > 1:\n",
    "            if spec[\"msLevel\"] == 3 and (ms2PrecMz - 0.1) < spec[\"precursorMz\"][-1][\"precursorMz\"] < (ms2PrecMz + 3):\n",
    "                ms3Num = spec[\"num\"]\n",
    "                output[str(ms2Num)] = ms3Num\n",
    "                break\n",
    "            i += 1\n",
    "            spec = reader[str(i)]\n",
    "\n",
    "    return output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys, re, numpy as np, pandas as pd\n",
    "from pyteomics import mzxml, mzml\n",
    "\n",
    "\n",
    "class progressBar:\n",
    "    def __init__(self, total):\n",
    "        self.total = total\n",
    "        self.barLength = 20\n",
    "        self.count = 0\n",
    "        self.progress = 0\n",
    "        self.block = 0\n",
    "        self.status = \"\"\n",
    "\n",
    "    def increment(self):\n",
    "        self.count += 1\n",
    "        self.progress = self.count / self.total\n",
    "        self.block = int(round(self.barLength * self.progress))\n",
    "        if self.progress == 1:\n",
    "            self.status = \"Done...\\r\\n\"\n",
    "        else:\n",
    "            self.status = \"\"\n",
    "        #         self.status = str(self.count) + \"/\" + str(self.total)\n",
    "        text = \"\\r  Progress: [{0}] {1}% {2}\".format(\"#\" * self.block + \"-\" * (self.barLength - self.block),\n",
    "                                                     int(self.progress * 100), self.status)\n",
    "        sys.stdout.write(text)\n",
    "        sys.stdout.flush()\n",
    "\n",
    "\n",
    "def getParams(paramFile):\n",
    "    parameters = dict()\n",
    "    with open(paramFile, 'r') as file:\n",
    "        for line in file:\n",
    "            if re.search(r'^#', line) or re.search(r'^\\s', line):\n",
    "                continue\n",
    "            line = re.sub(r'#.*', '', line)  # Remove comments (start from '#')\n",
    "            line = re.sub(r'\\s*', '', line)  # Remove all whitespaces\n",
    "\n",
    "            # Exception for \"feature_files\" parameter\n",
    "            if \"feature_files\" in parameters and line.endswith(\"feature\"):\n",
    "                parameters[\"feature_files\"].append(line)\n",
    "            else:\n",
    "                key = line.split('=')[0]\n",
    "                val = line.split('=')[1]\n",
    "                if key == \"feature_files\":\n",
    "                    parameters[key] = [val]\n",
    "                else:\n",
    "                    parameters[key] = val\n",
    "    return parameters\n",
    "\n",
    "\n",
    "def parseIdtxt(idTxt, params):\n",
    "    df = pd.read_table(idTxt, sep=\";\", skiprows=0, header=1)\n",
    "    df = df.groupby([\"Peptide\", \"Protein\", \"Outfile\"]).size().reset_index()\n",
    "    psms = {}\n",
    "    pep2psm = {}\n",
    "    prot2psm = {}\n",
    "    for i in range(df.shape[0]):\n",
    "        pep = df[\"Peptide\"][i]\n",
    "        prot = df[\"Protein\"][i]\n",
    "        psm = df[\"Outfile\"][i]\n",
    "        fraction = re.sub(r\"(.\\d+)$\", \".mzXML\", os.path.dirname(psm))\n",
    "        \n",
    "        \n",
    "        \n",
    "        fraction = \"NCI-11plex-1-F1-f10268.mzXML\"\n",
    "        \n",
    "        \n",
    "        if fraction in psms:\n",
    "            psms[fraction].append(psm)\n",
    "            psms[fraction] = list(set(psms[fraction]))\n",
    "        else:\n",
    "            psms[fraction] = [psm]\n",
    "        \n",
    "        if pep in pep2psm:\n",
    "            pep2psm[pep].append(psm)\n",
    "            pep2psm[pep] = list(set(pep2psm[pep]))\n",
    "        else:\n",
    "            pep2psm[pep] = [psm]\n",
    "        \n",
    "        if prot in prot2psm:\n",
    "            prot2psm[prot].append(psm)\n",
    "            prot2psm[prot] = list(set(prot2psm[prot]))\n",
    "        else:\n",
    "            prot2psm[prot] = [psm]\n",
    "        \n",
    "#         if fraction in pep2psm:\n",
    "#             if pep in pep2psm[fraction]:\n",
    "#                 pep2psm[fraction][pep].append(psm)\n",
    "#                 pep2psm[fraction][pep] = list(set(pep2psm[fraction][pep]))\n",
    "#             else:\n",
    "#                 pep2psm[fraction][pep] = [psm]\n",
    "#         else:\n",
    "#             pep2psm[fraction] = {}\n",
    "#             pep2psm[fraction][pep] = [psm]\n",
    "\n",
    "#         if fraction in prot2psm:\n",
    "#             if prot in prot2psm[fraction]:\n",
    "#                 prot2psm[fraction][prot].append(psm)\n",
    "#                 prot2psm[fraction][prot] = list(set(prot2psm[fraction][prot]))\n",
    "#             else:\n",
    "#                 prot2psm[fraction][prot] = [psm]\n",
    "#         else:\n",
    "#             prot2psm[fraction] = {}\n",
    "#             prot2psm[fraction][prot] = [psm]        \n",
    "\n",
    "    return psms, pep2psm, prot2psm\n",
    "\n",
    "\n",
    "# def matchMs2ToMs3(psms, reader):\n",
    "#     output = {} # Dictionary whose key = MS2 scan number (string), value = corresponding MS3 scan number (string)\n",
    "#     for psm in psms:\n",
    "#         ms2Num = int(os.path.basename(psm).split(\".\")[-3])  # It should be changed to \"int\" here to deal with numbers such as \"09999\"\n",
    "#         ms2PrecMz = reader[str(ms2Num)][\"precursorMz\"][0][\"precursorMz\"]\n",
    "\n",
    "#         # Find corresponding MS3 scan number\n",
    "#         i = ms2Num + 1\n",
    "#         spec = reader[str(i)]\n",
    "#         while spec[\"msLevel\"] > 1:\n",
    "#             if spec[\"msLevel\"] == 3 and (ms2PrecMz - 0.1) < spec[\"precursorMz\"][-1][\"precursorMz\"] < (ms2PrecMz + 3):\n",
    "#                 ms3Num = spec[\"num\"]\n",
    "#                 output[str(ms2Num)] = ms3Num\n",
    "#                 break\n",
    "#             i += 1\n",
    "#             spec = reader[str(i)]\n",
    "\n",
    "#     return output\n",
    "\n",
    "def matchMs2ToMs3(psms):\n",
    "    output = {} # Dictionary whose key = MS2 scan number (string), value = corresponding MS3 scan number (string)\n",
    "    for fraction in psms.keys():\n",
    "        output[fraction] = {}\n",
    "        reader = mzxml.MzXML(fraction)    \n",
    "        for psm in psms[fraction]:\n",
    "            ms2Num = int(os.path.basename(psm).split(\".\")[-3])  # It should be changed to \"int\" here to deal with numbers such as \"09999\"\n",
    "            ms2PrecMz = reader[str(ms2Num)][\"precursorMz\"][0][\"precursorMz\"]\n",
    "\n",
    "            # Find corresponding MS3 scan number\n",
    "            i = ms2Num + 1\n",
    "            spec = reader[str(i)]\n",
    "            while spec[\"msLevel\"] > 1:\n",
    "                if spec[\"msLevel\"] == 3 and (ms2PrecMz - 0.1) < spec[\"precursorMz\"][-1][\"precursorMz\"] < (ms2PrecMz + 3):\n",
    "                    ms3Num = spec[\"num\"]\n",
    "                    output[fraction][ms2Num] = ms3Num\n",
    "                    break\n",
    "                i += 1\n",
    "                spec = reader[str(i)]\n",
    "\n",
    "    return output\n",
    "\n",
    "\n",
    "\n",
    "def getReporterMz(name):\n",
    "    if name == \"sig126\":\n",
    "        return 126.127726\n",
    "    elif name == \"sig127\" or name == \"sig127N\":\n",
    "        return 127.124761\n",
    "    elif name == \"sig127C\":\n",
    "        return 127.131081\n",
    "    elif name == \"sig128N\":\n",
    "        return 128.128116\n",
    "    elif name == \"sig128\" or name == \"sig128C\":\n",
    "        return 128.134436\n",
    "    elif name == \"sig129\" or name == \"sig129N\":\n",
    "        return 129.131471\n",
    "    elif name == \"sig129C\":\n",
    "        return 129.137790\n",
    "    elif name == \"sig130N\":\n",
    "        return 130.134825\n",
    "    elif name == \"sig130\" or name == \"sig130C\":\n",
    "        return 130.141145\n",
    "    elif name == \"sig131\" or name == \"sig131N\":\n",
    "        return 131.138180\n",
    "    elif name == \"sig131C\":\n",
    "        return 131.144500\n",
    "    elif name == \"sig132N\":\n",
    "        return 132.141535\n",
    "    elif name == \"sig132C\":\n",
    "        return 132.147855\n",
    "    elif name == \"sig133N\":\n",
    "        return 133.144890\n",
    "    elif name == \"sig133C\":\n",
    "        return 133.151210\n",
    "    elif name == \"sig134N\":\n",
    "        return 134.148245\n",
    "\n",
    "\n",
    "def getReporterIntensity(scanNum, params, **kwargs):\n",
    "    tol = 10\n",
    "    spec = reader[scanNum]\n",
    "    reporterNames = params[\"tmt_reporters_used\"].split(\";\")\n",
    "\n",
    "    # Extract reporter ion m/z and intensity values\n",
    "    mzArray = []\n",
    "    intensityArray = []\n",
    "    for reporter in reporterNames:\n",
    "        if reporter in kwargs:\n",
    "            mz = getReporterMz(reporter) * (1 + kwargs[reporter] / 1e6)\n",
    "        else:\n",
    "            mz = getReporterMz(reporter)\n",
    "        lL = mz - mz * tol / 1e6\n",
    "        uL = mz + mz * tol / 1e6\n",
    "        ind = np.where((spec[\"m/z array\"] >= lL) & (spec[\"m/z array\"] <= uL))[0]\n",
    "        if len(ind) == 0:\n",
    "            reporterMz = 0\n",
    "        elif len(ind) == 1:\n",
    "            ind = ind[0]\n",
    "            reporterMz = spec[\"m/z array\"][ind]\n",
    "        elif len(ind) > 1:\n",
    "            ind2 = np.argmax(spec[\"intensity array\"][ind])\n",
    "            ind = ind[ind2]\n",
    "            reporterMz = spec[\"m/z array\"][ind]\n",
    "        if lL <= reporterMz < uL:\n",
    "            reporterIntensity = spec[\"intensity array\"][ind]\n",
    "        else:\n",
    "            reporterIntensity = 0\n",
    "        mzArray.append(reporterMz)\n",
    "        intensityArray.append(reporterIntensity)\n",
    "\n",
    "    # Extract K-TMT and R-y1 ion intensities\n",
    "    tmtNo = re.search(\"\\d+\", os.path.basename(params[\"impurity_matrix\"])).group(0)\n",
    "    if tmtNo == \"16\":\n",
    "        Kmz = 451.3199494907  # m/z of K-TMT-y1 ion for TMT16 (= 147.1128041645 + 304.2071453 (TMTpro-balancer))\n",
    "    else:\n",
    "        Kmz = 376.2757362992;  # m/z of K-TMT-y1 ion for TMT10, 11 (= 147.1128... + 229.1629321 (TMT-balancer))\n",
    "\n",
    "    Kintensity = 0\n",
    "    KlL = Kmz - Kmz * tol / 1e6\n",
    "    KuL = Kmz + Kmz * tol / 1e6\n",
    "    ind = np.where((spec[\"m/z array\"] >= KlL) & (spec[\"m/z array\"] <= KuL))[0]\n",
    "    if len(ind) > 1:\n",
    "        ind2 = np.argmax(spec[\"intensity array\"][ind])\n",
    "        ind = ind[ind2]\n",
    "        ind = ind[0]\n",
    "        Kintensity = spec[\"intensity array\"][ind]\n",
    "\n",
    "    Rintensity = 0\n",
    "    Rmz = 175.1189521741\n",
    "    RlL = Rmz - Rmz * tol / 1e6\n",
    "    RuL = Rmz + Rmz * tol / 1e6\n",
    "    ind = np.where((spec[\"m/z array\"] >= RlL) & (spec[\"m/z array\"] <= RuL))[0]\n",
    "    if len(ind) > 1:\n",
    "        ind2 = np.argmax(spec[\"intensity array\"][ind])\n",
    "        ind = ind[ind2]\n",
    "        ind = ind[0]\n",
    "        Rintensity = spec[\"intensity array\"][ind]\n",
    "\n",
    "    outArray = mzArray + intensityArray + [Kintensity, Rintensity]\n",
    "    return outArray\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mzxmlFile = \"NCI-11plex-1-F1-f10268.mzXML\"\n",
    "paramFile = \"jump_qc.params\"\n",
    "# idTxt = \"./sum_NCI-11Plex-1-F1-f10268/ID.txt\"\n",
    "idTxt = \"ID.txt\"\n",
    "\n",
    "\n",
    "params = getParams(paramFile)\n",
    "psms, pep2psm, prot2psm = parseIdtxt(idTxt, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms2ToMs3 = matchMs2ToMs3(psms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 697,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  sig126\t985 (98.50%) matched\n",
      "  sig127N\t975 (97.50%) matched\n",
      "  sig127C\t967 (96.70%) matched\n",
      "  sig128N\t985 (98.50%) matched\n",
      "  sig128C\t982 (98.20%) matched\n",
      "  sig129N\t981 (98.10%) matched\n",
      "  sig129C\t988 (98.80%) matched\n",
      "  sig130N\t973 (97.30%) matched\n",
      "  sig130C\t985 (98.50%) matched\n",
      "  sig131N\t960 (96.00%) matched\n",
      "  sig131C\t993 (99.30%) matched\n"
     ]
    }
   ],
   "source": [
    "#####################################\n",
    "# In case of MS3 quantification,    #\n",
    "# MS2 - MS3 scans need to be paired #\n",
    "#####################################\n",
    "ms2ToMs3 = matchMs2ToMs3(psms)\n",
    "\n",
    "####################################################\n",
    "# Extract TMT reporter ion intensities - 1st round #\n",
    "####################################################\n",
    "qdict = {}  # Dictionary; key = MS2 scan number, value = TMT reporter intensities (array)\n",
    "for fraction in sorted(ms2ToMs3.keys()):\n",
    "    print(\"  TMT reporter ion extraction from %s\" % os.path.basename(fraction))\n",
    "    for ms2, ms3 in ms2ToMs3[fraction].items():\n",
    "        reporterIntensity = getReporterIntensity(ms3, params)\n",
    "        qdict[ms2] = reporterIntensity\n",
    "\n",
    "# Create a dataFrame after the first extraction of reporter m/z and intensity values\n",
    "reporters = params[\"tmt_reporters_used\"].split(\";\")\n",
    "columnNames = [re.sub(\"sig\", \"mz\", i) for i in reporters] + reporters  + [\"K_y1\", \"R_y1\"]\n",
    "df = pd.DataFrame.from_dict(qdict, orient = 'index', columns = columnNames) # Pandas dataframe of TMT reporter intensities\n",
    "\n",
    "# Calculate m/z shift of each reporter ion (to be used for refining reporter intensities)\n",
    "mzShifts = {}\n",
    "nTot = df.shape[0]\n",
    "for reporter in reporters:\n",
    "    reporterMz = getReporterMz(reporter)\n",
    "    measuredMz = df[re.sub(\"sig\", \"mz\", reporter)][df[re.sub(\"sig\", \"mz\", reporter)] > 0]   # Exclude un-quantified ones\n",
    "    n = len(measuredMz)\n",
    "    print(\"  %s\\t%d (%.2f%%) matched\" % (reporter, n, n / nTot * 100))\n",
    "    mzShift = ((measuredMz - reporterMz) / reporterMz * 1e6).mean()\n",
    "    mzShifts[reporter] = mzShift\n",
    "\n",
    "##########################################################\n",
    "# Refine the reporter intensities considering m/z shifts #\n",
    "##########################################################\n",
    "qdict = {}\n",
    "for fraction in sorted(ms2ToMs3.keys()):\n",
    "    for ms2, ms3 in ms2ToMs3[fraction].items():\n",
    "        reporterIntensity = getReporterIntensity(ms3, params, **mzShifts)\n",
    "        qdict[ms2] = reporterIntensity\n",
    "df = pd.DataFrame.from_dict(qdict, orient = 'index', columns = columnNames)\n",
    "\n",
    "# TMT impurity correction\n",
    "dfImpurity = pd.read_table(\"TMT11.ini\", sep=\"\\t\", skiprows=1, header=None, index_col = 0)\n",
    "dfImpurity = pd.DataFrame(np.linalg.pinv(dfImpurity.values), dfImpurity.columns, dfImpurity.index)\n",
    "dfImpurity.columns = [str(\"sig\") + i for i in dfImpurity.columns]\n",
    "df2 = df[reporters].dot(dfImpurity.T)\n",
    "df2.columns = reporters\n",
    "df2[df2 <= 1] = 0 # Lower bound to 0 (1 is a margin to prevent underflow when taking log2)\n",
    "df[reporters] = pd.concat([df[reporters]/2, df2]).max(level = 0) # Upper bound to the half of original values (i.e. before correction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "# Generate raw_..._(non)zero.txt files\n",
    "# For this implementation, \n",
    "# 1. No zero-intensity filtering (keep all zero-intensity PSMs)\n",
    "# 2. No PPI filtering (PPI information is not available in Comet results)\n",
    "filterMethods = params[\"min_intensity_method\"].split(\",\")\n",
    "filterCutoffs = params[\"min_intensity_value\"].split(\",\")\n",
    "for i in range(len(filterMethods)):\n",
    "    if filterMethods[i] == \"1\":\n",
    "        df = df[df[reporters].min(axis = 1) >= float(filterCutoffs[i])]\n",
    "    elif filterMethods[i] == \"2\":\n",
    "        df = df[df[reporters].max(axis = 1) >= float(filterCutoffs[i])]\n",
    "    elif filterMethods[i] == \"3\":\n",
    "        df = df[df[reporters].mean(axis = 1) >= float(filterCutoffs[i])]\n",
    "    elif filterMethods[i] == \"4\":\n",
    "        df = df[df[reporters].median(axis = 1) >= float(filterCutoffs[i])]\n",
    "\n",
    "# For this implementation, min_intensity_..._psm_1_2_filtered is not considered\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 698,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Reporter\tMean[%]\tSD[%]\tSEM[%]\t#PSMs\n",
      "  sig126\t124.40\t39.64\t3.01\t173\n",
      "  sig127N\t104.01\t25.92\t1.97\t173\n",
      "  sig127C\t123.99\t21.20\t1.61\t173\n",
      "  sig128N\t145.34\t23.39\t1.78\t173\n",
      "  sig128C\t65.85\t19.85\t1.51\t173\n",
      "  sig129N\t62.71\t25.22\t1.92\t173\n",
      "  sig129C\t112.67\t19.60\t1.49\t173\n",
      "  sig130N\t72.80\t24.50\t1.86\t173\n",
      "  sig130C\t88.70\t26.61\t2.02\t173\n",
      "  sig131N\t44.70\t48.86\t3.71\t173\n",
      "  sig131C\t119.64\t12.18\t0.93\t173\n"
     ]
    }
   ],
   "source": [
    "# Loading-bias calculation\n",
    "# \"subDf\" is used to calculate the loading-bias informatio and normalization factors\n",
    "noiseLevel = 1000\n",
    "snRatio = params[\"SNratio_for_correction\"]\n",
    "\n",
    "snRatio = 0\n",
    "\n",
    "subDf = df[reporters][(df[reporters] > noiseLevel * snRatio).prod(axis = 1).astype(bool)] # Zero-intensity PSMs are excluded\n",
    "\n",
    "# Loading-bias calculation using trimmed-mean\n",
    "# - Divide each PSM by its row-wise mean and then take log2\n",
    "# - In this way, non-variant PSMs will have values close to 0\n",
    "subDf = np.log2(subDf.divide(subDf.mean(axis = 1), axis = 0))\n",
    "pctTrimmed = float(params[\"percentage_trimmed\"])\n",
    "n = 0\n",
    "for reporter in reporters:\n",
    "    if n == 0:\n",
    "        ind = ((subDf[reporter] > subDf[reporter].quantile(pctTrimmed / 200)) & (subDf[reporter] < subDf[reporter].quantile(1 - pctTrimmed / 200)))\n",
    "    else:\n",
    "        ind = ind & ((subDf[reporter] > subDf[reporter].quantile(pctTrimmed / 200)) & (subDf[reporter] < subDf[reporter].quantile(1 - pctTrimmed / 200)))\n",
    "    n += 1\n",
    "\n",
    "meanIntensity = 2 ** subDf.loc[ind].mean(axis = 0) * 100\n",
    "# To express SD as percentage in raw-scale, take the mean deviation of 2^sd and 2^(-sd)\n",
    "# For example, sdVal = 0.2 for a specific reporter which represent +/-0.2 in log2-scale\n",
    "# 2^0.2 = 1.1487, 2^(-0.2) = 0.87\n",
    "# positive variation = 2^0.2 - 1 = 0.1487 (assuming the center is ideally 1 (=2^0))\n",
    "# negative variation = 1 - 2^(-0.2) = 0.13\n",
    "# mean variation (in raw-scale) = (0.1487 + 0.13) / 2 = 0.1394\n",
    "# Let the standard deviation of the reporter be 13.94% (= 0.1394 * 100)\n",
    "sdVal = subDf.loc[ind].std(axis = 0)\n",
    "sdIntensity = ((2 ** sdVal - 1) + (1 - 2 ** (-sdVal))) / 2 * 100\n",
    "semIntensity = sdIntensity / np.sqrt(sum(ind))\n",
    "\n",
    "print(\"  Reporter\\tMean[%]\\tSD[%]\\tSEM[%]\\t#PSMs\")\n",
    "for i in range(len(reporters)):\n",
    "    print (\"  %s\\t%.2f\\t%.2f\\t%.2f\\t%d\" % (reporters[i], meanIntensity[i], sdIntensity[i], semIntensity[i], sum(ind)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 688,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            mz126      mz127N      mz127C      mz128N      mz128C      mz129N  \\\n",
      "1743   126.127235  127.124329  127.130867  128.127960  128.134018  129.131119   \n",
      "1785   126.127487  127.124008  127.130920  128.127579  128.133911  129.131271   \n",
      "1788   126.127571  127.124374  127.130913  128.127899    0.000000  129.131073   \n",
      "2060     0.000000    0.000000    0.000000    0.000000    0.000000    0.000000   \n",
      "2112   126.127518  127.124268  127.130753    0.000000  128.133911    0.000000   \n",
      "...           ...         ...         ...         ...         ...         ...   \n",
      "14921  126.127609  127.124680  127.130775  128.128006  128.134094  129.131531   \n",
      "14937  126.127510  127.124596  127.130760  128.127914  128.133987  129.131485   \n",
      "14941  126.127586  127.124825  127.130882  128.127975  128.133972  129.131393   \n",
      "14952  126.127480  127.124565  127.130623  128.127991  128.134109  129.131332   \n",
      "14973  126.127487  127.124657  127.130806  128.127914  128.134048  129.131302   \n",
      "\n",
      "           mz129C      mz130N      mz130C      mz131N  ...       sig128N  \\\n",
      "1743   129.137329  130.134293  130.140656  131.137726  ...   2044.381734   \n",
      "1785   129.137405  130.134293  130.140686  131.137970  ...    561.443993   \n",
      "1788   129.137360  130.134537  130.140717    0.000000  ...    543.972847   \n",
      "2060   129.137405    0.000000  130.140686    0.000000  ...      0.000000   \n",
      "2112   129.137268  130.134567    0.000000  131.137741  ...      0.000000   \n",
      "...           ...         ...         ...         ...  ...           ...   \n",
      "14921  129.137558  130.134735  130.140839  131.138351  ...  53206.414593   \n",
      "14937  129.137482  130.134659  130.140793  131.138092  ...  53390.740661   \n",
      "14941  129.137512  130.134735  130.140854  131.138290  ...  67331.770002   \n",
      "14952  129.137482  130.134293  130.140884  131.138138  ...   2416.674745   \n",
      "14973  129.137482  130.134689  130.140839  131.138184  ...  21906.558425   \n",
      "\n",
      "             sig128C        sig129N        sig129C        sig130N  \\\n",
      "1743     3379.231051    2928.529738    2043.698608    3092.957022   \n",
      "1785     1061.983259    1587.065454    1145.233827    4167.096802   \n",
      "1788       16.724380    2784.771606    1120.073456    2148.020972   \n",
      "2060        0.000000       0.000000     460.290168       0.000000   \n",
      "2112     1832.210870       0.000000     765.925995    1604.269402   \n",
      "...              ...            ...            ...            ...   \n",
      "14921  138808.358810  103290.967283  102023.664489  137268.377124   \n",
      "14937  107256.360374   70685.334735   87525.412386   81546.431070   \n",
      "14941   89093.745675  168159.281285   64398.802810  122558.553495   \n",
      "14952   40053.641758   69342.188015   30620.560695   12426.817921   \n",
      "14973   44382.944895   54536.832651   24802.260306   38743.748696   \n",
      "\n",
      "            sig130C        sig131N       sig131C  K_y1  R_y1  \n",
      "1743    2409.287769    4339.184921   1934.980214     0     0  \n",
      "1785    2597.538358    3372.293094   1359.644951     0     0  \n",
      "1788    1432.183007       0.000000   3500.339876     0     0  \n",
      "2060     733.320637       0.000000    317.830928     0     0  \n",
      "2112       0.000000    2886.668857   1444.806548     0     0  \n",
      "...             ...            ...           ...   ...   ...  \n",
      "14921  53633.894901  164450.375845  93280.097409     0     0  \n",
      "14937  49047.615365  131020.412994  42841.556746     0     0  \n",
      "14941  79347.911259  117312.699168  65494.197043     0     0  \n",
      "14952  86479.385619  210270.321222  65123.025830     0     0  \n",
      "14973  27579.374157   45738.106063  19187.698380     0     0  \n",
      "\n",
      "[1000 rows x 24 columns]\n"
     ]
    }
   ],
   "source": [
    "# Normalization (i.e. loading-bias correction)\n",
    "\n",
    "doNormalization = params[\"loading_bias_correction\"]\n",
    "normalizationMethod = params[\"loading_bias_correction_method\"]\n",
    "if doNormalization == \"1\":\n",
    "    if normalizationMethod == \"1\": # Trimmed-mean\n",
    "        normFactor = subDf.loc[ind].mean(axis = 0) - np.mean(subDf.loc[ind].mean())\n",
    "    elif normalizationMethod == \"2\": # Trimmed-median\n",
    "        normFactor = subDf.loc[ind].median(axis = 0) - np.mean(subDf.loc[ind].median())\n",
    "    \n",
    "    # Normalize the dataframe, df\n",
    "    psmMean = df[reporters].mean(axis = 1)\n",
    "    df[reporters] = np.log2(df[reporters].divide(psmMean, axis = 0).replace(0, np.nan))\n",
    "    df[reporters] = df[reporters] - normFactor\n",
    "    df[reporters] = 2 ** df[reporters]\n",
    "    df[reporters] = df[reporters].multiply(psmMean, axis = 0).replace(np.nan, 0)\n",
    "    \n",
    "print(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 738,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Summarization of PSMs to peptides/proteins\n",
    "pepDict = {}\n",
    "for pep, psms in pep2psm.items():\n",
    "    psms = [i.split(\".\")[-3] for i in psms]\n",
    "    subDf = df[reporters][df.index.isin(psms)]    \n",
    "    if subDf.empty:\n",
    "        continue\n",
    "    else:\n",
    "        # 1. Take log2 again for the summarization\n",
    "        # 2. Representative mean values of PSMs\n",
    "        # 3. Mean-centering across reporters and then take mean over PSMs\n",
    "        # 4. Add back the representative mean \n",
    "        subDf = np.log2(df[reporters][df.index.isin(psms)].replace(0, np.nan))\n",
    "        psmMean = subDf.mean(axis = 1).sort_values(ascending = False)\n",
    "        if len(psmMean) > 3:\n",
    "            repMean = np.mean(psmMean[0:3])\n",
    "        else:\n",
    "            repMean = np.mean(psmMean)\n",
    "        pepDict[pep] = (2 ** (subDf.sub(psmMean, axis = 0).mean(axis = 0) + repMean)).replace(np.nan, 0).tolist()\n",
    "\n",
    "protDict = {}\n",
    "for prot, psms in prot2psm.items():\n",
    "    psms = [i.split(\".\")[-3] for i in psms]\n",
    "    subDf = df[reporters][df.index.isin(psms)]    \n",
    "    if subDf.empty:\n",
    "        continue\n",
    "    else:\n",
    "        # 1. Take log2 again for the summarization\n",
    "        # 2. Representative mean values of PSMs\n",
    "        # 3. Mean-centering across reporters and then take mean over PSMs\n",
    "        # 4. Add back the representative mean \n",
    "        subDf = np.log2(df[reporters][df.index.isin(psms)].replace(0, np.nan))\n",
    "        psmMean = subDf.mean(axis = 1).sort_values(ascending = False)\n",
    "        if len(psmMean) > 3:\n",
    "            repMean = np.mean(psmMean[0:3])\n",
    "        else:\n",
    "            repMean = np.mean(psmMean)\n",
    "        protDict[prot] = (2 ** (subDf.sub(psmMean, axis = 0).mean(axis = 0) + repMean)).replace(np.nan, 0).tolist()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 735,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 698.4754145768168,\n",
       " 0.0,\n",
       " 689.6569468194513,\n",
       " 0.0,\n",
       " 543.8897821056173]"
      ]
     },
     "execution_count": 735,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psms = [\"1743\", \"1788\", \"2060\", \"2112\"]\n",
    "psms = [\"2060\"]\n",
    "subDf = np.log2(df[reporters][df.index.isin(psms)].replace(0, np.nan))\n",
    "psmMean = subDf.mean(axis = 1).sort_values(ascending = False)\n",
    "if len(psmMean) > 3:\n",
    "    repMean = np.mean(psmMean[0:3])\n",
    "else:\n",
    "    repMean = np.mean(psmMean)\n",
    "\n",
    "(2 ** (subDf.sub(psmMean, axis = 0).mean(axis = 0) + repMean)).replace(np.nan, 0).tolist()\n",
    "# subDf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 742,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[8672.600123957794,\n",
       " 42598.32396432284,\n",
       " 33058.185444339666,\n",
       " 10582.131903809353,\n",
       " 36174.17373936337,\n",
       " 11333.91503107109,\n",
       " 44119.07268432173,\n",
       " 90619.14183245438,\n",
       " 40836.19071467221,\n",
       " 20224.52576923634,\n",
       " 34820.11387530317]"
      ]
     },
     "execution_count": 742,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "protDict[\"sp|P14136|GFAP_HUMAN\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 712,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sig126</th>\n",
       "      <th>sig127N</th>\n",
       "      <th>sig127C</th>\n",
       "      <th>sig128N</th>\n",
       "      <th>sig128C</th>\n",
       "      <th>sig129N</th>\n",
       "      <th>sig129C</th>\n",
       "      <th>sig130N</th>\n",
       "      <th>sig130C</th>\n",
       "      <th>sig131N</th>\n",
       "      <th>sig131C</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1743</th>\n",
       "      <td>2491.761347</td>\n",
       "      <td>2853.252386</td>\n",
       "      <td>3846.551398</td>\n",
       "      <td>5162.429670</td>\n",
       "      <td>1751.503829</td>\n",
       "      <td>1376.757657</td>\n",
       "      <td>3101.246413</td>\n",
       "      <td>1959.612601</td>\n",
       "      <td>2265.832931</td>\n",
       "      <td>1036.492093</td>\n",
       "      <td>3311.244677</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1788</th>\n",
       "      <td>65951.405895</td>\n",
       "      <td>901.663526</td>\n",
       "      <td>2276.582764</td>\n",
       "      <td>1373.628768</td>\n",
       "      <td>8.668486</td>\n",
       "      <td>1309.174219</td>\n",
       "      <td>1699.675175</td>\n",
       "      <td>1360.927079</td>\n",
       "      <td>1346.907356</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>5989.974316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2060</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>698.475415</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>689.656947</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>543.889782</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2112</th>\n",
       "      <td>9497.996337</td>\n",
       "      <td>911.319322</td>\n",
       "      <td>830.762767</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>949.661123</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1162.267878</td>\n",
       "      <td>1016.421021</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>689.532597</td>\n",
       "      <td>2472.432513</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            sig126      sig127N      sig127C      sig128N      sig128C  \\\n",
       "1743   2491.761347  2853.252386  3846.551398  5162.429670  1751.503829   \n",
       "1788  65951.405895   901.663526  2276.582764  1373.628768     8.668486   \n",
       "2060      0.000000     0.000000     0.000000     0.000000     0.000000   \n",
       "2112   9497.996337   911.319322   830.762767     0.000000   949.661123   \n",
       "\n",
       "          sig129N      sig129C      sig130N      sig130C      sig131N  \\\n",
       "1743  1376.757657  3101.246413  1959.612601  2265.832931  1036.492093   \n",
       "1788  1309.174219  1699.675175  1360.927079  1346.907356     0.000000   \n",
       "2060     0.000000   698.475415     0.000000   689.656947     0.000000   \n",
       "2112     0.000000  1162.267878  1016.421021     0.000000   689.532597   \n",
       "\n",
       "          sig131C  \n",
       "1743  3311.244677  \n",
       "1788  5989.974316  \n",
       "2060   543.889782  \n",
       "2112  2472.432513  "
      ]
     },
     "execution_count": 712,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[reporters][df.index.isin(psms)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "psms = ['/research/rgs01/project_space/penggrp/BrainTumor_JP/penggrp/surendhar/MS3quan_JiHoon/NCI-11plex-1-F1-f10268/NCI-11plex-1-F1-f10268.1/NCI-11plex-1-F1-f10268.26125.26125.3.out']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['/research/rgs01/project_space/penggrp/BrainTumor_JP/penggrp/surendhar/MS3quan_JiHoon/NCI-11plex-1-F1-f10268/NCI-11plex-1-F1-f10268.1/NCI-11plex-1-F1-f10268.26125.26125.3.out']"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "can't assign to function call (<ipython-input-74-5d83a064a0dc>, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-74-5d83a064a0dc>\"\u001b[1;36m, line \u001b[1;32m1\u001b[0m\n\u001b[1;33m    a = [i for os.path.basename(i) in psms]\u001b[0m\n\u001b[1;37m              ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m can't assign to function call\n"
     ]
    }
   ],
   "source": [
    "a = [i for os.path.basename(i) in psms]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = [os.path.basename(i).split(\".\")[0] + \".mzXML\" + \"_\" + os.path.basename(i).split(\".\")[-3] for i in psms]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_table(\"ID.txt\", sep=\";\", skiprows=0, header=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/research/rgs01/project_space/penggrp/BrainTumor_JP/penggrp/surendhar/MS3quan_JiHoon/sum_Batch1_sample1/\n"
     ]
    }
   ],
   "source": [
    "with open('ID.txt') as f:\n",
    "    jumpfPath = f.readline()\n",
    "\n",
    "jumpfPath = re.sub(\"misc/idsum.db\", \"\", jumpfPath.split(\"=\")[1]).strip()\n",
    "print(jumpfPath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                Protein Group#  \\\n",
      "Protein Accession #                              \n",
      "co|CON_A2A4G1|Con               SJPG000001.001   \n",
      "co|CON_A2A5Y0|Con               SJPG000001.002   \n",
      "co|CON_ENSP00000377550|Con      SJPG000001.004   \n",
      "co|CON_P02533|Con               SJPG000001.007   \n",
      "co|CON_P05784|Con               SJPG000001.009   \n",
      "...                                        ...   \n",
      "sp|P04350|TBB4A_HUMAN           SJPG119534.001   \n",
      "sp|P68371|TBB4B_HUMAN           SJPG119535.001   \n",
      "sp|P16989|YBOX3_HUMAN           SJPG119841.001   \n",
      "sp|Q6ZNC4|ZN704_HUMAN           SJPG120261.001   \n",
      "tr|A0A087WXU7|A0A087WXU7_HUMAN  SJPG120301.006   \n",
      "\n",
      "                                                              Protein Description  \\\n",
      "Protein Accession #                                                                 \n",
      "co|CON_A2A4G1|Con               TREMBL:A2A4G1|Contaminant Tax_Id=10090 Gene_Sy...   \n",
      "co|CON_A2A5Y0|Con               TREMBL:A2A5Y0|Contaminant Tax_Id=10090 Gene_Sy...   \n",
      "co|CON_ENSP00000377550|Con           Tax_Id=9606 Gene_Symbol=KRT13 46 kDa protein   \n",
      "co|CON_P02533|Con               SWISS-PROT:P02533|Contaminant Tax_Id=9606 Gene...   \n",
      "co|CON_P05784|Con               SWISS-PROT:P05784|Contaminant Tax_Id=10090 Gen...   \n",
      "...                                                                           ...   \n",
      "sp|P04350|TBB4A_HUMAN           Tubulin beta-4A chain OS=Homo sapiens GN=TUBB4...   \n",
      "sp|P68371|TBB4B_HUMAN           Tubulin beta-4B chain OS=Homo sapiens GN=TUBB4...   \n",
      "sp|P16989|YBOX3_HUMAN           Y-box-binding protein 3 OS=Homo sapiens GN=YBX...   \n",
      "sp|Q6ZNC4|ZN704_HUMAN           Zinc finger protein 704 OS=Homo sapiens GN=ZNF...   \n",
      "tr|A0A087WXU7|A0A087WXU7_HUMAN  Zinc finger protein 761 OS=Homo sapiens GN=ZNF...   \n",
      "\n",
      "                                    GN  PSM#  Total Peptide#  Unique Peptide#  \\\n",
      "Protein Accession #                                                             \n",
      "co|CON_A2A4G1|Con                  NaN     4               4                1   \n",
      "co|CON_A2A5Y0|Con                  NaN     1               1                0   \n",
      "co|CON_ENSP00000377550|Con         NaN     3               2                0   \n",
      "co|CON_P02533|Con                  NaN     6               5                0   \n",
      "co|CON_P05784|Con                  NaN     1               1                0   \n",
      "...                                ...   ...             ...              ...   \n",
      "sp|P04350|TBB4A_HUMAN           TUBB4A    41              11                1   \n",
      "sp|P68371|TBB4B_HUMAN           TUBB4B    44              13                1   \n",
      "sp|P16989|YBOX3_HUMAN             YBX3     3               1                0   \n",
      "sp|Q6ZNC4|ZN704_HUMAN           ZNF704     1               1                0   \n",
      "tr|A0A087WXU7|A0A087WXU7_HUMAN  ZNF761     1               1                1   \n",
      "\n",
      "                                %Coverage  Peptide of the Highest Score  \\\n",
      "Protein Accession #                                                       \n",
      "co|CON_A2A4G1|Con                     9.47                  TRLEQEIATYR   \n",
      "co|CON_A2A5Y0|Con                     2.17                    ILDELTLCK   \n",
      "co|CON_ENSP00000377550|Con            5.49                AGLENTVAETECR   \n",
      "co|CON_P02533|Con                    11.46              GSCGIGGGIGGGSSR   \n",
      "co|CON_P05784|Con                     1.66                      SLETENR   \n",
      "...                                    ...                          ...   \n",
      "sp|P04350|TBB4A_HUMAN                28.15   SGPFGQIFRPDNFVFGQSGAGNNWAK   \n",
      "sp|P68371|TBB4B_HUMAN                32.58   SGPFGQIFRPDNFVFGQSGAGNNWAK   \n",
      "sp|P16989|YBOX3_HUMAN                 4.03              DGVPEGAQLQGPVHR   \n",
      "sp|Q6ZNC4|ZN704_HUMAN                 4.13            VTAAMVLTSLSTSPLVR   \n",
      "tr|A0A087WXU7|A0A087WXU7_HUMAN       27.78             DVM*LENYRNLVSLDR   \n",
      "\n",
      "                                                  Run#  Scan#  ...  ppm  \\\n",
      "Protein Accession #                                            ...        \n",
      "co|CON_A2A4G1|Con               NCI-11plex-1-F1-f10268  20608  ...  0.0   \n",
      "co|CON_A2A5Y0|Con               NCI-11plex-1-F1-f10268  37888  ...  0.0   \n",
      "co|CON_ENSP00000377550|Con      NCI-11plex-1-F1-f10268  17921  ...  0.0   \n",
      "co|CON_P02533|Con               NCI-11plex-1-F1-f10268  11223  ...  0.0   \n",
      "co|CON_P05784|Con               NCI-11plex-1-F1-f10268   5959  ...  0.0   \n",
      "...                                                ...    ...  ...  ...   \n",
      "sp|P04350|TBB4A_HUMAN           NCI-11plex-1-F1-f10268  43373  ...  0.0   \n",
      "sp|P68371|TBB4B_HUMAN           NCI-11plex-1-F1-f10268  43373  ...  0.0   \n",
      "sp|P16989|YBOX3_HUMAN           NCI-11plex-1-F1-f10268  18277  ...  0.0   \n",
      "sp|Q6ZNC4|ZN704_HUMAN           NCI-11plex-1-F1-f10268  46105  ...  0.0   \n",
      "tr|A0A087WXU7|A0A087WXU7_HUMAN  NCI-11plex-1-F1-f10268  45888  ...  0.0   \n",
      "\n",
      "                                Xcorr    dCn  Q-value(%)  GPCRs  TFs  \\\n",
      "Protein Accession #                                                    \n",
      "co|CON_A2A4G1|Con               3.645  0.650      1.0000      -    -   \n",
      "co|CON_A2A5Y0|Con               2.697  0.413      2.3906      -    -   \n",
      "co|CON_ENSP00000377550|Con      4.207  0.721      1.0000      -    -   \n",
      "co|CON_P02533|Con               4.001  0.642      1.0000      -    -   \n",
      "co|CON_P05784|Con               2.018  0.313      1.0000      -    -   \n",
      "...                               ...    ...         ...    ...  ...   \n",
      "sp|P04350|TBB4A_HUMAN           6.658  0.702      1.0000      -    -   \n",
      "sp|P68371|TBB4B_HUMAN           6.658  0.702      1.0000      -    -   \n",
      "sp|P16989|YBOX3_HUMAN           4.240  0.631      2.3906      -    -   \n",
      "sp|Q6ZNC4|ZN704_HUMAN           2.435  0.148      2.3906      -    -   \n",
      "tr|A0A087WXU7|A0A087WXU7_HUMAN  2.251  0.275      2.3906      -    -   \n",
      "\n",
      "                               epigenetic_factors kinases oncogenes  \\\n",
      "Protein Accession #                                                   \n",
      "co|CON_A2A4G1|Con                               -       -         -   \n",
      "co|CON_A2A5Y0|Con                               -       -         -   \n",
      "co|CON_ENSP00000377550|Con                      -       -         -   \n",
      "co|CON_P02533|Con                               -       -         -   \n",
      "co|CON_P05784|Con                               -       -         -   \n",
      "...                                           ...     ...       ...   \n",
      "sp|P04350|TBB4A_HUMAN                           -       -         -   \n",
      "sp|P68371|TBB4B_HUMAN                           -       -         -   \n",
      "sp|P16989|YBOX3_HUMAN                           -       -         -   \n",
      "sp|Q6ZNC4|ZN704_HUMAN                           -       -         -   \n",
      "tr|A0A087WXU7|A0A087WXU7_HUMAN                  -       -         -   \n",
      "\n",
      "                               spliceosomal_proteins  \n",
      "Protein Accession #                                   \n",
      "co|CON_A2A4G1|Con                                  -  \n",
      "co|CON_A2A5Y0|Con                                  -  \n",
      "co|CON_ENSP00000377550|Con                         -  \n",
      "co|CON_P02533|Con                                  -  \n",
      "co|CON_P05784|Con                                  -  \n",
      "...                                              ...  \n",
      "sp|P04350|TBB4A_HUMAN                              -  \n",
      "sp|P68371|TBB4B_HUMAN                              -  \n",
      "sp|P16989|YBOX3_HUMAN                              O  \n",
      "sp|Q6ZNC4|ZN704_HUMAN                              -  \n",
      "tr|A0A087WXU7|A0A087WXU7_HUMAN                     -  \n",
      "\n",
      "[2468 rows x 22 columns]\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "df = pd.read_table(\"id_uni_prot.txt\", sep = \"\\t\", skiprows = 0, header = 1, index_col = 1)\n",
    "print(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "infoColumns = [\"ProteinDescription\", \"GN\", r\"PSM#\", r\"TotalPeptide#\", r\"UniquePeptide#\", r\"%Coverage\",\n",
    "               \"PeptideoftheHighestScore\", r\"Run#\", r\"Scan#\", \"m/z\", \"z\", \"ppm\", \"Xcorr\", \"dCn\", r\"Q-value(%)\"]\n",
    "df.columns = df.columns.str.replace(' ', '')\n",
    "infoDf = df[infoColumns]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = {'col1': [1, 2, 4], 'col2': [3, 4, 7]}\n",
    "df1 = pd.DataFrame(data=d, index=[\"a\", \"c\", \"b\"])\n",
    "d = {'col1': [\"d\", \"e\", \"f\"], 'col2': [\"q\", \"w\", \"e\"]}\n",
    "df2 = pd.DataFrame(data=d, index=[\"c\", \"a\", \"e\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   col1  col2\n",
      "a     1     3\n",
      "c     2     4\n",
      "b     4     7\n",
      "  col1 col2\n",
      "c    d    q\n",
      "a    e    w\n",
      "e    f    e\n"
     ]
    }
   ],
   "source": [
    "print(df1)\n",
    "print(df2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True,  True, False])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df1.index.isin(df2.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "df3 = pd.concat([df1, df2], axis=1, join = \"inner\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>col1</th>\n",
       "      <th>col2</th>\n",
       "      <th>col1</th>\n",
       "      <th>col2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>a</th>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>e</td>\n",
       "      <td>w</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>c</th>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>d</td>\n",
       "      <td>q</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  col1 col2 col1 col2\n",
       "a    1    3    e    w\n",
       "c    2    4    d    q"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>col1_x</th>\n",
       "      <th>col2_x</th>\n",
       "      <th>col1_y</th>\n",
       "      <th>col2_y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>a</th>\n",
       "      <td>1.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>e</td>\n",
       "      <td>w</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b</th>\n",
       "      <td>4.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>c</th>\n",
       "      <td>2.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>d</td>\n",
       "      <td>q</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>e</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>f</td>\n",
       "      <td>e</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   col1_x  col2_x col1_y col2_y\n",
       "a     1.0     3.0      e      w\n",
       "b     4.0     7.0    NaN    NaN\n",
       "c     2.0     4.0      d      q\n",
       "e     NaN     NaN      f      e"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}