Post-Translational Modification (PTM)¶

This tutorial shows hwo to process and analyse PTM data using msmu. We will use a mock oxidation-site dataset for demonstration.

Load Required Libraries¶

In [1]:

Copied!

import msmu as mm
import msmu as mm

Data Preparation¶

In this tutorial, we will use peptide-level processed data from SAGE TMT workflow for both global protein and PTM data with code below.

More details about Sage TMT data processing can be found in DDA - TMT tutorial section.

In [2]:

Copied!





base_dir = "https://raw.githubusercontent.com/bertis-informatics/msmu/refs/heads/dev/data/sage_tmt"
sage_idents = f"{base_dir}/sage/results.sage.tsv"
sage_quants = f"{base_dir}/sage/tmt.tsv"
sdrf = f"{base_dir}/meta.sdrf.tsv"

ptm_mdata = mm.read_sage(identification_file=sage_idents, quantification_file=sage_quants, label="tmt")
# Sample metadata from the SDRF: TMT channels in obs already line up with comment[label]
ptm_mdata = mm.pp.attach_sdrf(ptm_mdata, sdrf)
ptm_mdata = mm.pp.apply_sdrf_to_obs(ptm_mdata)

# Without an SDRF, add_meta() joins any table onto obs instead -- here a CSV whose "tag" column
# holds the same channel names as obs:
# ptm_mdata = mm.pp.add_meta(ptm_mdata, f"{base_dir}/meta.csv", format="csv", metadata_on="tag")

ptm_mdata = mm.pp.add_filter(ptm_mdata, modality="psm", column="q_value", keep="lt", value=0.01)
ptm_mdata = mm.pp.apply_filter(ptm_mdata, modality="psm")

ptm_mdata = mm.pp.to_peptide(ptm_mdata)
ptm_mdata = mm.pp.log2_transform(ptm_mdata, modality="peptide")
ptm_mdata = mm.pp.normalise(ptm_mdata, modality="peptide", method="median")

ptm_mdata = mm.pp.add_filter(ptm_mdata, modality="peptide", column="q_value", keep="lt", value=0.01)
ptm_mdata = mm.pp.apply_filter(ptm_mdata, modality="peptide")

ptm_mdata
base_dir = "https://raw.githubusercontent.com/bertis-informatics/msmu/refs/heads/dev/data/sage_tmt"
sage_idents = f"{base_dir}/sage/results.sage.tsv"
sage_quants = f"{base_dir}/sage/tmt.tsv"
sdrf = f"{base_dir}/meta.sdrf.tsv"

ptm_mdata = mm.read_sage(identification_file=sage_idents, quantification_file=sage_quants, label="tmt")
# Sample metadata from the SDRF: TMT channels in obs already line up with comment[label]
ptm_mdata = mm.pp.attach_sdrf(ptm_mdata, sdrf)
ptm_mdata = mm.pp.apply_sdrf_to_obs(ptm_mdata)

# Without an SDRF, add_meta() joins any table onto obs instead -- here a CSV whose "tag" column
# holds the same channel names as obs:
# ptm_mdata = mm.pp.add_meta(ptm_mdata, f"{base_dir}/meta.csv", format="csv", metadata_on="tag")

ptm_mdata = mm.pp.add_filter(ptm_mdata, modality="psm", column="q_value", keep="lt", value=0.01)
ptm_mdata = mm.pp.apply_filter(ptm_mdata, modality="psm")

ptm_mdata = mm.pp.to_peptide(ptm_mdata)
ptm_mdata = mm.pp.log2_transform(ptm_mdata, modality="peptide")
ptm_mdata = mm.pp.normalise(ptm_mdata, modality="peptide", method="median")

ptm_mdata = mm.pp.add_filter(ptm_mdata, modality="peptide", column="q_value", keep="lt", value=0.01)
ptm_mdata = mm.pp.apply_filter(ptm_mdata, modality="peptide")

ptm_mdata

INFO - Validating SDRF metadata for https://raw.githubusercontent.com/bertis-informatics/msmu/refs/heads/dev/data/sage_tmt/meta.sdrf.tsv.
INFO - SDRF validation succeeded for https://raw.githubusercontent.com/bertis-informatics/msmu/refs/heads/dev/data/sage_tmt/meta.sdrf.tsv.
WARNING - apply_sdrf_to_obs: psm columns not projectable under 'comment[label]'; kept only in uns['sdrf']: ['assay name', 'comment[data file]', 'comment[fraction identifier]']
INFO - Applying var filters for psm: ['q_value_lt_0.01']
WARNING - Purity column not found in psm modality for TMT data. Skipping purity filtering.
INFO - Peptide-level identifications: 2242 (2189 at 1% FDR)
INFO - Building new peptide quantification data.
INFO - Applying var filters for peptide: ['q_value_lt_0.01']

Out[2]:

MuData object with n_obs × n_vars = 6 × 4499
  obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[biological replicate]', 'technology type', 'comment[proteomexchange accession number]', 'comment[proteomics data acquisition method]', 'comment[technical replicate]', 'comment[instrument]', 'comment[cleavage agent details]', 'factor value[condition]', 'factor value[time point]'
  uns:	'_cmd', 'sdrf'
  2 modalities
    psm:	6 × 2310
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[biological replicate]', 'technology type', 'comment[proteomexchange accession number]', 'comment[proteomics data acquisition method]', 'comment[technical replicate]', 'comment[instrument]', 'comment[cleavage agent details]', 'factor value[condition]', 'factor value[time point]'
      var:	'proteins', 'peptide', 'stripped_peptide', 'filename', 'scan_num', 'charge', 'peptide_length', 'expmass', 'calcmass', 'rt', 'missed_cleavages', 'semi_enzymatic', 'contaminant', 'PEP', 'score', 'q_value', 'decoy'
      uns:	'level', 'search_engine', 'quantification', 'label', 'acquisition', 'identification_file', 'quantification_file', 'decoy', 'filter', 'decoy_filter'
      varm:	'search_result', 'filter'
      layers:	None
    peptide:	6 × 2189
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[biological replicate]', 'technology type', 'comment[proteomexchange accession number]', 'comment[proteomics data acquisition method]', 'comment[technical replicate]', 'comment[instrument]', 'comment[cleavage agent details]', 'factor value[condition]', 'factor value[time point]'
      var:	'peptide', 'proteins', 'stripped_peptide', 'count_psm', 'PEP', 'q_value'
      uns:	'level', 'decoy', 'filter', 'decoy_filter'
      varm:	'filter'
      layers:	None

Global Protein Data Processing¶

For real data analysis, global protein data should be processed with other batch experiment for global protein quantification. Here, we will use the same SAGE TMT data for demonstration purpose.

In [3]:

Copied!





global_mdata = ptm_mdata.copy()
global_mdata = mm.pp.infer_protein(global_mdata)
global_mdata = mm.pp.to_protein(global_mdata)
global_mdata = mm.pp.add_filter(global_mdata, modality="protein", column="q_value", keep="lt", value=0.01)
global_mdata = mm.pp.apply_filter(global_mdata, modality="protein")
global_mdata = ptm_mdata.copy()
global_mdata = mm.pp.infer_protein(global_mdata)
global_mdata = mm.pp.to_protein(global_mdata)
global_mdata = mm.pp.add_filter(global_mdata, modality="protein", column="q_value", keep="lt", value=0.01)
global_mdata = mm.pp.apply_filter(global_mdata, modality="protein")

INFO - Starting protein inference
INFO - Initial proteins: 1772
INFO - Removed indistinguishable: 205
INFO - Removed subsettable: 68
INFO - Removed subsumable: 0
INFO - Total protein groups: 1499
INFO - Applying var filters for peptide: ['q_value_lt_0.01', 'peptide_type_eq_unique']
INFO - Protein-level identifications: 1479 (1446 at 1% FDR)
INFO - Applying var filters for protein: ['q_value_lt_0.01']

PTM Data Processing¶

Protein Inference with global protein data¶

PTM data uses the same peptide-protein mapping as global protein data, therefore we can directly use the peptide to protein mapping from global protein data processing step by using mm.pp.infer_protein function with progagated_from parameter.

In [4]:

Copied!

ptm_mdata = mm.pp.infer_protein(ptm_mdata, propagated_from=global_mdata)
ptm_mdata = mm.pp.infer_protein(ptm_mdata, propagated_from=global_mdata)

INFO - Starting protein inference

PTM Summarization¶

For PTM summarization, we need FASTA file that contains the protein sequences used for ptm site localization.

Then, PTM site is summarized by mm.pp.to_ptm with modi_name for specifying modality name (e.g. "oxidation" -> "oxidation_site") for PTM data and modification for specifying the modification string used for site localization in modified peptides.

In [5]:

Copied!

import urllib.request

fasta_url = "https://raw.githubusercontent.com/bertis-informatics/msmu/refs/heads/main/data/fasta/human_fasta.fasta"
fasta_file = "./human_fasta.fasta"

with urllib.request.urlopen(fasta_url) as response:
    with open(fasta_file, "wb") as f:
        fasta_ = response.read()
        f.write(fasta_)

ptm_mdata = mm.utils.attach_fasta(ptm_mdata, fasta_file=fasta_file)

ptm_mdata = mm.pp.to_ptm(ptm_mdata, modi_name="oxidation", modification="[+15.9949]")

ptm_mdata
import urllib.request

fasta_url = "https://raw.githubusercontent.com/bertis-informatics/msmu/refs/heads/main/data/fasta/human_fasta.fasta"
fasta_file = "./human_fasta.fasta"

with urllib.request.urlopen(fasta_url) as response:
    with open(fasta_file, "wb") as f:
        fasta_ = response.read()
        f.write(fasta_)

ptm_mdata = mm.utils.attach_fasta(ptm_mdata, fasta_file=fasta_file)

ptm_mdata = mm.pp.to_ptm(ptm_mdata, modi_name="oxidation", modification="[+15.9949]")

ptm_mdata

INFO - oxidation site level identifications: 355
INFO - Building new oxidation_site AnnData.

Out[5]:

MuData object with n_obs × n_vars = 6 × 4854
  obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[biological replicate]', 'technology type', 'comment[proteomexchange accession number]', 'comment[proteomics data acquisition method]', 'comment[technical replicate]', 'comment[instrument]', 'comment[cleavage agent details]', 'factor value[condition]', 'factor value[time point]'
  uns:	'_cmd', 'sdrf', 'protein_map', 'protein_info'
  3 modalities
    psm:	6 × 2310
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[biological replicate]', 'technology type', 'comment[proteomexchange accession number]', 'comment[proteomics data acquisition method]', 'comment[technical replicate]', 'comment[instrument]', 'comment[cleavage agent details]', 'factor value[condition]', 'factor value[time point]'
      var:	'proteins', 'peptide', 'stripped_peptide', 'filename', 'scan_num', 'charge', 'peptide_length', 'expmass', 'calcmass', 'rt', 'missed_cleavages', 'semi_enzymatic', 'contaminant', 'PEP', 'score', 'q_value', 'decoy'
      uns:	'level', 'search_engine', 'quantification', 'label', 'acquisition', 'identification_file', 'quantification_file', 'decoy', 'filter', 'decoy_filter'
      varm:	'search_result', 'filter'
      layers:	None
    peptide:	6 × 2189
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[biological replicate]', 'technology type', 'comment[proteomexchange accession number]', 'comment[proteomics data acquisition method]', 'comment[technical replicate]', 'comment[instrument]', 'comment[cleavage agent details]', 'factor value[condition]', 'factor value[time point]'
      var:	'peptide', 'proteins', 'stripped_peptide', 'count_psm', 'PEP', 'q_value', 'protein_group', 'peptide_type'
      uns:	'level', 'decoy', 'filter', 'decoy_filter'
      varm:	'filter'
      layers:	None
    oxidation_site:	6 × 355
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[biological replicate]', 'technology type', 'comment[proteomexchange accession number]', 'comment[proteomics data acquisition method]', 'comment[technical replicate]', 'comment[instrument]', 'comment[cleavage agent details]', 'factor value[condition]', 'factor value[time point]'
      var:	'count_psm', 'peptide', 'count_peptide', 'count_stripped_peptide', 'modified_protein', 'protein_group'
      uns:	'level'
      layers:	None

In [6]:

Copied!

ptm_mdata["oxidation_site"].to_df().T
ptm_mdata["oxidation_site"].to_df().T

Out[6]:

	TMT126	TMT127	TMT128	TMT129	TMT130	TMT131
A0A3B3IRV3\|M150,Q9ULC4\|M150	16.287181	15.531237	15.366565	15.043427	15.300467	15.785685
O00410\|M61	14.649085	14.665745	14.835433	14.423701	14.756286	15.716622
O00422\|M120	14.854015	15.049636	14.587649	14.446713	14.671983	14.602903
O14497\|M1273	13.220737	13.392085	13.252274	13.060778	12.828074	14.319759
O14950\|M40,P19105\|M39	15.827595	15.779716	15.652113	15.583274	15.564402	15.772828
...	...	...	...	...	...	...
Q9Y490\|M289	17.157233	17.170126	16.595825	16.162770	15.988541	16.272613
Q9Y490\|M72	14.405495	14.174609	13.832705	13.533238	13.454351	13.631721
Q9Y5B9\|M468	14.919705	15.107410	15.241040	14.694629	14.658880	15.913531
Q9Y639\|M379	14.107615	13.768589	13.299487	13.211008	13.219137	13.623166
Q9Y6G3\|M108	10.566706	12.021548	10.485821	12.230562	10.888014	11.915218

355 rows × 6 columns

PTM quantification adjustment by global protein levels¶

To adjust PTM quantification by global protein levels, we can use mm.pp.adjust_ptm_by_protein function. This function normalizes PTM quantification by corresponding protein levels to remove the confounding effect from protein abundance changes.

In [7]:

Copied!





ptm_mdata = mm.pp.adjust_ptm_by_protein(
    ptm_mdata,
    modality="oxidation_site",
    global_mdata=global_mdata,
    method="ridge",
    rescale=True,
)
ptm_mdata = mm.pp.adjust_ptm_by_protein(
    ptm_mdata,
    modality="oxidation_site",
    global_mdata=global_mdata,
    method="ridge",
    rescale=True,
)

In [8]:

Copied!

ptm_mdata["oxidation_site"].to_df().T
ptm_mdata["oxidation_site"].to_df().T

Out[8]:

	TMT126	TMT127	TMT128	TMT129	TMT130	TMT131
A0A3B3IRV3\|M150,Q9ULC4\|M150	14.935351	14.186535	14.023415	13.703324	13.957940	14.438584
O00410\|M61	14.025286	14.028702	14.196342	13.793879	14.128571	15.072368
O00422\|M120	14.359037	14.554203	14.093290	13.952682	14.177428	14.108508
O14497\|M1273	14.084276	14.253384	14.115401	13.926408	13.696746	15.168933
O14950\|M40,P19105\|M39	14.338397	14.290678	14.163049	14.094256	14.075232	14.283536
...	...	...	...	...	...	...
Q9Y490\|M289	14.799231	14.811959	14.245012	13.817503	13.645505	13.925939
Q9Y490\|M72	14.768749	14.537743	14.201190	13.905757	13.828493	14.003217
Q9Y5B9\|M468	14.037130	14.234196	14.362815	13.825602	13.781173	15.004232
Q9Y639\|M379	14.773293	14.436458	13.970387	13.882480	13.890556	14.291974
Q9Y6G3\|M108	13.433358	14.872647	13.348597	15.077745	13.746582	14.766219

338 rows × 6 columns