COMBSecretomics: A pragmatic methodological framework for higher-order drug combination analysis using secretomics

Abstract

Multi drug treatments are increasingly used in the clinic to combat complex and co-occurring diseases. However, most drug combination discovery efforts today are mainly focused on anticancer therapy and rarely examine the potential of using more than two drugs simultaneously. Moreover, there is currently no reported methodology for performing second- and higher-order drug combination analysis of secretomic patterns, meaning protein concentration profiles released by the cells. Here, we introduce COMBSecretomics (https://github.com/EffieChantzi/COMBSecretomics.git), the first pragmatic methodological framework designed to search exhaustively for second- and higher-order mixtures of candidate treatments that can modify, or even reverse malfunctioning secretomic patterns of human cells. This framework comes with two novel model-free combination analysis methods; a tailor-made generalization of the highest single agent principle and a data mining approach based on top-down hierarchical clustering. Quality control procedures to eliminate outliers and non-parametric statistics to quantify uncertainty in the results obtained are also included. COMBSecretomics is based on a standardized reproducible format and could be employed with any experimental platform that provides the required protein release data. Its practical use and functionality are demonstrated by means of a proof-of-principle pharmacological study related to cartilage degradation. COMBSecretomics is the first methodological framework reported to enable secretome-related second- and higher-order drug combination analysis. It could be used in drug discovery and development projects, clinical practice, as well as basic biological understanding of the largely unexplored changes in cell-cell communication that occurs due to disease and/or associated pharmacological treatment conditions.

Fig 1. COMBsecretomics conceptual workflow.

Disease associated (D) and healthy (H) cells are kept on the same experimental plate to avoid inter-plate variability. D cells are exposed to each and every treatment T from an exhaustive combination panel; here a panel of all 7 possible treatments using 3 pre-selected drugs T₁, T₂ and T₃ (at fixed concentrations) is shown as an example. D-treated and H cells are subsequently stimulated with each and every protein mixture S from a stimulation panel. Finally, release measurements for a protein panel of interest are collected for both cell types (f_{D,T, S} and f_H,S) using any technology that gives values proportional to the corresponding protein concentrations. The subsequent computational workflow include quality control procedures, normalization of the protein release differences (f_{D,T, S} − f_H,S), model-free higher-order combination analysis and non-parametric resampling statistics. The goal of all these methodological principles is to come up with an optimal (combination) treatment T* that reverses malfunctioning protein release patterns, meaning f_{D,T*, S} ≈ f_H,S.

Fig 2. Therapeutic need quantification.

Example showing the quantitative answers provided to question Q1b, using data from the OA case study presented later in this work. The external stimulations used, denoted S1; S2; S3, are shown in the y-axis, while the proteins measured are shown in the x-axis. The value in each patch is quantified by employing Eq (2) in the context of answering question Q1b.

Fig 3. Miniaturization of the resampling based leave-one-out validation approach per plate.

Left: original collected protein release d-dimensional vectors f for the six different cell states supported by COMBSecretomics ($H,T_o,S_o,H,T_o,S_y,D,T_o,S_o,D,T_o,S_y,D,T_x,S_o,D,T_x,S_y$). For each cell state, four intra-plate replicate measurements are shown as rows along with the corresponding replicate number. D cells are either untreated (T_o) or treated (T_x) while H cells are untreated (T_o). Both D and H cells are either stimulated (S_y) or not (S_o). Experimental wells of different cell states are colored differently. Right: N_v user-defined validation datasets are automatically created by employing a leave-one-out procedure among the four intra-plate replicate measurements for each cell state.

Fig 4. COMBSecretomics flowchart.

Experimental flow: disease associated (D) and healthy (H) cells are treated and stimulated in parallel on the same experimental plate p_i. An exhaustive combination panel is used for treating D cells, while a stimuli panel is employed for both D and H cells. An end point protein release assay of any sort can be used provided that it gives values proportional to the corresponding protein concentrations. Computational flow: a series of subsequent computational steps are employed for processing the protein release measurement values. Firstly, quality control procedures are employed. Secondly, protein release differences for stimulated and unstimulated cells are normalized per plate p_i. Then two model-free combination analysis methods are employed using the normalized protein release differences. Finally, non-parametric resampling statistics are used to quantify uncertainty for the obtained combination analysis results. Graphic and text files with all results are created and saved automatically.

Fig 5. GHSA analysis for stimulation S₁ and N_v = 10³ validation datasets.

All four combination treatments T₁₂, T₁₃, T₂₃, T₁₂₃ are shown on the x-axis. Each combination treatment is represented by a box plot showing the minimum value, 25^th, 50^th, 75^th percentiles and maximum value of all N_v GHSA indices obtained during a resampling-based leave-one-out validation approach.

Fig 6. Top-down hierarchical K-Means clustering for stimulation S₁ and N_v = 10³ validation datasets.

(a) Frequency/occurrence (%) for all unique partitions/clusters at the first hierarchical level. Annotations are provided for the three most frequent partitions. (b) Visualization of the clustering results after validation; the most dominant partitions at the first and second hierarchical levels are used. Each line corresponds to the centroid of each (sub-)cluster identified representing its restoration capacity defined by the normalized protein release differences between D, treated, stimulated with S₁ (D, T, S₁) and H, untreated, stimulated with S₁ (H, UT, S₁) cells. Solid and dotted lines are used for the first and second hierarchical level respectively. The gray solid line corresponds to the total therapeutic need, meaning the normalized protein release differences between D, untreated, stimulated with S₁ (D, UT, S₁) and H, untreated, stimulated with S₁ (H, UT, S₁) cells, which an ideal treatment should eradicate. The legend shows all treatments per (sub-)cluster identified without employing exhaustive subset search since only 7 treatments were used in total.