---
title: "Pattern Profiler: identifying shared protein profiles"
author: "Sergio Ciordia"
package: NADIA
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Pattern Profiler: identifying shared protein profiles}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  crop = NULL,
  message = FALSE,
  warning = FALSE
)

have_pattern_profiler <- all(vapply(
  c("Mfuzz", "Biobase", "e1071"),
  requireNamespace,
  logical(1),
  quietly = TRUE
))

# have_highcharter and prepare_vignette_highchart()
source("highcharts-setup.R")
have_pp_charts <- have_pattern_profiler && have_highcharter
```

# Introduction

Differential-expression analysis identifies proteins that change between pairs
of conditions. When an experiment contains three or more ordered conditions,
Pattern Profiler addresses a complementary question: **which proteins follow
the same expression pattern across the complete experiment?** Proteins that
increase, decrease, peak, or recover together may participate in related
biological processes even when their absolute abundances differ.

`pattern_profiler_analysis()` uses fuzzy c-means clustering, implemented with
the `Mfuzz` package (Kumar and Futschik, 2007). Unlike hard clustering, fuzzy
clustering does not force every protein into a single group. Instead, it assigns
a membership value between 0 and 1 to every protein--cluster pair. High
membership indicates a clear match to a pattern, whereas intermediate
membership identifies proteins with more ambiguous profiles.

This vignette shows how to select proteins, choose the number of clusters,
interpret memberships and expression profiles, and save the results for later
use.

# Setup and example data

```{r install, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("NADIA")
```

This vignette also uses `Mfuzz`, `highcharter`, `arrow` and `duckdb`, which
are optional;
`BiocManager::install("NADIA", dependencies = TRUE)` installs them together
with NADIA.

We use `nadia_dia`, the example DIA proteomics dataset distributed with NADIA.
The data are first processed to obtain both the intensity assays and the
differential-expression results required by Pattern Profiler.

```{r load-data}
library(NADIA)
data(nadia_dia)

res <- process_proteomics(nadia_dia, verbose = FALSE)
```

The processed experiment contains three ordered conditions (`A`, `B`, and `D`),
each represented by four samples. Replicates are retained during protein
quantification and are aggregated only when the condition-level profiles are
constructed.

<details class="code-details">
<summary>Show the code used to summarize the experimental design</summary>

```{r example-design, results='hide'}
design_table <- as.data.frame(SummarizedExperiment::colData(res$se_proc))
design_table <- as.data.frame(table(Condition = design_table$Condition))
names(design_table)[2] <- "Samples"

knitr::kable(
  design_table,
  caption = "Number of samples available for each condition."
)
```

</details>

```{r example-design-table, echo=FALSE}
knitr::kable(
  design_table,
  caption = "Number of samples available for each condition."
)
```

Pattern Profiler is most informative when the conditions have a meaningful
order, such as time points, doses, or disease stages. With only two conditions,
each profile contains a single transition and therefore provides little more
information than the direction of the corresponding fold change.

# Running Pattern Profiler

The function first selects proteins using the differential-expression table.
For each selected protein, it then aggregates replicates within conditions,
standardizes the resulting profile, selects or accepts a number of clusters,
and performs fuzzy c-means clustering. The example below evaluates two to five
clusters and uses the Xie--Beni index to choose among them.

<details class="function-parameters">
<summary>Arguments of <code>pattern_profiler_analysis()</code></summary>
<div class="parameter-content">
<dl>
<dt><code>se_proc</code></dt>
<dd>A processed <code>SummarizedExperiment</code> containing the protein-intensity assays and sample metadata.</dd>

<dt><code>DEPs_results</code></dt>
<dd>The differential-expression table returned by <code>process_proteomics()</code>.</dd>

<dt><code>assay_name</code></dt>
<dd>The assay used to build the profiles and to select the matching rows from <code>DEPs_results</code>.</dd>

<dt><code>filter_mode</code>, <code>alpha</code>, and <code>comparison</code></dt>
<dd>Control which proteins enter the analysis. The modes are explained in the next section; <code>comparison</code> is required only for <code>filter_mode = "specific"</code>.</dd>

<dt><code>condition_order</code></dt>
<dd>The order in which conditions appear in the profiles. If omitted, alphabetical order is used.</dd>

<dt><code>aggregate</code></dt>
<dd>How replicates are summarized within each condition: <code>"median"</code> (default) or <code>"mean"</code>.</dd>

<dt><code>auto_select_c</code>, <code>c_range</code>, and <code>c</code></dt>
<dd>Automatically evaluate candidate cluster numbers in <code>c_range</code>, or set a fixed <code>c</code> when <code>auto_select_c = FALSE</code>.</dd>

<dt><code>selection_method</code></dt>
<dd>The rule used for automatic selection: <code>"xb"</code>, <code>"consensus"</code>, or <code>"elbow"</code>.</dd>

<dt><code>min_membership</code></dt>
<dd>The minimum fuzzy membership required for a protein--cluster pair to be included in <code>long_output</code> (default: <code>0.25</code>).</dd>

<dt><code>n_starts</code></dt>
<dd>The number of random starts used to evaluate each candidate number of clusters; their metrics are averaged (default: <code>3</code>).</dd>

<dt><code>output_file</code></dt>
<dd>An optional path for writing <code>long_output</code> as a Parquet file. Nothing is written when it is <code>NULL</code>.</dd>

<dt><code>verbose</code></dt>
<dd>Controls progress messages.</dd>
</dl>
</div>
</details>

```{r run-pattern-profiler, eval=have_pattern_profiler, results='hide'}
set.seed(42)
pp <- pattern_profiler_analysis(
  se_proc          = res$se_proc,
  DEPs_results     = res$DEPs_results,
  assay_name       = "Impseqrob_min",
  filter_mode      = "any",
  alpha            = 0.05,
  condition_order  = c("A", "B", "D"),
  aggregate        = "median",
  c_range          = 2:5,
  auto_select_c    = TRUE,
  selection_method = "xb",
  min_membership   = 0.25,
  n_starts         = 3,
  verbose          = FALSE
)
```

The returned object contains the chosen cluster number, the estimated fuzzifier
(`m`), selection diagnostics, the fitted `Mfuzz` object, and a long-format table
for inspection and visualization. The overview below confirms that all selected
proteins passed the preprocessing checks. The output contains more rows than
proteins because a protein can meet the membership threshold in more than one
cluster.

<details class="code-details">
<summary>Show the code used to summarize the Pattern Profiler run</summary>

```{r run-overview, eval=have_pattern_profiler, results='hide'}
run_overview <- data.frame(
  Measure = c(
    "Proteins selected",
    "Proteins clustered",
    "Rows in long_output",
    "Selected clusters",
    "Fuzzifier (m)"
  ),
  Value = c(
    format(pp$n_features_input, scientific = FALSE),
    format(pp$n_features_final, scientific = FALSE),
    format(pp$n_rows_output, scientific = FALSE),
    format(pp$optimal_c, scientific = FALSE),
    format(round(pp$m, 3), nsmall = 3)
  )
)

knitr::kable(
  run_overview,
  caption = "Overview of the Pattern Profiler run."
)
```

</details>

```{r run-overview-table, eval=have_pattern_profiler, echo=FALSE}
knitr::kable(
  run_overview,
  caption = "Overview of the Pattern Profiler run."
)
```

# Selecting proteins for clustering

Clustering every quantified protein can obscure biologically responsive
patterns with a large number of stable profiles. `filter_mode` therefore lets
the user define the protein set before clustering.

| filter_mode | Proteins included |
|:---|:---|
| `"any"` | Proteins significant in at least one comparison at the selected `alpha`. |
| `"specific"` | Proteins significant in the comparison supplied through `comparison`. |
| `"all"` | All quantified proteins; no significance filter is applied. |

<div class="pattern-note">
<strong>Important:</strong> <code>filter_mode = "all"</code> does not mean
"significant in every comparison". It disables significance filtering and
includes every protein available in the selected assay.
</div>

For this dataset, the three modes produce substantially different input sizes.
The counts below are obtained before missing-value and variance checks are
applied by Pattern Profiler.

<details class="code-details">
<summary>Show the code used to compare the filtering strategies</summary>

```{r filter-counts, results='hide'}
dep_assay <- as.data.frame(res$DEPs_results)
dep_assay <- dep_assay[dep_assay$Assay == "Impseqrob_min", ]

filter_counts <- data.frame(
  Setting = c(
    'filter_mode = "any"',
    'filter_mode = "specific", comparison = "B-A"',
    'filter_mode = "all"'
  ),
  Proteins = c(
    length(unique(dep_assay$Protein.IDs[dep_assay$adj.P.Val <= 0.05])),
    length(unique(dep_assay$Protein.IDs[
      dep_assay$Comparison == "B-A" & dep_assay$adj.P.Val <= 0.05
    ])),
    nrow(res$se_proc)
  )
)

knitr::kable(
  filter_counts,
  caption = "Number of proteins selected by three filtering strategies."
)
```

</details>

```{r filter-counts-table, echo=FALSE}
knitr::kable(
  filter_counts,
  caption = "Number of proteins selected by three filtering strategies."
)
```

The `"any"` mode retains 794 proteins that respond in at least one comparison,
whereas a `"specific"` analysis of `B-A` would retain 497. The `"all"` mode
would include all 1,997 quantified proteins and is useful when coordinated
variation among non-significant proteins is also of interest, but it may make
the dominant clusters less specific to the biological response.

# Choosing the number of clusters

The number of clusters determines the resolution of the analysis. Too few
clusters can merge distinct profiles, whereas too many can divide one pattern
into small, weakly supported groups. When `auto_select_c = TRUE`, NADIA fits each
candidate value in `c_range` from several random starts and reports four complementary
diagnostics.

| Metric | What it evaluates | Preferred value |
|:---|:---|:---|
| XB | Within-cluster compactness relative to separation between centroids. | Lower |
| FPC | How clearly proteins are assigned rather than shared among clusters. | Higher |
| AMM | Average maximum membership across proteins. | Higher |
| Dmin | Minimum distance between any two cluster centroids. | Higher |

The default `selection_method = "xb"` minimizes the Xie--Beni (XB) index.
`"consensus"` combines the ranks of all four diagnostics and can be useful when
they disagree. `"elbow"` selects the largest change in centroid separation and
should be treated as an exploratory alternative rather than as a universal
optimum.

<details class="code-details">
<summary>Show the code used to prepare the cluster-selection diagnostics</summary>

```{r selection-metrics, eval=have_pattern_profiler, results='hide'}
selection_table <- pp$selection_metrics
selection_table[c("XB", "XB_sd", "FPC", "AMM", "Dmin")] <- lapply(
  selection_table[c("XB", "XB_sd", "FPC", "AMM", "Dmin")],
  round,
  digits = 4
)

knitr::kable(
  selection_table,
  row.names = FALSE,
  caption = "Cluster-selection diagnostics for c = 2 to 5."
)
```

</details>

```{r selection-metrics-table, eval=have_pattern_profiler, echo=FALSE}
knitr::kable(
  selection_table,
  row.names = FALSE,
  caption = "Cluster-selection diagnostics for c = 2 to 5."
)
```

All four diagnostics favor two clusters in this example: `c = 2` has the
lowest XB value and the highest FPC, AMM, and Dmin values. The selected solution
therefore provides two well-separated expression patterns without adding
poorly supported subdivisions. These metrics are guides rather than biological
proof; candidate solutions should still be checked visually and interpreted in
the context of the experiment.

<div class="pattern-note">
Fuzzy c-means starts from a random partition, drawn from R's random number
generator. Call <code>set.seed()</code> immediately before
<code>pattern_profiler_analysis()</code> to make a run reproducible, and raise
<code>n_starts</code> for a more stable choice of the number of clusters. Cluster
numbers are arbitrary labels and may be swapped between runs even when the
underlying patterns are equivalent.
</div>

Before clustering, NADIA removes profiles that cannot be standardized reliably,
including constant profiles and features with excessive missingness. Remaining
isolated missing values are filled by k-nearest-neighbor imputation. Failed
candidate fits are stored as missing diagnostic values rather than as zeros, so
they cannot be mistaken for exceptionally good solutions.

# Understanding the soft-clustering output

The main result is `long_output`. It contains the protein identifier, cluster,
membership, and one standardized value for each condition. Because expression
is converted to a row-wise z-score, these values describe the **relative profile
of each protein across the ordered conditions**; they must not be interpreted as
absolute protein abundance.

<details class="code-details">
<summary>Show the code used to preview the long-format output</summary>

```{r long-output, eval=have_pattern_profiler, results='hide'}
knitr::kable(
  head(pp$long_output, 6),
  digits = 3,
  caption = "First rows of the long-format Pattern Profiler output."
)
```

</details>

```{r long-output-table, eval=have_pattern_profiler, echo=FALSE}
knitr::kable(
  head(pp$long_output, 6),
  digits = 3,
  caption = "First rows of the long-format Pattern Profiler output."
)
```

Each row represents one protein--cluster association that satisfies
`min_membership`. `summarize_pattern_profiler()` provides an overview of these
associations and reports how many proteins occur in more than one cluster.

<details class="code-details">
<summary>Show the code used to summarize soft-clustering assignments</summary>

```{r membership-summary, eval=have_pattern_profiler, results='hide'}
pp_summary <- summarize_pattern_profiler(pp$long_output)

membership_overview <- data.frame(
  Measure = c(
    "Unique proteins",
    "Protein--cluster associations",
    "Proteins associated with multiple clusters"
  ),
  Value = c(
    pp_summary$n_unique_features,
    pp_summary$n_total_entries,
    pp_summary$n_multi_cluster_features
  )
)

knitr::kable(
  membership_overview,
  caption = "Summary of the soft-clustering assignments."
)
```

</details>

```{r membership-summary-table, eval=have_pattern_profiler, echo=FALSE}
knitr::kable(
  membership_overview,
  caption = "Summary of the soft-clustering assignments."
)
```

At the analysis threshold of 0.25, 73 proteins are associated with both
clusters. These overlapping assignments are not duplicates or errors: they
identify profiles that share characteristics with both centroids. The cluster
summary shows the size and strength of each set of assignments.

<div class="cluster-summary-table">

<details class="code-details">
<summary>Show the code used to summarize cluster membership</summary>

```{r cluster-summary, eval=have_pattern_profiler, results='hide'}
cluster_table <- as.data.frame(pp_summary$cluster_summary)
cluster_table[c("mean_membership", "min_membership", "max_membership")] <-
  lapply(
    cluster_table[c("mean_membership", "min_membership", "max_membership")],
    round,
    digits = 3
  )

knitr::kable(
  cluster_table,
  caption = "Membership summary for each cluster."
)
```

</details>

```{r cluster-summary-table, eval=have_pattern_profiler, echo=FALSE}
knitr::kable(
  cluster_table,
  caption = "Membership summary for each cluster."
)
```

</div>

Both clusters have high average memberships, indicating that the two-pattern
solution is well defined overall. A stricter display threshold can be used when
only the clearest representatives are needed.

<details class="code-details">
<summary>Show the code used to compare membership thresholds</summary>

```{r membership-thresholds, eval=have_pattern_profiler, results='hide'}
thresholds <- c(0.25, 0.50, 0.70, 0.90)
threshold_table <- data.frame(
  Minimum_membership = thresholds,
  Associations_retained = vapply(
    thresholds,
    function(x) sum(pp$long_output$Membership >= x),
    numeric(1)
  )
)

knitr::kable(
  threshold_table,
  caption = "Effect of applying stricter membership thresholds."
)
```

</details>

```{r membership-thresholds-table, eval=have_pattern_profiler, echo=FALSE}
knitr::kable(
  threshold_table,
  caption = "Effect of applying stricter membership thresholds."
)
```

Raising the threshold from 0.25 to 0.50 removes low-membership secondary
associations while retaining one best-supported assignment for every protein in
this two-cluster solution. At 0.90, only 554 highly characteristic associations
remain. The appropriate threshold depends on whether the aim is broad discovery
or a conservative list for downstream interpretation.

The threshold supplied to `pattern_profiler_analysis()` determines which rows
are retained in `long_output`. Plotting functions can apply an additional,
stricter threshold, but they cannot recover associations that were removed
during the analysis.

# Visualizing the expression patterns

Visualization connects the numerical memberships to their biological meaning.
NADIA provides interactive Highcharts plots for inspecting all protein profiles
within individual clusters and for comparing cluster centroids.

## Protein profiles within each cluster

`cluster_profile_highchart_list()` creates one interactive figure per cluster.
Individual protein lines reveal within-cluster heterogeneity, while the centroid
line, the mean or median of the proteins drawn (not the center fitted by the
clustering), summarizes the dominant pattern. Here, a membership threshold of 0.70 is used as
a conservative starting point for displaying reliable cluster profiles.

```{r individual-profile-code, eval = if (have_pp_charts) 1 else FALSE}
profile_plots <- cluster_profile_highchart_list(
  pp$long_output,
  conditions = pp$conditions,
  min_membership = 0.70,
  centroid_summary = "median"
)

profile_plots[[1]]
```

```{r cluster-shapes, include = FALSE}
# Size and shape of each cluster at the 0.70 threshold, so that the text below
# follows the clustering rather than describing one particular run of it: the
# numbering of the clusters is arbitrary and changes with the random starts.
if (have_pattern_profiler) {
  strong <- pp$long_output[pp$long_output$Membership >= 0.70, ]
  n_cluster <- tapply(strong$FeatureID, strong$Cluster,
                      function(x) length(unique(x)))
  centroid <- sapply(split(strong[, pp$conditions], strong$Cluster),
                     function(d) apply(d, 2, median))
  falls <- centroid["D", ] < (centroid["A", ] + centroid["B", ]) / 2
} else {
  n_cluster <- c("1" = NA, "2" = NA)
  falls <- c("1" = TRUE, "2" = FALSE)
}
shape <- ifelse(
  falls,
  "remain relatively stable between `A` and `B` and then decrease in `D`",
  "rise from `A` to `D`"
)
cluster_down <- names(falls)[falls][1]
cluster_up <- names(falls)[!falls][1]
```

The first cluster is shown below.

::: {.interactive-chart-compact .pattern-profiler-chart}
```{r cluster-profile-plot, echo = FALSE, eval = have_pp_charts}
prepare_vignette_highchart(profile_plots[[1]])
```
:::

The thin red lines represent the `r n_cluster[["1"]]` individual protein
profiles retained at this threshold, and the thick dark-red line is their median
centroid. Most proteins `r shape[["1"]]`. Their concentration around the
centroid indicates a coherent cluster pattern.

## Comparing cluster centroids

The centroid plot provides a compact overview of the patterns. Here the median
profile is used because it is less sensitive than the mean to proteins at the
edges of a cluster. The same membership threshold of 0.70 is applied so that
this summary and the cluster-level figure above represent the same protein
assignments.

```{r centroid-highchart-code, eval = if (have_pp_charts) 1 else FALSE}
centroid_plot <- cluster_centroids_highchart(
  pp$long_output,
  conditions = pp$conditions,
  min_membership = 0.70,
  centroid_summary = "median",
  title = "Pattern Profiler cluster centroids"
)

centroid_plot
```

::: {.interactive-chart-compact .pattern-profiler-chart}
```{r centroid-plot, echo = FALSE, eval = have_pp_charts}
prepare_vignette_highchart(centroid_plot)
```
:::

The centroids summarize `r n_cluster[["1"]]` proteins in cluster 1 and
`r n_cluster[["2"]]` in cluster 2. Cluster `r cluster_down` remains high in `A`
and `B` before decreasing in `D`, whereas cluster `r cluster_up` rises from `A`
to `D`. These opposing trajectories are expressed as within-protein
z-scores and therefore describe relative patterns rather than absolute
abundances.

<div class="figure-space"></div>

# Saving and reusing the results

Pattern Profiler results can be saved either as a standalone Parquet table or
inside a NADIA project. Both approaches avoid repeating the clustering when the
profiles need to be inspected or plotted again.

## Saving the long-format table as Parquet

Set `output_file` when running the analysis to write `long_output` directly.
Use the `.parquet` extension because the function writes a Parquet file.

```{r parquet-roundtrip, eval=have_pattern_profiler && requireNamespace("arrow", quietly=TRUE)}
parquet_file <- file.path(tempdir(), "pattern_profiler.parquet")

set.seed(42)
invisible(capture.output(pp_saved <- pattern_profiler_analysis(
  se_proc          = res$se_proc,
  DEPs_results     = res$DEPs_results,
  assay_name       = "Impseqrob_min",
  filter_mode      = "any",
  condition_order  = c("A", "B", "D"),
  auto_select_c    = FALSE,
  c                = pp$optimal_c,
  min_membership   = 0.25,
  output_file      = parquet_file,
  verbose          = FALSE
)))

reloaded <- read_pattern_profiler_data(
  parquet_file,
  min_membership = 0.70
)

data.frame(
  Rows_saved = nrow(pp_saved$long_output),
  Rows_reloaded_at_0.70 = nrow(reloaded)
)

unlink(parquet_file)
```

`read_pattern_profiler_data()` validates the required columns and can apply a
stricter membership threshold while loading the file. In this example, the
saved table contains all associations at 0.25, whereas the reloaded object keeps
only those at or above 0.70. This stricter filter retains 748 of the 867 saved
protein--cluster associations.

## Storing Pattern Profiler in a NADIA project

A `.nadia` project can store the preprocessing, differential-expression
analysis, and Pattern Profiler result together. This is convenient when the
complete analysis must be transferred or reopened as one reproducible object.

```{r nadia-roundtrip, eval=have_pattern_profiler && requireNamespace("duckdb", quietly=TRUE)}
project_file <- file.path(tempdir(), "pattern_profiler_example.nadia")

write_nadia(
  project_file,
  preprocessing = nadia_dia,
  result = res,
  verbose = FALSE
)

nadia_add_pattern_profiler(
  project_file,
  pattern_profiler = pp,
  verbose = FALSE
)

stored_project <- read_nadia(project_file)
stored_profiles <- nadia_pattern_profiler(stored_project)

data.frame(
  Rows_original = nrow(pp$long_output),
  Rows_recovered = nrow(stored_profiles)
)

unlink(c(project_file, paste0(project_file, ".wal")))
```

The matching row counts confirm that the complete long-format clustering table
was recovered. Alternatively, `pattern_profiler = pp` can be supplied directly
to `write_nadia()` when the project and clustering are saved at the same time.

To combine Pattern Profiler assignments with differential-abundance statistics
in a single table, for example as input to a downstream functional or enrichment
analysis, use `deps_with_clusters()` on the two results in memory or
`nadia_deps_with_clusters()` on a stored project. Both functions are described
in `vignette("results-and-export", package = "NADIA")`.

# Practical recommendations

Pattern Profiler is most useful as an exploratory layer after a well-controlled
differential-expression analysis. The following practices make the results
easier to reproduce and interpret:

- provide an explicit `condition_order` that reflects the experimental design;
- begin with `filter_mode = "any"`, then broaden or restrict the protein set to
  answer a specific biological question;
- inspect all cluster-selection diagnostics rather than treating one index as
  definitive;
- call `set.seed()` before each run, and check that the main profiles remain
  stable across seeds and plausible cluster numbers;
- use membership to distinguish representative proteins from ambiguous ones;
- interpret standardized trajectories as relative shapes, not abundance
  differences; and
- perform pathway or functional enrichment only after verifying that the
  corresponding cluster is coherent and biologically interpretable.

With three or more ordered conditions, this workflow complements pairwise
differential expression by revealing groups of proteins that share a response
across the experiment. It does not establish regulation or causality by itself,
but it provides focused protein sets for subsequent biological validation.

# References

Kumar L, Futschik ME (2007). “Mfuzz: a software package for soft clustering of
microarray data.” *Bioinformation*, **2**(1), 5--7.

# Session information

The package versions used to build this vignette are recorded below to support
reproducibility.

```{r session-information}
sessionInfo()
```