DDA - TMT¶

This tutorial demonstrates how to analyze Tandem Mass Tag (TMT) labeled data from data-dependent acquisition (DDA) using msmu package.

Data Preparation¶

Original dataset is from PXD013361 (Magnusson et al., 2019) and the search was performed with Sage v0.14.7.

For demonstration purposes, the example dataset has been trimmed to include only 5,000 PSMs.

In [1]:

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"

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"

Load Required Packages¶

If you haven't installed the msmu package yet, please follow the installation guide.

In [2]:

import msmu as mm
import plotly.io as pio

pio.renderers.default = "png"

import msmu as mm import plotly.io as pio pio.renderers.default = "png"

Read Data¶

You can read data from various proteomics software outputs. Below are examples for Sage, MaxQuant, and FragPipe formats.

For this tutorial, we will use the Sage output as an example.

read_sage() function reads the Sage output files (tmt.tsv, results.sage.tsv) and creates modalities at MuData object.

In [3]:

# Sage format
mdata = mm.read_sage(
    identification_file=sage_idents,
    quantification_file=sage_quants,
    sdrf_file=sdrf,
    label="tmt",
)

# # MaxQuant format
# mdata = mm.read_maxquant(
#     identification_file="path_to_maxquant_output",
#     sdrf_file=sdrf,
#     label="tmt",
#     acquisition="dda",
# )

# # FragPipe format
# mdata = mm.read_fragpipe(
#     identification_file="path_to_fragpipe_output",
#     sdrf_file=sdrf,
#     label="tmt",
#     acquisition="dda",
# )

# Sage format mdata = mm.read_sage( identification_file=sage_idents, quantification_file=sage_quants, sdrf_file=sdrf, label="tmt", ) # # MaxQuant format # mdata = mm.read_maxquant( # identification_file="path_to_maxquant_output", # sdrf_file=sdrf, # label="tmt", # acquisition="dda", # ) # # FragPipe format # mdata = mm.read_fragpipe( # identification_file="path_to_fragpipe_output", # sdrf_file=sdrf, # label="tmt", # acquisition="dda", # )

In [4]:

mdata

mdata

Out[4]:

MuData object with n_obs × n_vars = 6 × 3805
  obs:	'source_name', 'characteristics_organism', 'characteristics_organism_part', 'characteristics_disease', 'characteristics_biological_replicate', 'assay_name', 'technology_type', 'comment_proteomexchange_accession_number', 'comment_proteomics_data_acquisition_method', 'comment_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', 'comment_instrument', 'comment_cleavage_agent_details', 'factor_value_condition', 'factor_value_time_point'
  uns:	'sdrf', '_cmd'
  1 modality
    psm:	6 x 3805
      obs:	'source_name', 'characteristics_organism', 'characteristics_organism_part', 'characteristics_disease', 'characteristics_biological_replicate', 'assay_name', 'technology_type', 'comment_proteomexchange_accession_number', 'comment_proteomics_data_acquisition_method', 'comment_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', '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', 'sdrf'
      varm:	'search_result'

Removing Blank Channels¶

If your TMT data contains blank channels, you can remove them by subsetting MuData object.

In [5]:

mdata = mm.pp.add_filter(
    mdata,
    modality="psm",
    column="factor_value_condition",
    keep="not_contains",
    value="BLANK",
    on="obs",
)
mdata = mm.pp.apply_filter(mdata, modality="psm")

mdata = mm.pp.add_filter( mdata, modality="psm", column="factor_value_condition", keep="not_contains", value="BLANK", on="obs", ) mdata = mm.pp.apply_filter(mdata, modality="psm")

INFO - Applying obs filters for psm: ['factor_value_condition_not_contains_BLANK']

Saving Raw MuData Object¶

After reading the search output, saving the MuData object as an h5mu file is recommended for future use.

In [6]:

mdata.write_h5mu("dda_tmt_PXD013361_raw.h5mu")

mdata.write_h5mu("dda_tmt_PXD013361_raw.h5mu")

Handling PSM level¶

Filtering - PSM¶

You can filter the data based on the column values, such as q-value. You can also filter the data based on string containment, which can be useful for removing contaminants or decoys.

Filtering is split into two steps: first, you mark a filter condition using mm.pp.add_filter(), and then you apply the filter using mm.pp.apply_filter().

Here, we keep protein groups with q-value < 0.01 and remove contaminants (protein IDs containing "contam_").

In [7]:

mdata = mm.pp.add_filter(mdata, modality="psm", column="q_value", keep="lt", value=0.01)
mdata = mm.pp.add_filter(mdata, modality="psm", column="proteins", keep="not_contains", value="contam_")
mdata = mm.pp.apply_filter(mdata, modality="psm")

mdata

mdata = mm.pp.add_filter(mdata, modality="psm", column="q_value", keep="lt", value=0.01) mdata = mm.pp.add_filter(mdata, modality="psm", column="proteins", keep="not_contains", value="contam_") mdata = mm.pp.apply_filter(mdata, modality="psm") mdata

INFO - Applying var filters for psm: ['q_value_lt_0.01', 'proteins_not_contains_contam_']
INFO - Applying obs filters for psm: ['factor_value_condition_not_contains_BLANK']

Out[7]:

MuData object with n_obs × n_vars = 6 × 2264
  obs:	'source_name', 'characteristics_organism', 'characteristics_organism_part', 'characteristics_disease', 'characteristics_biological_replicate', 'assay_name', 'technology_type', 'comment_proteomexchange_accession_number', 'comment_proteomics_data_acquisition_method', 'comment_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', 'comment_instrument', 'comment_cleavage_agent_details', 'factor_value_condition', 'factor_value_time_point'
  uns:	'sdrf', '_cmd'
  1 modality
    psm:	6 x 2264
      obs:	'source_name', 'characteristics_organism', 'characteristics_organism_part', 'characteristics_disease', 'characteristics_biological_replicate', 'assay_name', 'technology_type', 'comment_proteomexchange_accession_number', 'comment_proteomics_data_acquisition_method', 'comment_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', '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', 'sdrf', 'filter', 'decoy_filter'
      obsm:	'filter'
      varm:	'search_result', 'filter'

Normalization¶

Here, we log2 transform and normalize the data at the PSM level.

Median centering normalization is applied using mm.pp.normalize() function.

Optionally, if you need to correct batch effects using Global Internal Standard (GIS) channels, you can scale the data using mm.pp.scale_feature() function. Make sure to have GIS channels in each TMT batch. The example code below assumes that the GIS channels are named with prefix "POOLED_".

In [8]:

mdata = mm.pp.log2_transform(mdata, modality="psm")
mdata = mm.pp.normalize(mdata, modality="psm", method="median", fraction=True)

# GIS-based scaling (if needed)
# mdata = mm.pp.correct_batch_effect(mdata=mdata, modality="psm", method="gis", gis_prefix="POOLED_")

mdata = mm.pp.log2_transform(mdata, modality="psm") mdata = mm.pp.normalize(mdata, modality="psm", method="median", fraction=True) # GIS-based scaling (if needed) # mdata = mm.pp.correct_batch_effect(mdata=mdata, modality="psm", method="gis", gis_prefix="POOLED_")

Computing precursor isolation purity¶

For TMT data, it is recommended to compute precursor isolation purity to exclude low-purity features from quantification aggregation. compute_precursor_isolation_purity() function calculates precursor isolation purity using mzML files.

For demonstration purposes, we will skip this step in the tutorial.

In [9]:

# Computing precursor isolation purity
# mdata = mm.pp.compute_precursor_isolation_purity(mdata, mzml_paths=["path/to/sample1.mzML", "path/to/sample2.mzML", ...])

# Computing precursor isolation purity # mdata = mm.pp.compute_precursor_isolation_purity(mdata, mzml_paths=["path/to/sample1.mzML", "path/to/sample2.mzML", ...])

Handling peptide level¶

Summarization - peptide¶

You can summarize psm-level data to peptide-level data using the mm.pp.to_peptide() function.

If isolation purity was calculated and purity_threshold parameter is set, features below this cutoff are excluded from quantification aggregation.

In [10]:

mdata = mm.pp.to_peptide(mdata, purity_threshold=0.7)

# If precursor purity is computed, filter low-purity features
# mdata = mm.pp.to_peptide(mdata, purity_threshold=0.7)

mdata = mm.pp.to_peptide(mdata, purity_threshold=0.7) # If precursor purity is computed, filter low-purity features # mdata = mm.pp.to_peptide(mdata, purity_threshold=0.7)

WARNING - Purity column not found in psm modality for TMT data. Skipping purity filtering.
INFO - Peptide-level identifications: 2204 (2151 at 1% FDR)
INFO - Building new peptide quantification data.

Filtering - peptide¶

In [11]:

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

mdata

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

INFO - Applying var filters for peptide: ['q_value_lt_0.01']

Out[11]:

MuData object with n_obs × n_vars = 6 × 4415
  obs:	'source_name', 'characteristics_organism', 'characteristics_organism_part', 'characteristics_disease', 'characteristics_biological_replicate', 'assay_name', 'technology_type', 'comment_proteomexchange_accession_number', 'comment_proteomics_data_acquisition_method', 'comment_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', 'comment_instrument', 'comment_cleavage_agent_details', 'factor_value_condition', 'factor_value_time_point'
  uns:	'sdrf', '_cmd'
  2 modalities
    psm:	6 x 2264
      obs:	'source_name', 'characteristics_organism', 'characteristics_organism_part', 'characteristics_disease', 'characteristics_biological_replicate', 'assay_name', 'technology_type', 'comment_proteomexchange_accession_number', 'comment_proteomics_data_acquisition_method', 'comment_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', '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', 'sdrf', 'filter', 'decoy_filter'
      obsm:	'filter'
      varm:	'search_result', 'filter'
    peptide:	6 x 2151
      obs:	'source_name', 'characteristics_organism', 'characteristics_organism_part', 'characteristics_disease', 'characteristics_biological_replicate', 'assay_name', 'technology_type', 'comment_proteomexchange_accession_number', 'comment_proteomics_data_acquisition_method', 'comment_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', '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'

Protein inference¶

You can infer protein-level data from peptide-level data using the mm.pp.infer_protein() function.

In [12]:

mdata = mm.pp.infer_protein(mdata)

mdata = mm.pp.infer_protein(mdata)

INFO - Starting protein inference
INFO - Initial proteins: 1722
INFO - Removed indistinguishable: 175
INFO - Removed subsettable: 62
INFO - Removed subsumable: 0
INFO - Total protein groups: 1485

Handling protein level¶

Summarization - protein¶

You can summarize peptide-level data to protein-level data using the mm.pp.to_protein() function.

As default, top 3 peptides within protein group can be used for protein group quantification aggregation. If top_n is None, all peptides will be used.

In [13]:

mdata = mm.pp.to_protein(mdata, top_n=3, rank_method="median_intensity")

mdata = mm.pp.to_protein(mdata, top_n=3, rank_method="median_intensity")

INFO - Applying var filters for peptide: ['q_value_lt_0.01', 'peptide_type_eq_unique']
INFO - Protein-level identifications: 1465 (1432 at 1% FDR)

Filtering - protein¶

In [14]:

mdata = mm.pp.add_filter(mdata, modality="protein", column="q_value", keep="lt", value=0.01)
mdata = mm.pp.apply_filter(mdata, modality="protein", on="var")

mdata

mdata = mm.pp.add_filter(mdata, modality="protein", column="q_value", keep="lt", value=0.01) mdata = mm.pp.apply_filter(mdata, modality="protein", on="var") mdata

INFO - Applying var filters for protein: ['q_value_lt_0.01']

Out[14]:

MuData object with n_obs × n_vars = 6 × 5847
  obs:	'source_name', 'characteristics_organism', 'characteristics_organism_part', 'characteristics_disease', 'characteristics_biological_replicate', 'assay_name', 'technology_type', 'comment_proteomexchange_accession_number', 'comment_proteomics_data_acquisition_method', 'comment_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', 'comment_instrument', 'comment_cleavage_agent_details', 'factor_value_condition', 'factor_value_time_point'
  uns:	'sdrf', '_cmd'
  3 modalities
    psm:	6 x 2264
      obs:	'source_name', 'characteristics_organism', 'characteristics_organism_part', 'characteristics_disease', 'characteristics_biological_replicate', 'assay_name', 'technology_type', 'comment_proteomexchange_accession_number', 'comment_proteomics_data_acquisition_method', 'comment_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', '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', 'sdrf', 'filter', 'decoy_filter'
      obsm:	'filter'
      varm:	'search_result', 'filter'
    peptide:	6 x 2151
      obs:	'source_name', 'characteristics_organism', 'characteristics_organism_part', 'characteristics_disease', 'characteristics_biological_replicate', 'assay_name', 'technology_type', 'comment_proteomexchange_accession_number', 'comment_proteomics_data_acquisition_method', 'comment_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', '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'
    protein:	6 x 1432
      obs:	'source_name', 'characteristics_organism', 'characteristics_organism_part', 'characteristics_disease', 'characteristics_biological_replicate', 'assay_name', 'technology_type', 'comment_proteomexchange_accession_number', 'comment_proteomics_data_acquisition_method', 'comment_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', 'comment_instrument', 'comment_cleavage_agent_details', 'factor_value_condition', 'factor_value_time_point'
      var:	'count_psm', 'count_stripped_peptide', 'PEP', 'q_value'
      uns:	'level', 'decoy', 'filter', 'decoy_filter'
      varm:	'filter'

Visualization¶

ID plot¶

In [15]:

mm.pl.plot_id(mdata, modality="protein")

mm.pl.plot_id(mdata, modality="protein")

Intensity distribution plot¶

In [16]:

mm.pl.plot_intensity(mdata, modality="protein", groupby="factor_value_condition")

mm.pl.plot_intensity(mdata, modality="protein", groupby="factor_value_condition")

Saving Processed MuData Object¶

You can save MuData object into an h5mu file. This allows you to easily reload the processed data in future sessions without repeating the entire analysis pipeline.

In [17]:

mdata.write_h5mu("dda_tmt_PXD013361.h5mu")

mdata.write_h5mu("dda_tmt_PXD013361.h5mu")

Citation¶

Magnusson, R., Rundquist, O., Kim, M. J., Hellberg, S., Na, C. H., Benson, M., ... & Gustafsson, M. (2019). A validated strategy to infer protein biomarkers from RNA-Seq by combining multiple mRNA splice variants and time-delay. BioRxiv, 599373.

Lazear, M. R. (2023). Sage: an open-source tool for fast proteomics searching and quantification at scale. Journal of Proteome Research, 22(11), 3652-3659.