01 Process proteomics data¶

In [1]:

import msmu as mm
from pathlib import Path
import requests
import tarfile

import msmu as mm from pathlib import Path import requests import tarfile

In [2]:

import plotly.io as pio

pio.renderers["notebook_connected"].global_init = False
pio.renderers.default = "plotly_mimetype+notebook_connected"

import plotly.io as pio pio.renderers["notebook_connected"].global_init = False pio.renderers.default = "plotly_mimetype+notebook_connected"

In [3]:

url = "https://github.com/bertis-informatics/msmu/releases/download/v0.4.3/MSV000089280.tar.gz"
sdrf = f"https://github.com/bertis-informatics/msmu/releases/download/v0.4.3/sdrf.tsv"
base_dir = Path(url).name.split(".")[0]
idents = f"{base_dir}/results.sage.tsv"
quants = f"{base_dir}/QuantifiedPeptides.tsv"

r = requests.get(url)
r.raise_for_status()

with open(Path(url).name, "wb") as f:
    f.write(r.content)

with tarfile.open(Path(url).name, "r:gz") as tar:
    members = [m for m in tar.getmembers() if not Path(m.name).name.startswith("._")]
    tar.extractall(members=members, filter="data")

url = "https://github.com/bertis-informatics/msmu/releases/download/v0.4.3/MSV000089280.tar.gz" sdrf = f"https://github.com/bertis-informatics/msmu/releases/download/v0.4.3/sdrf.tsv" base_dir = Path(url).name.split(".")[0] idents = f"{base_dir}/results.sage.tsv" quants = f"{base_dir}/QuantifiedPeptides.tsv" r = requests.get(url) r.raise_for_status() with open(Path(url).name, "wb") as f: f.write(r.content) with tarfile.open(Path(url).name, "r:gz") as tar: members = [m for m in tar.getmembers() if not Path(m.name).name.startswith("._")] tar.extractall(members=members, filter="data")

In [4]:

mdata = mm.read_sage(identification_file=idents, label="label_free")
mdata = mm.io.add_quant(mdata, quant_data=quants, quant_tool="flashlfq")

mdata = mm.read_sage(identification_file=idents, label="label_free") mdata = mm.io.add_quant(mdata, quant_data=quants, quant_tool="flashlfq")

INFO - Reading SAGE Identification data: 1 file(s)
INFO - Added quantification modality 'peptide' using flashlfq data.

In [5]:

mdata = mm.pp.attach_sdrf(mdata, sdrf=sdrf)
mdata = mm.pp.apply_sdrf_to_obs(mdata)

mdata = mm.pp.attach_sdrf(mdata, sdrf=sdrf) mdata = mm.pp.apply_sdrf_to_obs(mdata)

INFO - Validating SDRF metadata for https://github.com/bertis-informatics/msmu/releases/download/v0.4.3/sdrf.tsv.
INFO - SDRF validation succeeded for https://github.com/bertis-informatics/msmu/releases/download/v0.4.3/sdrf.tsv.

PSM-level processing¶

Filter PSMs¶

In [6]:

mdata = mm.pp.add_filter(mdata, modality="psm", column="q_value", keep="lt", value=0.01, on="var")
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, on="var") mdata = mm.pp.apply_filter(mdata, modality="psm") mdata

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

Out[6]:

MuData object with n_obs × n_vars = 70 × 85857
  obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[cell line]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomics data acquisition method]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'comment[modification parameters]', 'comment[modification parameters].1', 'comment[modification parameters].2', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[file uri]', 'factor value[sample]', 'factor value[condition]', 'factor value[rna data file]'
  uns:	'_log', 'sdrf'
  2 modalities
    psm:	70 × 80375
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[cell line]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomics data acquisition method]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'comment[modification parameters]', 'comment[modification parameters].1', 'comment[modification parameters].2', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[file uri]', 'factor value[sample]', 'factor value[condition]', 'factor value[rna data file]'
      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:	70 × 5482
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[cell line]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomics data acquisition method]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'comment[modification parameters]', 'comment[modification parameters].1', 'comment[modification parameters].2', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[file uri]', 'factor value[sample]', 'factor value[condition]', 'factor value[rna data file]'
      uns:	'level'
      layers:	None

Peptide-level processing¶

Summarize PSMs to peptides¶

In [7]:

mdata = mm.pp.to_peptide(mdata)
mdata

mdata = mm.pp.to_peptide(mdata) mdata

INFO - Peptide-level identifications: 5547 (3910 at 1% FDR)
INFO - Using existing peptide quantification data.

Out[7]:

MuData object with n_obs × n_vars = 70 × 80375
  obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[cell line]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomics data acquisition method]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'comment[modification parameters]', 'comment[modification parameters].1', 'comment[modification parameters].2', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[file uri]', 'factor value[sample]', 'factor value[condition]', 'factor value[rna data file]'
  uns:	'_log', 'sdrf'
  2 modalities
    psm:	70 × 80375
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[cell line]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomics data acquisition method]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'comment[modification parameters]', 'comment[modification parameters].1', 'comment[modification parameters].2', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[file uri]', 'factor value[sample]', 'factor value[condition]', 'factor value[rna data file]'
      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:	70 × 5547
      var:	'peptide', 'proteins', 'stripped_peptide', 'count_psm', 'PEP', 'q_value'
      uns:	'level', 'decoy'
      layers:	None

Filter peptide identifications¶

In [8]:

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

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") mdata

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

Out[8]:

MuData object with n_obs × n_vars = 70 × 84285
  obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[cell line]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomics data acquisition method]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'comment[modification parameters]', 'comment[modification parameters].1', 'comment[modification parameters].2', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[file uri]', 'factor value[sample]', 'factor value[condition]', 'factor value[rna data file]'
  uns:	'_log', 'sdrf'
  2 modalities
    psm:	70 × 80375
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[cell line]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomics data acquisition method]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'comment[modification parameters]', 'comment[modification parameters].1', 'comment[modification parameters].2', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[file uri]', 'factor value[sample]', 'factor value[condition]', 'factor value[rna data file]'
      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:	70 × 3910
      var:	'peptide', 'proteins', 'stripped_peptide', 'count_psm', 'PEP', 'q_value'
      uns:	'level', 'decoy', 'filter', 'decoy_filter'
      varm:	'filter'
      layers:	None

Transform and normalize peptide intensities¶

In [9]:

mdata = mm.dt.save_layer(mdata, modality="psm", layer="raw")

mdata = mm.pp.log2_transform(mdata, modality="peptide")
mdata = mm.pp.normalize(mdata, modality="peptide", method="pairwise_median")

mdata = mm.dt.save_layer(mdata, modality="psm", layer="raw") mdata = mm.pp.log2_transform(mdata, modality="peptide") mdata = mm.pp.normalize(mdata, modality="peptide", method="pairwise_median")

WARNING - normalise[pairwise_median] peptide block X|all|all: samples 'nanoSPLITs_C10_211005J_1cell_B2' and 'nanoSPLITs_SVEC_211007J_1cell_A12' share only 62 features (median pair shares 383); their pairwise median rests on few features. Check whether one of these runs failed.
INFO - normalise[pairwise_median] peptide layer X block all|all: 70 samples; observed features per sample 206..1916 (median 691); shift -1.577..+1.591 log2

Check normalized peptide intensities¶

In [10]:

peptide_intensity = mm.pl.plot_intensity(
    mdata, modality="peptide", ptype="box",
    colorby="factor value[condition]", width=900, height=500,
)
peptide_intensity.show()

peptide_intensity = mm.pl.plot_intensity( mdata, modality="peptide", ptype="box", colorby="factor value[condition]", width=900, height=500, ) peptide_intensity.show()

Protein inference¶

In [11]:

mdata = mm.pp.infer_protein(mdata)

mdata = mm.pp.infer_protein(mdata)

INFO - Starting protein inference
INFO - Initial proteins: 1862
INFO - Removed indistinguishable: 524
INFO - Removed subsettable: 200
INFO - Removed subsumable: 11
INFO - Total protein groups: 1127

Protein-level processing¶

Summarize peptides to proteins¶

In [12]:

mdata = mm.pp.to_protein(mdata, agg_method="median_polish")

mdata = mm.pp.to_protein(mdata, agg_method="median_polish")

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

Filter proteins¶

In [13]:

mdata = mm.pp.add_filter(mdata, modality="protein", column="q_value", keep="lt", value=0.01)
mdata = mm.pp.calculate_missingness(mdata, modality="protein")
mdata = mm.pp.add_filter(mdata, modality="protein", column="missing_ratio", keep="lt", value=1)
mdata = mm.pp.apply_filter(mdata, modality="protein")

mdata

mdata = mm.pp.add_filter(mdata, modality="protein", column="q_value", keep="lt", value=0.01) mdata = mm.pp.calculate_missingness(mdata, modality="protein") mdata = mm.pp.add_filter(mdata, modality="protein", column="missing_ratio", keep="lt", value=1) mdata = mm.pp.apply_filter(mdata, modality="protein") mdata

INFO - protein: 14 of 1089 features have no value in any sample
INFO - Applying var filters for protein: ['q_value_lt_0.01', 'missing_ratio_lt_1']

Out[13]:

MuData object with n_obs × n_vars = 70 × 85218
  obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[cell line]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomics data acquisition method]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'comment[modification parameters]', 'comment[modification parameters].1', 'comment[modification parameters].2', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[file uri]', 'factor value[sample]', 'factor value[condition]', 'factor value[rna data file]'
  uns:	'_log', 'sdrf'
  3 modalities
    psm:	70 × 80375
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[cell line]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomics data acquisition method]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'comment[modification parameters]', 'comment[modification parameters].1', 'comment[modification parameters].2', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[file uri]', 'factor value[sample]', 'factor value[condition]', 'factor value[rna data file]'
      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, 'raw'
    peptide:	70 × 3910
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[cell line]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomics data acquisition method]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'comment[modification parameters]', 'comment[modification parameters].1', 'comment[modification parameters].2', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[file uri]', 'factor value[sample]', 'factor value[condition]', 'factor value[rna data file]'
      var:	'peptide', 'proteins', 'stripped_peptide', 'count_psm', 'PEP', 'q_value', 'protein_group', 'peptide_type'
      uns:	'level', 'decoy', 'filter', 'decoy_filter', 'normalisation'
      varm:	'filter'
      layers:	None
    protein:	70 × 933
      obs:	'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[disease]', 'characteristics[cell line]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomics data acquisition method]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'comment[modification parameters]', 'comment[modification parameters].1', 'comment[modification parameters].2', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[file uri]', 'factor value[sample]', 'factor value[condition]', 'factor value[rna data file]'
      var:	'count_psm', 'count_stripped_peptide', 'PEP', 'q_value', 'n_missing', 'missing_ratio'
      uns:	'level', 'decoy', 'filter', 'decoy_filter'
      varm:	'filter'
      layers:	None

Check protein identifications¶

In [14]:

protein_ids = mm.pl.plot_id(
    mdata, modality="protein",
    colorby="factor value[condition]", width=900, height=500,
)
protein_ids.show()

protein_ids = mm.pl.plot_id( mdata, modality="protein", colorby="factor value[condition]", width=900, height=500, ) protein_ids.show()

Inspect protein overlap¶

In [15]:

protein_overlap = mm.pl.plot_upset(
    mdata, modality="protein", groupby="factor value[condition]",
    width=700, height=450,
)
protein_overlap.show()

protein_overlap = mm.pl.plot_upset( mdata, modality="protein", groupby="factor value[condition]", width=700, height=450, ) protein_overlap.show()

Inspect protein correlations¶

In [16]:

protein_correlation = mm.pl.plot_correlation(
    mdata, modality="protein", width=800, height=700,
)
protein_correlation.show()

protein_correlation = mm.pl.plot_correlation( mdata, modality="protein", width=800, height=700, ) protein_correlation.show()

Save processed data¶

In [17]:

mdata.write_h5mu("MSV000089280.h5mu")

mdata.write_h5mu("MSV000089280.h5mu")

Inspect differential expression¶

In [18]:

de = mm.tl.run_de(
    mdata, modality="protein", category="factor value[condition]",
    ctrl="C10", expr="SVEC",
)
volcano = mm.pl.plot_volcano(
    de.to_df(), ctrl=de.ctrl, expr=de.expr, log2fc_threshold=de.fc_pct_5,
)
volcano.show()

de = mm.tl.run_de( mdata, modality="protein", category="factor value[condition]", ctrl="C10", expr="SVEC", ) volcano = mm.pl.plot_volcano( de.to_df(), ctrl=de.ctrl, expr=de.expr, log2fc_threshold=de.fc_pct_5, ) volcano.show()

WARNING - mm.tl.run_de: the default DE engine is now 'limma' (previously a permutation test). Pass stat_method explicitly to choose an engine. This notice will be removed in a future release.