Tissue-specific proteoform inference across five mouse organs

This notebook shows how to use ProteoPy for proteoform group inference following a previously established approach that we published in Bludau et al. (2021):

Bludau I, Frank M, Dörig C, Cai Y, Heusel M, Rosenberger G, Picotti P, Collins BC, Röst H, and Aebersold R. Systematic detection of functional proteoform groups from bottom-up proteomic datasets. Nature Communications, 12:3810, 2021. doi:10.1038/s41467-021-24030-x

In this paper, we introduced COPF (COrrelation-based functional ProteoForm assessment), a data-driven strategy that assigns peptides to co-varying proteoform groups based on peptide correlation patterns. It was applied to SWATH-MS (DIA) data from five mouse tissues (brain, brown adipose tissue, heart, liver, and quadriceps) across eight BXD mice, originally published by Williams et al. (2018):

Williams EG, Wu Y, Wolski W, Kim JY, Lan J, Hasan M, Halter C, Jha P, Ryu D, Auwerx J, and Aebersold R. Quantifying and Localizing the Mitochondrial Proteome Across Five Tissues in A Mouse Population. Molecular & Cellular Proteomics, 17(9):1766–1777, 2018. doi:10.1074/mcp.RA118.000554

In this notebook, we show how the COPF implementation in ProteoPy can be used to perform the core proteoform inference workflow.

What this notebook reproduces

This is a faithful reproduction, not an approximation. Running it end to end recovers the same 63 proteins that Bludau et al. report as carrying multiple proteoform groups — the same count and the same accessions, including the two the paper singles out, Ldb3 (Q9JKS4) and Sorbs2 (Q3UTJ2) — together with the proteoform score plot of Figure 6B.

ProteoPy’s COPF implementation is itself tested against the intermediate outputs of the original R code on this very dataset. Given the same peptides, its correlations, dendrograms, cluster assignments and proteoform scores agree with the reference to between 1e-14 and 1e-12, which leaves the peptide set as the one thing this notebook has to get right.

Setup

[1]:

import random
from pathlib import Path
import numpy as np
import scanpy as sc
import matplotlib as mpl
from matplotlib.pyplot import rc_context

import proteopy as pr  # Convention: import proteopy as pr
from proteopy.utils import is_proteodata

# Set random seed for reproducibility
random.seed(42)

# Create a data directory in your current working directory
# to store files downloaded in this notebook.
cwd = Path(".").resolve()
(cwd / "data").mkdir(parents=True, exist_ok=True)

# Protein sequences, needed to locate each peptide in its protein.
# Shipped alongside this notebook so it runs straight after cloning.
FASTA_PATH = "uniprot_mus-musculus_reviewed_bludau-2021.fasta"

/home/ifichtner/miniforge3/envs/proteoform-repro_exp-04/lib/python3.10/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Reading in the data

The Williams et al. (2018) dataset contains peptide-level SWATH-MS intensities from five tissues of eight BXD mice. ProteoPy provides a built-in download function (pr.download.williams_2018()) that fetches the dataset from the original supplementary archive, processes it, and saves it as three separate files:

  • Intensities — long-format table with columns sample_id, peptide_id, and intensity (one row per measurement)

  • Sample annotation — maps each sample_id to its tissue and mouse_id

  • Peptide annotation — maps each peptide_id to its protein_id and gene_id

This separation into three files mirrors a common data delivery format in proteomics. The files are then read into an AnnData object using pr.read.long().

[2]:

# Define paths to the data files. These are the same paths that will be used in the download function below.
intensities_path = "data/williams-2018_mouse-tissue_intensities.tsv"
sample_annotation_path = (
    "data/williams-2018_mouse-tissue_sample_annotation.tsv"
)
peptide_annotation_path = (
    "data/williams-2018_mouse-tissue_peptide_annotation.tsv"
)

pr.download.williams_2018(
    intensities_path=intensities_path,
    sample_annotation_path=sample_annotation_path,
    var_annotation_path=peptide_annotation_path,
    force=True,  # Overwrite files if already exist
)

A quick look at the first rows of each file shows the data format ProteoPy expects for each information layer:

[3]:

# intensities
!head -n5 {intensities_path}

sample_id       peptide_id      intensity
C57_Brain       AAAAAAAAAAAAAAAGAAGK    26302.0
DBA_Brain       AAAAAAAAAAAAAAAGAAGK    15460.0
45_Brain        AAAAAAAAAAAAAAAGAAGK    24631.0
66_Brain        AAAAAAAAAAAAAAAGAAGK    26554.0

[4]:

# sample_annotation
!head -n5 {sample_annotation_path}

sample_id       tissue  mouse_id
C57_Brain       Brain   C57
DBA_Brain       Brain   DBA
45_Brain        Brain   45
66_Brain        Brain   66

[5]:

# var_annotation
!head -n5 {peptide_annotation_path}

peptide_id      protein_id      gene_id
AAAAAAAAAAAAAAAGAAGK    P55012  Slc12a2
AAAAADLANR      Q9JHS4  Clpx
AAAADGEPLHNEEER Q80WW9  Ddrgk1
AAAAEGARPLER    Q80UM7  Mogs

With the three files on disk, pr.read.long() reads the long-format intensities, pivots them into a dense matrix, and merges the sample and peptide annotations into a single AnnData object.

[6]:

# Create a peptide-level AnnData object from the downloaded files.
# `fill_na` is deliberately left unset: see above.
adata = pr.read.long(
    intensities=intensities_path,
    sample_annotation=sample_annotation_path,
    var_annotation=peptide_annotation_path,
    level="peptide",
)

adata

[6]:

AnnData object with n_obs × n_vars = 40 × 32690
    obs: 'sample_id', 'tissue', 'mouse_id'
    var: 'peptide_id', 'protein_id', 'protein_id_annotation', 'gene_id'

[7]:

print("Is proteodata: ", is_proteodata(adata))

Is proteodata:  (True, 'peptide')

[8]:

adata.obs.head(n=10)

[8]:
sample_id tissue mouse_id
101_Brain 101_Brain Brain 101
45_Brain 45_Brain Brain 45
66_Brain 66_Brain Brain 66
68_Brain 68_Brain Brain 68
73_Brain 73_Brain Brain 73
80_Brain 80_Brain Brain 80
BAT_101 BAT_101 BAT 101
BAT_45 BAT_45 BAT 45
BAT_66 BAT_66 BAT 66
BAT_68 BAT_68 BAT 68 
[9]:

adata.var.head(n=10)

[9]:
peptide_id protein_id protein_id_annotation gene_id
AAAAAAAAAAAAAAAGAAGK AAAAAAAAAAAAAAAGAAGK P55012 P55012 Slc12a2
AAAAADLANR AAAAADLANR Q9JHS4 Q9JHS4 Clpx
AAAADGEPLHNEEER AAAADGEPLHNEEER Q80WW9 Q80WW9 Ddrgk1
AAAAEGARPLER AAAAEGARPLER Q80UM7 Q80UM7 Mogs
AAAAKEEAPK AAAAKEEAPK Q91XV3 Q91XV3 Basp1
AAAANLC(UniMod:4)PGDVILAIDGFGTESMTHADAQDR AAAANLC(UniMod:4)PGDVILAIDGFGTESMTHADAQDR O70209 O70209 Pdlim3
AAAAYALGR AAAAYALGR Q6ZQ73 Q6ZQ73 Cand2
AAADLM(UniMod:35)AYC(UniMod:4)EAHAKEDPLLTPVPASENPFR AAADLM(UniMod:35)AYC(UniMod:4)EAHAKEDPLLTPVPAS... P63213 P63213 Gng2
AAADLMAYC(UniMod:4)EAHAK AAADLMAYC(UniMod:4)EAHAK P63213 P63213 Gng2
AAADLMAYC(UniMod:4)EAHAKEDPLLTPVPASENPFR AAADLMAYC(UniMod:4)EAHAKEDPLLTPVPASENPFR P63213 P63213 Gng2

You can define consistent color schemes and category orders for the different tissues and BXD mice used for all figures created by ProteoPy.

[10]:

# Tissues
n_tissues = adata.obs["tissue"].nunique()
cmap = mpl.colormaps["Set2"]
adata.uns["colors_tissue"] = cmap(range(n_tissues)).tolist()

adata.uns["order_tissue"] = ["Brain", "BAT", "Heart", "Liver", "Quad"]

# BXD mice
n_mice = adata.obs["mouse_id"].nunique()
cmap = mpl.colormaps["tab10"]
adata.uns["colors_mouse_id"] = cmap(range(n_mice)).tolist()

adata.uns["order_mouse_id"] = sorted(adata.obs["mouse_id"].unique().tolist())

Looking at the sample distribution shows successful data import in agreement with the study design. The dataset includes 5 tissues derived from 8 BXD mice:

[11]:

pr.pl.n_samples_per_category(
    adata,
    category_key="tissue",
    order=adata.uns["order_tissue"],
    color_scheme=adata.uns["colors_tissue"],
)
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_17_0.png 
[12]:

pr.pl.n_samples_per_category(
    adata,
    category_key="mouse_id",
    order=adata.uns["order_mouse_id"],
    color_scheme=adata.uns["colors_mouse_id"],
)
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_18_0.png

Quality control and preprocessing

As per the original analysis, we remove one protein by hand: A2ASS6 (titin), which carries far more peptides than any other and distorts the analysis.

[13]:

A2ASS6_mask = adata.var["protein_id"] == "A2ASS6"
print(f"N peptides for protein A2ASS6: {A2ASS6_mask.sum()}")

N peptides for protein A2ASS6: 919

[14]:

adata = adata[:, ~A2ASS6_mask]

Summarising redundant peptides

Several peptides can report on the same stretch of a protein: missed cleavages produce overlapping sequences, and one sequence can appear in several modified forms introduced during sample preparation or MS analysis. We therefore group peptides by where they sit in the protein sequence and keep the most abundant member of each group.

[15]:

pr.pp.summarize_peptides_by_neighbourhood_union(
    adata,
    annotator=FASTA_PATH,
    top_n=1,
    on_unknown_protein="keep",
    on_unlocated_peptide="keep",
    keep_var_cols=["gene_id"],
    verbose=True,
)

summarize_peptides_by_neighbourhood_union: 31771 -> 25613 peptides across 25613 neighbourhood group(s)

Filtering

In preparation for the COPF algorithm we filter the dataset on three criteria:

  1. Peptides must be quantified in every sample.

  2. Peptides must have a positive variance.

  3. A protein needs at least 2 peptides to be scorable.

[16]:

pr.pp.filter_var_completeness(adata, min_fraction=1.0)
pr.pp.remove_zero_variance_vars(adata)
pr.pp.filter_proteins_by_peptide_count(adata, min_count=2)

print(
    f"{adata.n_vars:,} peptides across "
    f"{adata.var['protein_id'].nunique():,} proteins"
)

137 var removed
Removed 928 proteins and 928 peptides.
24,534 peptides across 2,885 proteins

This leaves 24,534 peptides across 2,885 proteins. Every peptide is now quantified in all 40 samples, so the exploratory analysis below needs no imputation.

[17]:

pr.pl.n_peptides_per_sample(
    adata,
    zero_to_na=True,
    order_by="tissue",
    order=adata.uns["order_tissue"],
    color_scheme=adata.uns["colors_tissue"],
    print_stats=True,
)

Global:
 mean_count  std_count  median_count  min_count  max_count  mean_pct  std_pct  median_pct  min_pct  max_pct
    24322.5      158.4       24405.5      23828      24471      99.1      0.6        99.5     97.1     99.7

Per tissue:
tissue  mean_count  std_count  median_count  min_count  max_count  mean_pct  std_pct  median_pct  min_pct  max_pct
   BAT     24317.8      125.0       24259.5      24188      24471      99.1      0.5        98.9     98.6     99.7
 Brain     24425.5        5.2       24424.5      24420      24435      99.6      0.0        99.6     99.5     99.6
 Heart     24053.1      104.8       24050.0      23828      24178      98.0      0.4        98.0     97.1     98.5
 Liver     24424.9       10.4       24427.0      24409      24441      99.6      0.0        99.6     99.5     99.6
  Quad     24391.4       23.4       24389.0      24368      24439      99.4      0.1        99.4     99.3     99.6
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_28_1.png 
[17]:

<Axes: ylabel='#'>

[18]:

pr.pl.n_proteins_per_sample(
    adata,
    zero_to_na=True,
    order_by="tissue",
    order=adata.uns["order_tissue"],
    color_scheme=adata.uns["colors_tissue"],
    print_stats=True,
)

Global:
 mean_count  std_count  median_count  min_count  max_count  mean_pct  std_pct  median_pct  min_pct  max_pct
     2884.8        0.4        2885.0       2883       2885     100.0      0.0       100.0     99.9    100.0

Per tissue:
tissue  mean_count  std_count  median_count  min_count  max_count  mean_pct  std_pct  median_pct  min_pct  max_pct
   BAT      2884.9        0.4        2885.0       2884       2885     100.0      0.0       100.0    100.0    100.0
 Brain      2885.0        0.0        2885.0       2885       2885     100.0      0.0       100.0    100.0    100.0
 Heart      2884.4        0.7        2884.5       2883       2885     100.0      0.0       100.0     99.9    100.0
 Liver      2885.0        0.0        2885.0       2885       2885     100.0      0.0       100.0    100.0    100.0
  Quad      2885.0        0.0        2885.0       2885       2885     100.0      0.0       100.0    100.0    100.0
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_29_1.png 
[18]:

<Axes: ylabel='#'>

Since proteoform group inference depends on high peptide coverage, we also inspect the distribution of peptides per protein:

[19]:

pr.pl.n_peptides_per_protein(adata, print_stats=True)

    mean  median  mode  variance  min  max
8.503986     6.0     2 80.414426    2  117
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_31_1.png 
[19]:

<Axes: xlabel='Number of peptide_id per protein_id', ylabel='# protein_id'>

Note that no log transformation, data normalization or missing value imputation is performed here. This is because of the overall very high data consistency across the dataset and because COPF is based on covariation rather than direct intensity differences across conditions. This is also in agreement with the analysis by Bludau et al. (2021). However, other datasets might require additional preprocessing steps.

Exploratory analysis of peptide-level data

Before proteoform inference, ProteoPy can be used to perform classical exploratory data analysis steps, for example to investigate the primary sources of variation in the dataset.

Sample correlation matrix

High intra-tissue correlation and lower inter-tissue correlation show that peptide intensities are strongly driven by the tissue of origin and not by the mouse strain.

[20]:

pr.pl.sample_correlation_matrix(
    adata,
    method="pearson",
    margin_color="tissue",
    color_scheme=adata.uns["colors_tissue"],
)
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_35_0.png

Dimensionality reduction

The AnnData-based proteodata object created and used by ProteoPy is directly compatible with other AnnData based python packages such as scanpy (import scanpy as sc), which can for example be used for dimensionality reduction and clustering.

Here, the PCA and UMAP analyses further confirm that tissue origin is the dominant source of variation in the dataset.

PCA analysis

[21]:

sc.tl.pca(adata)

with rc_context({"figure.figsize": (5, 3)}):
    sc.pl.pca_variance_ratio(adata, n_pcs=50, log=True)
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_38_0.png 
[22]:

sc.pl.pca(
    adata,
    color=["tissue", "tissue", "tissue"],
    dimensions=[(0, 1), (1, 2), (2, 3)],
    ncols=3,
    size=90,
    palette=adata.uns["colors_tissue"],
)

sc.pl.pca(
    adata,
    color=["mouse_id", "mouse_id", "mouse_id"],
    dimensions=[(0, 1), (1, 2), (2, 3)],
    ncols=3,
    size=90,
    palette=adata.uns["colors_mouse_id"],
)
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_39_0.png
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_39_1.png

UMAP

[23]:

sc.pp.neighbors(adata, n_neighbors=4)
sc.tl.umap(adata)

sc.pl.umap(
    adata,
    color=["tissue"],
    size=100,
    palette=adata.uns["colors_tissue"],
)

sc.pl.umap(
    adata,
    color=["mouse_id"],
    size=100,
    palette=adata.uns["colors_mouse_id"],
)

/home/ifichtner/miniforge3/envs/proteoform-repro_exp-04/lib/python3.10/site-packages/sklearn/manifold/_spectral_embedding.py:328: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
  warnings.warn(
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_41_1.png
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_41_2.png

Proteoform inference

ProteoPy implements the COPF workflow’s four main steps:

  1. compute pairwise correlations between all peptides of a protein

  2. cluster the correlation distances hierarchically

  3. cut each dendrogram into two candidate proteoform groups

  4. score the difference between within- and between-cluster correlation, and compute Benjamini–Hochberg adjusted p-value.

These four steps are covered by ProteoPy’s test suite against the original R implementation’s own intermediate outputs on this dataset, agreeing to between 1e-14 and 1e-12.

[24]:

# COPF pipeline: correlate, cluster, score
pr.tl.pairwise_peptide_correlations(adata)
pr.tl.peptide_dendograms_by_correlation(
    adata, method="agglomerative-hierarchical-clustering"
)
pr.tl.peptide_clusters_from_dendograms(
    adata, n_clusters=2, min_peptides_per_cluster=2
)
pr.tl.proteoform_scores(adata, min_score=0.1, min_pval_adj=0.1)

The pseudo-volcano plot shows proteoform scores vs. adjusted p-values (related to Figure 6B in Bludau et al. (2021)). The two highlighted proteins — Ldb3 (Q9JKS4) and Sorbs2 (Q3UTJ2) — are both called significant, and are examined in detail below.

[25]:

# The two proteins the paper singles out, highlighted by gene symbol below.
bludauI_proteoforms = ["Ldb3", "Sorbs2"]

adata.var.loc[
    adata.var["gene_id"].isin(bludauI_proteoforms),
    ["protein_id", "gene_id"],
].drop_duplicates()

[25]:
protein_id gene_id
VLPGPSQPR Q9JKS4 Ldb3
VGIFPISYVEK Q3UTJ2 Sorbs2 
[26]:

pr.pl.proteoform_scores(
    adata,
    adj=True,
    pval_threshold=0.1,
    score_threshold=0.1,
    highlight_prots=bludauI_proteoforms,
    protein_id_key="gene_id",
)

Looks like you are using a tranform that doesn't support FancyArrowPatch, using ax.annotate instead. The arrows might strike through texts. Increasing shrinkA in arrowprops might help.
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_46_1.png 
[26]:

<Axes: xlabel='Proteoform Score', ylabel='adj. p-value'>

[27]:

pf_df = pr.get.proteoforms_df(adata, only_proteins=True)
sig_df = pf_df[pf_df["is_proteoform"] == 1]
n_proteoforms = sig_df.shape[0]
print(f"{n_proteoforms} significant proteoforms inferred.")

63 significant proteoforms inferred.

[ ]:

# The 63 accessions Bludau et al. (2021) report for this dataset, taken
# from a run of the authors' original R workflow on the same input
# https://github.com/ibludau/ProteoformAnanlysis/blob/main/MouseTissue/GetMouseTissueProteoforms_paper.R.
PUBLISHED_63 = {
    "O08601",
    "O54724",
    "O54749",
    "O54754",
    "O88844",
    "O88986",
    "P02088",
    "P07356",
    "P08226",
    "P09671",
    "P10852",
    "P13541",
    "P16125",
    "P16332",
    "P19096",
    "P37040",
    "P41216",
    "P50431",
    "P51660",
    "P58281",
    "P63030",
    "P63318",
    "Q00623",
    "Q04447",
    "Q3THK7",
    "Q3UTJ2",
    "Q5IRJ6",
    "Q5RKZ7",
    "Q60598",
    "Q61207",
    "Q62261",
    "Q71LX4",
    "Q7TNG8",
    "Q80X90",
    "Q8BMS1",
    "Q8BWF0",
    "Q8BYI9",
    "Q8C165",
    "Q8CI94",
    "Q8VDD5",
    "Q8VEM8",
    "Q8VIJ6",
    "Q91VN4",
    "Q91WD5",
    "Q924X2",
    "Q99J39",
    "Q99JB8",
    "Q99KK7",
    "Q99LC5",
    "Q9CPQ1",
    "Q9CQR4",
    "Q9D0K2",
    "Q9D5T0",
    "Q9DBH5",
    "Q9DBM2",
    "Q9ES97",
    "Q9ET01",
    "Q9JI39",
    "Q9JKS4",
    "Q9QYR6",
    "Q9QZ47",
    "Q9Z204",
    "Q9Z2Z6",
}

significant = set(sig_df["protein_id"])

print(f"peptides analysed         : {adata.n_vars:,}")
print(f"proteins analysed         : {adata.var['protein_id'].nunique():,}")
print(f"proteins scored (BH n)    : {len(pf_df):,}")
print()
print(
    f"significant proteins      : {len(significant)}   (published: {len(PUBLISHED_63)})"
)
print(f"recovered of the published: {len(significant & PUBLISHED_63)}")
print(f"missed                    : {len(PUBLISHED_63 - significant)}")
print(f"false positives           : {len(significant - PUBLISHED_63)}")
print(f"IDENTICAL SET             : {significant == PUBLISHED_63}")
print()
print(f"Sorbs2 (Q3UTJ2) recovered : {'Q3UTJ2' in significant}")
print(f"Ldb3   (Q9JKS4) recovered : {'Q9JKS4' in significant}")

peptides analysed         : 24,534
proteins analysed         : 2,885
proteins scored (BH n)    : 1,272

significant proteins      : 63   (published: 63)
recovered of the published: 63
missed                    : 0
false positives           : 0
IDENTICAL SET             : True

Sorbs2 (Q3UTJ2) recovered : True
Ldb3   (Q9JKS4) recovered : True

Cleaning up before the downstream sections

COPF marks peptides that fit neither cluster as outliers. They played their part in the scoring above and are removed here so the per-proteoform sections that follow work with cluster members only. Note this also drops the proteins that could not be scored at all, so the object is smaller from this point on.

[29]:

# Remove COPF outlier peptides (cluster_id == 1000000)
copf_outliers = adata.var["cluster_id"] == 1000000
adata = adata[:, ~copf_outliers]

Proteoform exploration

ProteoPy also recovered the two exemplary proteins with proteoform groups highlighted in Bludau et al (2021): LIM domain-binding protein 3 (Ldb3, Q9JKS4) and Sorbin and SH3 domain-containing protein 2 (Sorbs2, Q3UTJ2).

[30]:

# Summarisation keeps each surviving peptide's own identifier, modifications
# included, as the reference does. The sequence maps below need the bare
# amino-acid sequence to locate peptides on the protein, so derive it into a
# separate column rather than rewriting the identifiers.
adata.var["peptide_sequence"] = adata.var["peptide_id"].str.replace(
    r"\(UniMod:\d+\)", "", regex=True
)

/tmp/ipykernel_3009267/1804092523.py:5: ImplicitModificationWarning: Trying to modify attribute `.var` of view, initializing view as actual.
  adata.var["peptide_sequence"] = adata.var["peptide_id"].str.replace(

LIM domain-binding protein 3 — Ldb3 (Q9JKS4)

Ldb3 is a muscle-specific protein. COPF assigns its peptides to two tissue-specific proteoform groups, consistent with known alternative splice variants (reproduces Figure 7A).

[31]:

pr.pl.proteoform_intensities(
    adata,
    protein_ids="Q9JKS4",
    order_by="tissue",
    order=adata.uns["order_tissue"],
    xlab_rotation=45,
    color_scheme=["#ff7f0e", "#1f77b4"],
)
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_54_0.png 
[32]:

pr.get.proteoforms_df(adata, proteins="Q9JKS4")

[32]:
protein_id peptide_id cluster_id proteoform_score proteoform_score_pval proteoform_score_pval_adj is_proteoform
0 Q9JKS4 VLPGPSQPR 0.0 0.56289 0.00166 0.055571 1.0
1 Q9JKS4 QYNNPIGLYSAETLR 0.0 0.56289 0.00166 0.055571 1.0
2 Q9JKS4 EMAQMYQMSLR 0.0 0.56289 0.00166 0.055571 1.0
3 Q9JKS4 ASSEGAQGSVSPK 0.0 0.56289 0.00166 0.055571 1.0
4 Q9JKS4 ASGAGLLGGSLPVK 0.0 0.56289 0.00166 0.055571 1.0
5 Q9JKS4 TSLADVC(UniMod:4)FVEEQNNVYC(UniMod:4)ER 1.0 0.56289 0.00166 0.055571 1.0
6 Q9JKS4 TQSKPEDEADEWAR 1.0 0.56289 0.00166 0.055571 1.0
7 Q9JKS4 TPLC(UniMod:4)GHC(UniMod:4)NNVIR 1.0 0.56289 0.00166 0.055571 1.0
8 Q9JKS4 SWHPEEFNC(UniMod:4)AYC(UniMod:4)K 1.0 0.56289 0.00166 0.055571 1.0
9 Q9JKS4 SKRPIPISTTAPPIQSPLPVIPHQK 1.0 0.56289 0.00166 0.055571 1.0
10 Q9JKS4 SASYNLSLTLQK 1.0 0.56289 0.00166 0.055571 1.0
11 Q9JKS4 QTWHTTC(UniMod:4)FVC(UniMod:4)AAC(UniMod:4)K 1.0 0.56289 0.00166 0.055571 1.0
12 Q9JKS4 LQGGKDFNMPLTISR 1.0 0.56289 0.00166 0.055571 1.0
13 Q9JKS4 IMGEVMHALR 1.0 0.56289 0.00166 0.055571 1.0
14 Q9JKS4 ILAQMTGTEYMQDPDEEALRR 1.0 0.56289 0.00166 0.055571 1.0
15 Q9JKS4 GPFLVAMGR 1.0 0.56289 0.00166 0.055571 1.0
16 Q9JKS4 GAPAYNPTGPQVTPLAR 1.0 0.56289 0.00166 0.055571 1.0
17 Q9JKS4 C(UniMod:4)YEQFFAPIC(UniMod:4)AK 1.0 0.56289 0.00166 0.055571 1.0
18 Q9JKS4 AAQSQLSQGDLVVAIDGVNTDTMTHLEAQNK 1.0 0.56289 0.00166 0.055571 1.0

Ldb3 peptide sequence map (reproduces Figure 7C)

Mapping detected peptides onto the canonical protein sequence and annotated UniProt alternative sequences confirms that the two proteoform groups correspond to known splice variants.

[33]:

# Canonical sequence and UniProt alternative sequences for Q9JKS4
# Source: https://www.uniprot.org/uniprotkb/Q9JKS4/entry#sequences
seq = "MSYSVTLTGPGPWGFRLQGGKDFNMPLTISRITPGSKAAQSQLSQGDLVVAIDGVNTDTMTHLEAQNKIKSASYNLSLTLQKSKRPIPISTTAPPIQSPLPVIPHQKDPALDTNGSLATPSPSPEARASPGALEFGDTFSSSFSQTSVCSPLMEASGPVLPLGSPVAKASSEGAQGSVSPKVLPGPSQPRQYNNPIGLYSAETLREMAQMYQMSLRGKASGAGLLGGSLPVKDLAVDSASPVYQAVIKTQSKPEDEADEWARRSSNLQSRSFRILAQMTGTEYMQDPDEEALRRSSTPIEHAPVCTSQATSPLLPASAQSPAAASPIAASPTLATAAATHAAAASAAGPAASPVENPRPQASAYSPAAAASPAPSAHTSYSEGPAAPAPKPRVVTTASIRPSVYQPVPASSYSPSPGANYSPTPYTPSPAPAYTPSPAPTYTPSPAPTYSPSPAPAYTPSPAPNYTPTPSAAYSGGPSESASRPPWVTDDSFSQKFAPGKSTTTVSKQTLPRGAPAYNPTGPQVTPLARGTFQRAERFPASSRTPLCGHCNNVIRGPFLVAMGRSWHPEEFNCAYCKTSLADVCFVEEQNNVYCERCYEQFFAPICAKCNTKIMGEVMHALRQTWHTTCFVCAACKKPFGNSLFHMEDGEPYCEKDYINLFSTKCHGCDFPVEAGDKFIEALGHTWHDTCFICAVCHVNLEGQPFYSKKDKPLCKKHAHAINV"
alt_seqs = {
    "Q9JKS4-2_0": {"seq_coord": (107, 227), "group": "alt_2"},
    "Q9JKS4-3_0": {"seq_coord": (295, 357), "group": "alt_3"},
    "Q9JKS4-4_0": {"seq_coord": (107, 227), "group": "alt_4"},
    "Q9JKS4-4_1": {"seq_coord": (295, 357), "group": "alt_4"},
    "Q9JKS4-5_0": {"seq_coord": (295, 327), "group": "alt_5"},
    "Q9JKS4-5_1": {"seq_coord": (328, 723), "group": "alt_5"},
    "Q9JKS4-6_0": {"seq_coord": (107, 227), "group": "alt_6"},
    "Q9JKS4-6_1": {"seq_coord": (295, 327), "group": "alt_6"},
    "Q9JKS4-6_2": {"seq_coord": (328, 723), "group": "alt_6"},
}

pr.pl.peptides_on_prot_sequence(
    adata,
    protein_id="Q9JKS4",
    group_by="proteoform_id",
    alt_pep_sequence_key="peptide_sequence",
    ref_sequence=seq,
    add_sequences=alt_seqs,
    title="LIM domain binding protein 3 (Q9JKS4)",
    figsize=(8, 4),
)
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_57_0.png 
[33]:

<Axes: title={'center': 'LIM domain binding protein 3 (Q9JKS4)'}, xlabel='Position'>

Sorbin and SH3 domain-containing protein 2 — Sorbs2 (Q3UTJ2)

Sorbs2 peptides are assigned to two proteoform groups with distinct tissue expression patterns: one group is abundant across brain, heart, and liver, while the other is brain-specific (reproduces Figure 7D).

[34]:

pr.pl.proteoform_intensities(
    adata,
    protein_ids="Q3UTJ2",
    order_by="tissue",
    order=adata.uns["order_tissue"],
    xlab_rotation=45,
)
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_59_0.png 
[35]:

pr.get.proteoforms_df(adata, proteins="Q3UTJ2")

[35]:
protein_id peptide_id cluster_id proteoform_score proteoform_score_pval proteoform_score_pval_adj is_proteoform
0 Q3UTJ2 LAFLVSPVPFR 0.0 0.619637 0.000343 0.015562 1.0
1 Q3UTJ2 ASVVEALDSALKDIC(UniMod:4)DQIK 0.0 0.619637 0.000343 0.015562 1.0
2 Q3UTJ2 VGIFPISYVEK 1.0 0.619637 0.000343 0.015562 1.0
3 Q3UTJ2 SYSSTLTDLGR 1.0 0.619637 0.000343 0.015562 1.0
4 Q3UTJ2 SIFEYEPGK 1.0 0.619637 0.000343 0.015562 1.0
5 Q3UTJ2 SFISSSPSSPSR 1.0 0.619637 0.000343 0.015562 1.0
6 Q3UTJ2 QGIFPVSYVEVVKR 1.0 0.619637 0.000343 0.015562 1.0
7 Q3UTJ2 GLGDQSSSR 1.0 0.619637 0.000343 0.015562 1.0
8 Q3UTJ2 AQPARPPPPVQPGEIGEAIAK 1.0 0.619637 0.000343 0.015562 1.0
9 Q3UTJ2 APHYPGIGPVDESGIPTAIR 1.0 0.619637 0.000343 0.015562 1.0
10 Q3UTJ2 ADLPGSSSTFTK 1.0 0.619637 0.000343 0.015562 1.0

Sorbs2 peptide sequence map (reproduces Figure 7F)

The brain-specific proteoform group maps to a region covered by alternative sequences 3, 4 and 5 in UniProt (10 in Bludau et. al. Figure 7F). The blue peptides assigned to proteoform group 0 map to a region previously described as brain-specific and exclusive to neurons.

[36]:

# Canonical sequence and UniProt alternative sequences for Q3UTJ2
# Source: https://www.uniprot.org/uniprotkb/Q3UTJ2/entry#sequences
seq = "MNTDSGGCARKRAAMSVTLTSVKRVQSSPNLLAAGRESQSPDSAWRSYNDRNPETLNGDATYSSLAAKGFRSVRPNLQDKRSPTQSQITINGNSGGAVSPVSYYQRPFSPSAYSLPASLNSSIIMQHGRSLDSAETYSQHAQSLDGTMGSSIPLYRSSEEEKRVTVIKAPHYPGIGPVDESGIPTAIRTTVDRPKDWYKTMFKQIHMVHKPGLYNSPYSAQSHPAAKTQTYRPLSKSHSDNGTDAFKEVPSPVPPPHVPPRPRDQSSTLKHDWDPPDRKVDTRKFRSEPRSIFEYEPGKSSILQHERPVSIYQSSIDRSLERPSSSASMAGDFRKRRKSEPAVGPLRGLGDQSSSRTSPGRADLPGSSSTFTKSFISSSPSSPSRAQGGDDSKMCPPLCSYSGLNGTPSGELECCNAYRQHLDVPGDSQRAITFKNGWQMARQNAEIWSSTEETVSPKIKSRSCDDLLNDDCDSFPDPKTKSESMGSLLCEEDSKESCPMTWASPYIQEVCGNSRSRLKHRSAHNAPGFLKMYKKMHRINRKDLMNSEVICSVKSRILQYEKEQQHRGLLHGWSQSSTEEVPRDVVPTRISEFEKLIQKSKSMPNLGDEMLSPITLEPPQNGLCPKRRFSIESLLEEETQVRHPSQGQRSCKSNTLVPIHIEVTSDEQPRTHMEFSDSDQDGVVSDHSDYVHVEGSSFCSESDFDHFSFTSSESFYGSSHHHHHHHHHHRHLISSCKGRCPASYTRFTTMLKHERAKHENMDRPRRQEMDPGLSKLAFLVSPVPFRRKKILTPQKQTEKAKCKASVVEALDSALKDICDQIKAEKRRGSLPDNSILHRLISELLPQIPERNSSLHALKRSPMHQPFHPLPPDGASHCPLYQNDCGRMPHSASFPDVDTTSNYHAQDYGSALSLQDHESPRSYSSTLTDLGRSASRERRGTPEKEKLPAKAVYDFKAQTSKELSFKKGDTVYILRKIDQNWYEGEHHGRVGIFPISYVEKLTPPEKAQPARPPPPVQPGEIGEAIAKYNFNADTNVELSLRKGDRIILLKRVDQNWYEGKIPGTNRQGIFPVSYVEVVKRNAKGAEDYPDPPLPHSYSSDRIYTLSSNKPQRPGFSHENIQGGGEPFQALYNYTPRNEDELELRESDVVDVMEKCDDGWFVGTSRRTKFFGTFPGNYVKRL"
alt_seqs = {
    "Q3UTJ2-2_0": {"seq_coord": (210, 211), "group": "alt_2"},
    "Q3UTJ2-2_1": {"seq_coord": (307, 308), "group": "alt_2"},
    "Q3UTJ2-3_0": {"seq_coord": (3, 34), "group": "alt_3"},
    "Q3UTJ2-3_1": {"seq_coord": (210, 211), "group": "alt_3"},
    "Q3UTJ2-3_2": {"seq_coord": (387, 914), "group": "alt_3"},
    "Q3UTJ2-3_3": {"seq_coord": (1125, 1180), "group": "alt_3"},
    "Q3UTJ2-4_0": {"seq_coord": (3, 34), "group": "alt_4"},
    "Q3UTJ2-4_1": {"seq_coord": (387, 914), "group": "alt_4"},
    "Q3UTJ2-5_0": {"seq_coord": (44, 67), "group": "alt_5"},
    "Q3UTJ2-5_1": {"seq_coord": (210, 211), "group": "alt_5"},
    "Q3UTJ2-5_2": {"seq_coord": (307, 308), "group": "alt_5"},
    "Q3UTJ2-5_3": {"seq_coord": (387, 914), "group": "alt_5"},
    "Q3UTJ2-5_4": {"seq_coord": (1125, 1135), "group": "alt_5"},
    "Q3UTJ2-6_0": {"seq_coord": (44, 67), "group": "alt_6"},
    "Q3UTJ2-6_1": {"seq_coord": (210, 211), "group": "alt_6"},
    "Q3UTJ2-6_2": {"seq_coord": (308, 316), "group": "alt_6"},
    "Q3UTJ2-6_3": {"seq_coord": (316, 1180), "group": "alt_6"},
    "Q3UTJ2-7_0": {"seq_coord": (310, 313), "group": "alt_7"},
    "Q3UTJ2-7_1": {"seq_coord": (313, 1180), "group": "alt_7"},
}

pr.pl.peptides_on_prot_sequence(
    adata,
    protein_id="Q3UTJ2",
    group_by="proteoform_id",
    alt_pep_sequence_key="peptide_sequence",
    ref_sequence=seq,
    add_sequences=alt_seqs,
    title="Sorbin and SH3 domain-containing protein 2 (Q3UTJ2)",
    figsize=(8, 4),
)
../_images/tutorials_bludau-2021_tissue-specific-proteoform-inference-across-five-mouse-organs_62_0.png 
[36]:

<Axes: title={'center': 'Sorbin and SH3 domain-containing protein 2 (Q3UTJ2)'}, xlabel='Position'>

Proteoform quantification and statistical analysis

To test if the inferred proteoform groups are indeed tissue-specific, a one-way ANOVA analysis can be performed directly using ProteoPy functions. In strong agreement with Bludau et al. (2021), our analysis revealed that 58 out of 63 proteoform-containing proteins (92.1%) are expressed in a tissue-specific manner.

[37]:

# Aggregate peptide intensities to proteoform level
adata_pfs = adata.copy()
pr.pp.quantify_proteoforms(adata_pfs, group_by="proteoform_id")

[38]:

# Retain only significant proteoforms (score >= 0.1, adj. p-value <= 0.1)
pf_mask = (adata_pfs.var["proteoform_score"].astype(float) >= 0.1) & (
    adata_pfs.var["proteoform_score_pval_adj"].astype(float) <= 0.1
)
adata_pfs = adata_pfs[:, pf_mask].copy()

[39]:

# Log2 transform for statistical testing
adata_pfs.layers["raw"] = adata_pfs.X
adata_pfs.X[adata_pfs.X == 0] = np.nan
adata_pfs.X = np.log2(adata_pfs.X)
adata_pfs.X[np.isnan(adata_pfs.X)] = 0

[40]:

# One-way ANOVA across tissues
pr.tl.differential_abundance(
    adata_pfs,
    method="anova_oneway",
    multitest_correction="bonferroni",
    group_by="tissue",
    alpha=0.01,
    space="log",
)

Saved test results in .varm['anova_oneway;tissue;all']

[41]:

anova_results = pr.get.differential_abundance_df(
    adata_pfs, keys="anova_oneway;tissue;all"
)
anova_results.rename(columns={"var_id": "proteoform_id"}, inplace=True)
anova_results

[41]:
proteoform_id test_type group_by design fstat pval mean_Brain mean_BAT mean_Heart mean_Liver mean_Quad pval_adj is_diff_abundant
0 O08601_0 anova_oneway tissue all 25.846223 5.018276e-10 14.604930 16.265425 12.130302 14.243520 11.643799 6.323028e-08 True
1 O08601_1 anova_oneway tissue all 280.923208 4.532385e-26 17.073096 16.195556 16.293561 19.305650 16.558551 5.710805e-24 True
2 O54724_0 anova_oneway tissue all 12.603239 1.878084e-06 15.879478 13.639894 8.489357 10.545068 11.828292 2.366386e-04 True
3 O54724_1 anova_oneway tissue all 179.779144 8.197678e-23 16.042224 19.217817 18.051066 15.596972 16.924093 1.032907e-20 True
4 O54749_0 anova_oneway tissue all 5.288221 1.938883e-03 13.162980 12.879554 14.428646 13.568651 13.071253 2.442993e-01 False
... ... ... ... ... ... ... ... ... ... ... ... ... ...
121 Q9QZ47_1 anova_oneway tissue all 91.893665 4.637443e-18 15.004601 15.819740 14.530142 14.795794 19.128501 5.843178e-16 True
122 Q9Z204_0 anova_oneway tissue all 75.275766 1.071660e-16 16.117296 12.403909 12.737527 13.374945 12.716591 1.350292e-14 True
123 Q9Z204_1 anova_oneway tissue all 152.715690 1.224770e-21 15.790335 14.611674 13.502690 16.977432 13.218752 1.543210e-19 True
124 Q9Z2Z6_0 anova_oneway tissue all 15.374666 2.380911e-07 13.995002 14.484431 13.383277 13.692943 12.821110 2.999948e-05 True
125 Q9Z2Z6_1 anova_oneway tissue all 112.326552 1.853570e-19 16.554941 19.275779 18.310678 18.047013 16.098757 2.335499e-17 True

126 rows × 13 columns

[42]:

# Count proteins with all proteoforms being significantly tissue-specific
protein_id_map = adata_pfs.var[["protein_id", "protein_id_old"]].set_index(
    "protein_id"
)["protein_id_old"]
anova_results["protein_id"] = anova_results["proteoform_id"].map(
    protein_id_map
)
n_tissue_specific_pfs = (
    anova_results.groupby("protein_id")["is_diff_abundant"].all().sum()
)

print(
    f"{n_tissue_specific_pfs} tissue-specific proteoform groups found via ANOVA."
)

58 tissue-specific proteoform groups found via ANOVA.

Summary

This notebook reproduced the proteoform inference workflow of Bludau et al. (2021) with ProteoPy:

  • Proteoform detection: COPF called 63 proteins as carrying multiple proteoform groups - the same count and the same 63 accessions as the paper, checked above: 0 missed, 0 false positives.

  • Biological verification: the tissue-specific assignments for Ldb3 (Q9JKS4) and Sorbs2 (Q3UTJ2) match the published findings (Figures 7A, 7C, 7D, 7F), with peptide-to-sequence mappings consistent with known alternative splice variants.

  • Tissue specificity: ANOVA recovered 58 tissue-specific proteoform groups (92.1 %), against the 56 (88.9 %) reported by Bludau et al.

Underneath all of this, ProteoPy’s COPF implementation is tested directly against the original R implementation’s intermediate outputs on this dataset — pairwise correlations, dendrograms, cluster assignments and proteoform scores — agreeing to between 1e-14 and 1e-12. That is what lets this notebook treat the peptide set as the only thing it needs to get right.

[43]:

!pip freeze

adjustText==1.3.0
alabaster==1.0.0
anndata==0.11.4
anyio==4.14.2
argon2-cffi==25.1.0
argon2-cffi-bindings==25.1.0
array-api-compat==1.15.0
arrow==1.4.0
asttokens==3.0.2
async-lru==2.3.0
attrs==26.1.0
babel==2.18.0
beautifulsoup4==4.15.0
biopython==1.86
bleach==6.4.0
certifi==2026.7.22
cffi==2.1.0
charset-normalizer==3.4.9
comm==0.2.3
contourpy==1.3.2
coverage==7.15.2
cycler==0.12.1
debugpy==1.8.21
decorator==5.3.1
defusedxml==0.7.1
docutils==0.21.2
et_xmlfile==2.0.0
exceptiongroup==1.3.1
executing==2.2.1
fastjsonschema==2.22.1
flake8==7.3.0
fonttools==4.63.0
fqdn==1.5.1
h11==0.16.0
h5py==3.16.0
httpcore==1.0.9
httpx==0.28.1
idna==3.18
igraph==1.0.0
imagesize==2.0.0
iniconfig==2.3.0
ipykernel==7.2.0
ipython==8.39.0
isoduration==20.11.0
jedi==0.20.0
Jinja2==3.1.6
joblib==1.5.3
json5==0.15.0
jsonpointer==3.1.1
jsonschema==4.26.0
jsonschema-specifications==2025.9.1
jupyter-events==0.12.1
jupyter-lsp==2.3.1
jupyter_client==8.9.1
jupyter_core==5.9.1
jupyter_server==2.20.0
jupyter_server_terminals==0.5.4
jupyterlab==4.5.6
jupyterlab_pygments==0.3.0
jupyterlab_server==2.28.0
kiwisolver==1.5.0
lark==1.3.1
latexcodec==3.0.1
legacy-api-wrap==1.5
llvmlite==0.46.0
markdown-it-py==3.0.0
MarkupSafe==3.0.3
matplotlib==3.10.8
matplotlib-inline==0.2.2
mccabe==0.7.0
mdit-py-plugins==0.6.1
mdurl==0.1.2
mistune==3.3.4
myst-parser==4.0.1
natsort==8.4.0
nbclient==0.11.0
nbconvert==7.17.1
nbformat==5.10.4
nbsphinx==0.9.8
nest-asyncio==1.6.0
networkx==3.4.2
notebook_shim==0.2.4
numba==0.64.0
numpy==2.2.6
numpydoc==1.10.0
openpyxl==3.1.5
overrides==7.7.0
packaging @ file:///home/conda/feedstock_root/build_artifacts/bld/rattler-build_packaging_1777103621/work
pandas==2.3.3
pandoc==2.4
pandocfilters==1.5.1
parso==0.8.7
patsy==1.0.2
pexpect==4.9.0
pillow==12.3.0
platformdirs==4.11.0
pluggy==1.6.0
plumbum==2.0.2
ply==3.11
pooch==1.9.0
prometheus_client==0.26.0
prompt_toolkit==3.0.53
-e git+https://github.com/UKHD-NP/proteopy.git@cbec4335d4607bcc0eb0a3b4b6341ca916325fb6#egg=proteopy
psutil==7.2.2
ptyprocess==0.7.0
pure_eval==0.2.3
pyarrow==23.0.1
pybtex==0.26.1
pybtex-docutils==1.0.3
pycodestyle==2.14.0
pycparser==3.0
pyflakes==3.4.0
Pygments==2.20.0
pynndescent==0.6.0
pyparsing==3.3.2
pytest==9.1.1
pytest-cov==7.1.0
python-dateutil==2.9.0.post0
python-json-logger==4.1.0
pytz==2026.3.post1
PyYAML==6.0.3
pyzmq==27.1.0
referencing==0.37.0
requests==2.34.2
rfc3339-validator==0.1.4
rfc3986-validator==0.1.1
rfc3987-syntax==1.1.0
rpds-py==0.30.0
scanpy==1.11.5
scikit-learn==1.7.2
scipy==1.15.3
seaborn==0.13.2
Send2Trash==2.1.0
session-info2==0.4.1
six==1.17.0
snowballstemmer==3.1.1
soupsieve==2.9.1
Sphinx==8.1.3
sphinx-autodoc-typehints==3.0.1
sphinx-copybutton==0.5.2
sphinx_design==0.6.1
sphinx_rtd_theme==3.1.0
sphinxcontrib-applehelp==2.0.0
sphinxcontrib-bibtex==2.7.0
sphinxcontrib-devhelp==2.0.0
sphinxcontrib-htmlhelp==2.1.0
sphinxcontrib-jquery==4.1
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==2.0.0
sphinxcontrib-serializinghtml==2.0.0
stack-data==0.6.3
statsmodels==0.14.6
terminado==0.18.1
texttable==1.7.0
threadpoolctl==3.6.0
tinycss2==1.5.1
tomli==2.4.1
tornado==6.5.7
tqdm==4.70.0
traitlets==5.16.0
typing_extensions==4.16.0
tzdata==2026.3
umap-learn==0.5.12
uri-template==1.3.0
urllib3==2.7.0
wcwidth==0.8.2
webcolors==25.10.0
webencodings==0.5.1
websocket-client==1.9.0