Post-translational modification (PTM)¶

This tutorial shows how to process and analyze 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

Load a recorded result¶

Complete DDA-TMT with the current msmu version, then run this notebook from the same working directory. Its processed result, dda_tmt_PXD013361.h5mu, contains peptide and protein modalities on a log2 scale.

In [2]:

Copied!

ptm_mdata = mm.read_h5mu("dda_tmt_PXD013361.h5mu")

ptm_mdata
ptm_mdata = mm.read_h5mu("dda_tmt_PXD013361.h5mu")

ptm_mdata

Out[2]:

MuData object with n_obs × n_vars = 6 × 5945
  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:	'_log', 'sdrf'
  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:	'acquisition', 'decoy', 'decoy_filter', 'filter', 'label', 'level', 'normalisation', 'quantification', 'search_engine'
      obsm:	'filter'
      varm:	'filter', 'search_result'
      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:	'decoy', 'decoy_filter', 'filter', 'level'
      varm:	'filter'
      layers:	None
    protein:	6 × 1446
      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', 'count_stripped_peptide', 'PEP', 'q_value'
      uns:	'decoy', 'decoy_filter', 'filter', 'level'
      varm:	'filter'
      layers:	None

Reuse global protein data¶

For this demonstration, use the processed protein modality from the same DDA-TMT result as the global protein reference. In a real PTM experiment, process the matching global protein dataset separately and use that result for adjustment.

In [3]:

Copied!

global_mdata = ptm_mdata.copy()
global_mdata = ptm_mdata.copy()

PTM data processing¶

PTM site summarization¶

PTM sites are localized from each peptide's own protein accessions and the attached FASTA, so no protein inference is needed on the PTM side. Site IDs therefore depend only on the peptide and the FASTA -- the same PTM data yields the same sites whether or not a global dataset is involved.

We need the FASTA file that contains the protein sequences used for 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 [4]:

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 - 2 contaminant accessions are not in the attached FASTA (Cont_P13645, Cont_P62937). Search engines add their own contaminant entries, so this is expected unless the same contaminants were part of the search database.

INFO - oxidation site level identifications: 326

INFO - Building new oxidation_site AnnData.

Out[4]:

MuData object with n_obs × n_vars = 6 × 6271
  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:	'_log', 'sdrf'
  4 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:	'acquisition', 'decoy', 'decoy_filter', 'filter', 'label', 'level', 'normalisation', 'quantification', 'search_engine'
      obsm:	'filter'
      varm:	'filter', 'search_result'
      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:	'decoy', 'decoy_filter', 'filter', 'level'
      varm:	'filter'
      layers:	None
    protein:	6 × 1446
      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', 'count_stripped_peptide', 'PEP', 'q_value'
      uns:	'decoy', 'decoy_filter', 'filter', 'level'
      varm:	'filter'
      layers:	None
    oxidation_site:	6 × 326
      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', 'count_site'
      uns:	'level'
      layers:	None

In [5]:

Copied!

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

Out[5]:

	TMT126	TMT127	TMT128	TMT129	TMT130	TMT131
A0A2R8Y4L2\|M186;P09651\|M186;Q32P51\|M186	NaN	NaN	NaN	NaN	9.791745	NaN
A0A3B3IRV3\|M150;Q9ULC4\|M150	16.146860	15.515804	15.319607	14.977872	15.435319	15.841854
A6NNZ2\|M257;P04350\|M257;P07437\|M257;P68371\|M257;Q13509\|M257;Q13885\|M257;Q3ZCM7\|M257;Q9BUF5\|M257;Q9BVA1\|M257;Q9H4B7\|M257	10.556790	11.907977	11.633446	11.897772	12.119013	13.536121
O00410\|M61	14.561682	14.649298	14.851874	14.467341	14.707265	15.723757
O00422\|M120	14.766611	15.033188	14.604090	14.490353	14.622963	14.610038
...	...	...	...	...	...	...
Q9Y490\|M289	17.238530	17.288507	16.707468	16.144944	15.994228	16.350998
Q9Y490\|M72	14.318091	14.158161	13.849147	13.576878	13.405330	13.638856
Q9Y5B9\|M468	14.942548	15.090657	15.185116	14.713860	14.679580	15.903320
Q9Y639\|M379	14.139623	13.679764	13.345722	13.290486	13.325169	13.463748
Q9Y6G3\|M108	10.491659	12.002548	10.412292	12.359474	10.773333	11.961554

326 rows × 6 columns

PTM quantification adjustment by global protein levels¶

mm.pp.adjust_ptm_by_protein() computes differential PTM usage: each site's log intensity minus its parent protein's, so a site's change reads relative to whatever its protein did.

A site's denominator is found through the accessions it was localized on, translated into a global protein group via the global dataset's protein_map. Nothing is looked up by peptide, so a PTM peptide the global run never observed is still adjustable whenever its protein was quantified there -- the common case under enrichment.

The adjusted values replace the quantification that was read, as with every other msmu transform. Sites that could not be adjusted become NaN, and .var["adjustment_status"] records why for each one. Use mm.dt.save_layer() first if you want to keep the unadjusted values alongside.

In [6]:

Copied!





# Keep the unadjusted values so both analyses stay available.
ptm_mdata = mm.dt.save_layer(ptm_mdata, modality="oxidation_site", layer="unadjusted")

ptm_mdata = mm.pp.adjust_ptm_by_protein(
    ptm_mdata,
    modality="oxidation_site",
    global_mdata=global_mdata,
)
# Keep the unadjusted values so both analyses stay available.
ptm_mdata = mm.dt.save_layer(ptm_mdata, modality="oxidation_site", layer="unadjusted")

ptm_mdata = mm.pp.adjust_ptm_by_protein(
    ptm_mdata,
    modality="oxidation_site",
    global_mdata=global_mdata,
)

INFO - PTM denominator resolution: 315 / 326 sites adjustable (96.6%)

INFO -   unadjusted [not_quantified]: 6

INFO -   unadjusted [shared_groups]: 5

INFO -   unadjusted [no_estimate]: 2 (denominator found, but the estimator produced no value)

INFO - Protein-adjusted quantification written to oxidation_site.X; 13 of 326 sites are NaN because they could not be adjusted.

In [7]:

Copied!

# Adjusted values are in .X; the pre-adjustment ones remain in the 'unadjusted' layer.
ptm_mdata["oxidation_site"].to_df().T
# Adjusted values are in .X; the pre-adjustment ones remain in the 'unadjusted' layer.
ptm_mdata["oxidation_site"].to_df().T

Out[7]:

	TMT126	TMT127	TMT128	TMT129	TMT130	TMT131
A0A2R8Y4L2\|M186;P09651\|M186;Q32P51\|M186	NaN	NaN	NaN	NaN	11.970668	NaN
A0A3B3IRV3\|M150;Q9ULC4\|M150	14.325371	14.325371	14.325371	14.325371	14.325371	14.325371
A6NNZ2\|M257;P04350\|M257;P07437\|M257;P68371\|M257;Q13509\|M257;Q13885\|M257;Q3ZCM7\|M257;Q9BUF5\|M257;Q9BVA1\|M257;Q9H4B7\|M257	NaN	NaN	NaN	NaN	NaN	NaN
O00410\|M61	17.594138	16.653357	16.482979	17.078261	17.657197	17.007952
O00422\|M120	14.325371	14.325371	14.325371	14.325371	14.325371	14.325371
...	...	...	...	...	...	...
Q9Y490\|M289	18.278125	18.235015	18.008291	17.634992	17.700149	17.893160
Q9Y490\|M72	15.357687	15.104669	15.149970	15.066927	15.111251	15.181018
Q9Y5B9\|M468	13.671135	14.072903	14.030695	14.133666	13.727967	13.843130
Q9Y639\|M379	14.325371	14.325371	14.325371	14.325371	14.325371	14.325371
Q9Y6G3\|M108	12.078847	12.531111	11.700796	12.558135	11.730972	12.362716

326 rows × 6 columns