Structural Bioinformatics Library
Template C++ / Python API for developping structural bioinformatics applications.
User Manual


Authors: A. Lheritier and F. Cazals


Comparing two sets of multivariate samples is a central problem in data analysis. From a statistical standpoint, one way to perform such a comparison is to resort to a non-parametric two-sample test (TST), which checks whether the two sets can be seen as i.i.d. samples of an identical unknown distribution (the null hypothesis, denoted H0).

Accepting or rejecting H0 requires comparing the value obtained for the test statistic used against a pre-defined threshold, usually denoted $ \alpha $. In rejecting H0, further analysis are called for, and in particular:

  • Effect size calculation: one wishes to assess the magnitude of the difference.

  • Feedback analysis: one wishes to identify regions in parameter space accounting for this difference.

Effect size and feedback analysis are easy for one-dimensional data, as one may use the effect size associated with a Mann-Whitney U-test, or may perform a density estimation and compute some kind of distance (say the Mallows distance) between the obtained distributions. However, such analysis are more complex for multivariate data.

This package fills this gap, proposing a method based upon a combination of statistical learning (regression) and computational topology methods [24] . More formally, consider two populations, each given as a point cloud in $ \Rnt^d $. Our analysis proceeds in three steps:

  • Step 1: a discrepancy measure is estimated at each data point.

  • Step 2: a clustering of samples is performed. A cluster reported involves samples with significant discrepancy (threshold user-defined), which in addition are in spatial proximity.

  • Step 3: the effect size is decomposed in a sum of contributions per cluster.

(A) (B)
(C) (D)
(E) (F)
Method: input and main output. (A) Point clouds–here uniformly sampling two ellipsis. (B) Local Discrepancy on a scale 0-1 computed at each sample point. (C) Persistence Diagram associated with the sublevel sets of the estimated discrepancy (see text for details). (D) Simplified Disconnectivity Forest illustrating the merge events between sublevel sets. Note that four clusters stand out (the four legs of the tree). (E) Clustering: clusters identified on the disconnectivity forest, mapped back into the ambient space. Note that the green and blue clusters are associated with the red population, while the remaining two come from the blue population. (F) Effect Size Decomposition. The area under the dashed line represents the total effect size. The area under the continuous line is the maximum effect size possible. The area of each bar represents the contribution of each cluster to the effect size. (G) Discrepancy Cumulative Distribution.

Using the programs


We aim at modeling the discrepancy between two datasets $ \seqLS{x}{\nbsamplesX} \equiv {x_1,\dots,x_{\nbsamplesX} } $ and $ \seqLS{y}{\nbsamplesY} \equiv {y_1,\dots,y_{\nbsamplesY} } $, in some fixed dimension Euclidean space $ \Rnt^d $. We denote $ n=\nbsamplesX+\nbsamplesY $ the total number of samples. We view these data as coming from two unknown densities $ \densd{X} $ and $ \densd{Y} $, et denote $ \dens $ the average density $ \dens \equiv \nicefrac{(\densdX+\densdY)}{2} $.

For the sake of conciseness, the X and Y populations are respectively assigned the labels 0 and 1, and are respectively represented in blue and red.

Step 1

In the first step, we assign a label to each set and we compute, for each sample point, a discrepancy measure based on comparing an estimate of the conditional probability distribution of the label given a position versus the global unconditional label distribution.

Divergence measures between distributions. To define our discrepancy measure, recall that the Kullback-Leibler divergence (KLD) between two probability densities $ f $ and $ g $ defined as

$ \divKL{f}{g} \equiv \int_{-\infty}^\infty f(x) \log \frac{f(x)}{g(x)} dx $

with the conventions $ 0 \log 0 = 0 $ and $ 0 \log \frac{0}{0}=0 $.

For two discrete distributions $ P $ and $ Q $ over a finite set $ \alphabet $, the the Kullback-Leibler divergence is defined by:

$ \divKL{P}{Q} \equiv \sum_{\rvl\in\alphabet} P(\rvl) \log \frac{P(\rvl)}{Q(\rvl)} $

The Jensen-Shannon divergence (JSD) [79] allows one to symmetrize and smooth the KLD by taking the average KLD of $ \densdX $ and $ \densdY $ to the average density $ \dens $, that is:

$ \JSD{\densdX}{\densdY} \equiv \frac{1}{2}\left(\divKL{\densdX}{\dens} + \divKL{\densdY}{\dens}\right) $

Discrepancy measure. In taking the average in Eq. (eq-JSD), two random variables are implicity defined: a position variable $ \rvP $ with density $ \densdZ \equiv \dens $ and a binary label $ \rvL $ that indicates from which original density (i.e. $ \densdX $ or $ \densdY $) an instance of $ \rvP $ is obtained. Formally, considering the alphabet $ \mathcal{A} = {\labelX,\labelY} $ and $ X\sim F_X,Y\sim F_Y $, the following pair of random variables is defined:

$ (\rvL,\rvP)= \begin{cases} (\labelX,X) & \text{ with prob. } \frac{1}{2} \\ (\labelY,Y) & \text{ with prob. } \frac{1}{2} \end{cases} $

In the sequel, we will consider the conditional and unconditional mass functions $ \pmflz=\probaCond{\rvL=\rvl}{\rvP=\rvp} $ and $ \pmfl=\proba{\rvL=\rvl}=\frac{1}{2} $ respectively, as well as the joint probability density $ \densdLZ $. We will also use the notation $ \densd{\rvl} $ to denote $ \densdX $ (resp.~ $ \densdY $) when $ \rvl=\labelX $ (resp.~ $ \rvl=\labelY $).

The connexion between the JS divergence and the previous (conditional) mass functions is given by [24] :

One has:
$ \JSD{\densdX}{\densdY} = \int_{\R^d} \densdZ{z} \divKL{\cpmf[\cdot][\rvp]}{\cpmf[\cdot]} dz $

The previous lemma shows that the JSD can be seen as the average, over $ z\in\R^d $, of the KLD between the conditional and unconditional distribution of labels. More formally, we define:

The discrepancy at location $ z $ is defined as the KL divergence:
$ \deltaz \equiv \divKL{\cpmf[\cdot][\rvp]}{\cpmf[\cdot]}. $

The following comments are in order:

  • $ \deltaz $ ranges between 0 and 1 and is 0 iff $ \densdX{\rvp}=\densdY{\rvp} $.

  • Since $ \pmfl $ is known but $ \pmflz $ is not, the problem we consider now is the one of estimating $ \pmflz $ at each given location $ \rvp $. In the sequel, such an estimator is denoted $ \cestimi{n}{\rvl}{\rvp} $.

Taking for granted such an estimator (see section Discrepancy estimation for its effective calculation), we finally obtain an estimator for $ \deltaz $ :

$ \deltazestim{n} \equiv \divKL{\cestimi{n}{\cdot}{\rvp}}{\cpmf[\cdot]} $

In practice, this estimate is computed using a regressor based on nearest neighbors, see section Discrepancy estimation . The number of neighbors used for each samples is set to $ n^{\kofnpower} $, with $ \kofnpower $ a user defined parameter.

Joint distribution compatible sampling. Practically, the sizes $ n_0 $ and $ n_1 $ of the samples processed may be unbalanced, which would not be compatible with Eq. (eq-pair-LZ).

To ensure that balanced data sets are used in the estimation, we resort to random multiplexer generating samples $ (L,Z) $, using a Bernoulli distribution to pick from the two data sets (Fig. fig-random-gadget ; see also details in [24] ).

The multiplexer halts when the data on one of the two data sets has been exhausted, and the resulting pairs are used to perform the estimation. Note that selected samples from the non exhausted data set may not be used, yielding some information loss. Therefore, we may repeat the process $ B $ times, so as to eventually use all samples from the larger set. Upon performing these $ B $ estimates, the discrepancy estimate for a given sample is taken as the median of the $ B $ estimates.

For a given estimation fold, we denote $ n' \leq n_0+n_1 $ the number of samples selected by the random multiplexer. The regressor used to estimate the discrepancy uses these $ n' $ samples only. On the other hand, using this model, the discrepancy is estimated on all $ n=n_0+n_1 $ samples. See also Fig. fig-workflow-ddbc for the articulation between the three steps.

Step 2

To describe the second step in simple terms–see [24] for the details, we refer to the so-called watershed transform used in image processing or topographical analysis, see Watershed image processing and Drainage basin. In a nutshell, given a height function, called the landscape for short in the sequel, the watershed transform partitions its into its catchment basins, with one basin for each local minimum.

To transpose these notions in our setting, we proceed in two steps. First, we build a nearest neighbor graph connecting the samples—e.g. we connect each sample to its $ k_n $ nearest neighbors. Second, we lift this graph into $ \mathbb{R}^{d+1} $ by adding a new coordinate to each sample, namely its estimated discrepancy. This lifted graph provides a discrete representation for the aforementioned landscape. We now describe our method using this landscape (Fig. fig-samples-to-height-function ).

On the landscape, we focus on samples whose estimated discrepancy is less than $ -\deltaMax $, with $ \deltaMax $ a user defined threshold. The region of the landscape matching this specification is called the sublevel set associated with $ -\deltaMax $. We define a cluster as a connected component of this sublevel set. Note that each such cluster may be associated with the lowest local minimum of this connected component. However, in case of noisy estimates, many small cluster may be reported. To ease the interpretation of results, we therefore smooth the landscape, and focus on c.c. of the smoothed landscape. Smoothing is carried out using topological persistence [54] , which consists of canceling local minima associated with non significant basins. More precisely, recall that the persistence of a local minimum is the height difference to the lowest saddle connecting this local minimum to a deeper local minimum. We retain local minima whose persistence is above a threshold $ \rho $.

To recap, our analysis aims at:

  • Identifying samples with significant discrepancy.

  • Grouping such samples into clusters, the samples in one cluster being associated with significant local maxima of the estimated discrepancy.

Method: algorithm illustrated on a toy 1D example
(A, discrepancy estimation) A discrepancy is estimated $ \hat{\delta}(z_i) $ at each sample $ z_i $ (see text for details). (B, clusters as connected components) Upon reverting the sign of the discrepancy, clusters are defined as connected components (c.c.) of the region of the previous curve below a user defined threshold $ -\deltaMax $. On this example, four such components are observed, two large ones and two small ones. Note that each cluster can be charged to one local minimum of the function $ -\hat{\delta} $. (C, smoothed clusters upon applying topological persistence) To get rid of small clusters, the curve $ -\hat{\delta} $ is smoothed using topological persistence. This simplification yields two clusters. Note that the blue sample squeezed between $ C_1 $ and $ C_2 $ does not appear in the red cluster since its discrepancy is below $ \deltaMax $ (and likewise for the red point with respect to $ C_3 $ and $ C_4 $).

Step 3

This step consists of adding up the discrepancies on a per sample basis in a cluster $ C $:

$ \JSDi{C}{\densdX}{\densdY} \equiv \frac{1}{n}\sum_{z \in \seqLS{z}{n} \cap C}\deltazestim{n} $
$ \JSDi{C}{\densdX}{\densdY} \equiv \frac{1}{\nbsamplesX+\nbsamplesY}\sum_{z \in (\seqLS{x}{\nbsamplesX} \cup \seqLS{y}{\nbsamplesY}) \cap C} \hat{\delta}(\rvp). $


Putting everything together, the program $ \text{\feedback} $ performs the following plots:

  • Persistence diagram: A plot showing for each minima a red cross with coordinates $ (x,y) $ corresponding to its birth and death dates respectively, while analyzing the landscape whose elevation is the negative estimated discrepancy.

  • Discrepancy cumulative distribution function: The cumulative distribution function of the estimated discrepancies.

  • Disconnectivity forest: a hierarchical representation of the merge events associated with sublevel sets of the estimated discrepancy.

  • Clusters: A plot similar to the raw data plot, with one color per cluster. The points not belonging to any cluster are colored in gray.

  • JSD decomposition plot: A 1D plot presenting a synthetic view without relying on the MDS embedding. The $ x $ coordinate represents the sample space and always ranges from 0 to
    1. The total estimated JSD is represented by the area under the dashed line. And the maximum possible JSD which is always 1 is represented by the area under the continuous line. One bar is depicted for each cluster plus another one (the last one) for the points not belonging to any cluster. The area of each bar represents the contribution of the corresponding cluster to the total JSD and its color corresponds to the proportion of samples $ \labelX $ in the cluster.

For high dimensional date, we also provide the following 2D embeddings:

  • Raw data embedding: For samples embedded in 2D or 3D space, a plot of the points with a color to indicate the label (blue: $ \labelX $; red: $ \labelY $). For samples embedded in a higher dimensional space, a 2D embedding of these samples obtained using multi-dimensional scaling (MDS). In any case, the goal of this plot is to intuitively visualize the distributions of the two populations.

  • Discrepancy shaded data embedding — aka discrepancy plot: A plot similar to the raw data plot, except that each sample is color coded using a heat palette from fully transparent white to red across yellow, as a function of the value of the estimated discrepancy $ \deltazestim{n} $. That is, a point with $ \deltazestim{n} $ equal to zero (resp.~one) is colored fully transparent white (resp.~non transparent red).

Random multiplexer generating pairs (label, position).

Input: Specifications and File Types

In the sequel, we specify the main steps and the options offered, see also Fig. fig-workflow-ddbc .

Step 1: using and

The first step is carried out by the python program $ \text{\feedback} $. Its main input are two files listing points. A basic calculation is launched as follows:

> -f data/ellipse_w_2_h_10_n_1000__X.dat  -g data/ellipse_w_2_h_10_n_1000__Y.dat  --theta0 0.5  --kofn-power  0.66 -B 1 --directory results
The main options of the program $ \text{\feedback} $ are:
-f string: First input file in matrix format (one row = one observation/sample). NB: these samples are assigned the label 0.
-g string: Second input file in matrix format, as above. NB: these samples are assigned the label 1.
–directory string(= "."): Directory for the output files
–theta0 float(= 0.5): theta0 parameter for the random multiplexer
–kofn-power float(= 0.66): Power p in kofn=k^p, to set the number of neighbors used by the regressor
-B float(= 1): Number of iterations of the resampling process

File Name


Points text file A list of points in a given dimension.
Points text file A list of points in a given dimension.
Input files for the first step of the run described in section Input: Specifications and File Types .

Optional step. This step aims to provide plots that may be used in the third step. It is carried out by the python program $ \text{\ddbcPlot} $. Its main input is the prefix of the output files of $ \text{\feedback} $. A basic calculation is launched as follows:

> -f results/ellipse_w_2_h_10_n_1000__X_ellipse_w_2_h_10_n_1000__Y
The main options of the program $ \text{\ddbcPlot} $ are:
-f string: File prefix

See Table table-ddbc-output-step-1 to see the input files of $ \text{\ddbcPlot} $ .

Step 2: using

The second step is carried out using the python program $ \text{\ddbcClustering} $. It uses the executable $ \text{\sblMTBAnngeuclid} $, from the SBL, and its main input is the set of points with their discrepancy, output of $ \text{\feedback} $. A basic calculation is launched as follows:

> -p results/ellipse_w_2_h_10_n_1000__X_ellipse_w_2_h_10_n_1000__Y.pointsWdim -w results/ellipse_w_2_h_10_n_1000__X_ellipse_w_2_h_10_n_1000__Y.heightDiv -n 6 -P --rho 0.1 --deltamax 0.1 -d results -o ellipse_w_2_h_10_n_1000__X_ellipse_w_2_h_10_n_1000__Y -v
The main options of the program $ \text{\ddbcClustering} $ are:
-p string: Height function: samples i.e. observations, in point d format (.pointsWdim file format, see step 1)
-w string: Height function: elevation i.e. discrepancy computed at step 1 (.heightDiv file format, see step 1)
-d string(= "."): Directory for output files (default is working directory)
-o string: Prefix for output files
-n int: Number of neighbors used to build the Nearest Neighbors Graph
-P: Order preserving perturbation of heights to break ties between neighboring samples
–rho float: 'Persistence threshold rho (default: no persistence)
–deltamax float: Sublevelset threshold deltaMax (default: no sublevelset)
-v: Verbose mode

Note that if the executable $ \text{\sblMTBAnngeuclid} $ is not in the path, the path to $ \text{\sblMTBAnngeuclid} $ may be specified using the options -e. See Table table-ddbc-output-step-1 to see the input files of $ \text{\ddbcClustering} $ .

Step 3: using

The third step is carried out using the python program $ \text{\ddbcPlotClustering} $. Its main input is the clustering made by $ \text{\ddbcClustering} $, with the output points of the program $ \text{\feedback} $ . A basic calculation is launched as follows:

> -f results/ellipse_w_2_h_10_n_1000__X_ellipse_w_2_h_10_n_1000__Y.vertices -c results/ellipse_w_2_h_10_n_1000__X_ellipse_w_2_h_10_n_1000__Y__sublevelsets_connected_components.xml -p results/ellipse_w_2_h_10_n_1000__X_ellipse_w_2_h_10_n_1000__Y.2Dpoints
The main options of the program $ \text{\ddbcPlotClustering} $ are:
-f string: Vertices filename generated in step 1 (.vertices file format, see step 1)
-c string: Connected components filename generated in step 2 (*__sublevelsets_connected_components.xml file format, see step 2)
-p string: Optional 2D embedding generated by (.2Dpoints file format)

See Tables table-ddbc-output-step-1 , table-ddbc-output-optional-step and table-ddbc-output-step-2 to see the input files of $ \text{\ddbcPlotClustering} $ .

Output: Specifications and File Types

Step 1.

PreviewFile Name Description
Points text file A list of discrepancy, one per point.
Weights text file A list of discrepancy, one per point.
Text file Classification of points...
Output files for the first step of the run described in section Input: Specifications and File Types .

Optional step.

PreviewFile Name


Points text file A list of 2D points from the input points.
Click itPopulations image 2D representation of the populations.
Click itDiscrepancy image 2D representation of the discrepancy of the points.
Output files for the optional step of the run described in section Input: Specifications and File Types .

Step 2.

PreviewFile Name


Disconnectivity Forest as EPS file Image of the disconnectivity forest.
Persistence Diagram as PDF file Image of the persistence diagram.
Clusters as XML file Clusters of points.
NNG as XML file Nearest Neighbors Graph.
Output files for the second step of the run described in section Input: Specifications and File Types .

Note that there are other files printed by the program $ \text{\ddbcClustering} $ that are not of interest for this application. We refer the reader to the package Morse_theory_based_analyzer for more details about these files.

Step 3.

PreviewFile Name Description
Click itClustering image 2D representation of the clusters.
Click itBarplots image Barplot representation of the clustering.
Click itDiscrepancy Cumulative Distribution image Curve representing the discrepancy cumulative distribution function.
Output files for the third step of the run described in section Input: Specifications and File Types .


Mixture of Gaussians

An example illustrating the processing of a mixture of Gaussians is provided on Fig. fig-gaussians-in-out .

In this example – see [24] for the detailed specification, the first population is drawn according to a mixture of four spherical Gaussians with equal probability, while the second one is drawn according to a mixture of four non-spherical Gaussians with equal probability

The different steps are run as following: -f data/gaussianMixture6__d_2_n_2000_X.dat -g data/gaussianMixture6__d_2_n_2000_Y.dat --theta0 0.5 --kofn-power 0.66 -B 1 --directory results -f results/gaussianMixture6__d_2_n_2000_X_gaussianMixture6__d_2_n_2000_Y -p results/gaussianMixture6__d_2_n_2000_X_gaussianMixture6__d_2_n_2000_Y.pointsWdim -w results/gaussianMixture6__d_2_n_2000_X_gaussianMixture6__d_2_n_2000_Y.heightDiv -n 6 -P --rho 0.1 --deltamax 0.25 -d results -o gaussianMixture6__d_2_n_2000_X_gaussianMixture6__d_2_n_2000_Y -v -f results/gaussianMixture6__d_2_n_2000_X_gaussianMixture6__d_2_n_2000_Y.vertices -c results/gaussianMixture6__d_2_n_2000_X_gaussianMixture6__d_2_n_2000_Y__sublevelsets_connected_components.xml -p results/gaussianMixture6__d_2_n_2000_X_gaussianMixture6__d_2_n_2000_Y.2Dpoints

(A) (B)
(C) (D)
(E) (F)
Processing a mixture of Gaussians. The panels are identical to those in Fig. fig-ellipsis-in-out . The following values were used: persistence simplification at threshold $ \rho=0.1 $; clustering with $ \deltaMax=0.25 $.

Higher dimensional case

An example illustrating the processing of real data in high dimension, we process a mixture of handwritten digits, each represented as an $ 28\times 28 $ gray-scale image – see the MNIST dataset [75] .

The two populations see also fig-digits-in-out , were chosen as follows:

  • blue population: 100 digits 3, 500 digits 6, and 1000 digits 8.
  • red population: 1000 digits 3, 500 digits 6, ans 100 digits 8.

The different steps are run as following: -f data/mnist368_100_500_1000_X.dat -g data/mnist368_1000_500_100_Y.dat --theta0 0.5 --kofn-power 0.66 -B 1 --directory results -f results/mnist368_100_500_1000_X_mnist368_1000_500_100_Y -p results/mnist368_100_500_1000_X_mnist368_1000_500_100_Y.pointsWdim -w results/mnist368_100_500_1000_X_mnist368_1000_500_100_Y.heightDiv -n 30 -P --rho 0.1 --deltamax 0.34 -d results -o mnist368_100_500_1000_X_mnist368_1000_500_100_Y -v -f results/mnist368_100_500_1000_X_mnist368_1000_500_100_Y.vertices -c results/mnist368_100_500_1000_X_mnist368_1000_500_100_Y__sublevelsets_connected_components.xml -p results/mnist368_100_500_1000_X_mnist368_1000_500_100_Y.2Dpoints

(A) (B)
(C) (D)
(E) (F)
Processing a mixture of handwritten digits represented as images (digits: 3, 6, 8). The panels are identical to those in Fig. fig-ellipsis-in-out . The following values were used: persistence simplification at threshold $ \rho=0.1 $; clustering with $ \deltaMax=0.34 $.

Algorithms and Methods

Discrepancy estimation

In this section, we explain how to compute the estimator of Eq. (eq-deltaz), using the framework of regression [63] , and more precisely regressors based on $ k_n $ nearest neighbors. We define:

Given a random vector $ (Z, R) $, where $ Z\in\Rd $ and the response variable $ R\in\R $, the $ \text{\emph{regression function}} $ is defined as
$ m(z)=\condexpecd{R}{Z=z}. $

In the regression problem, the goal is to build an estimator $ m_n(z) $ of $ m(z) $ using a set of $ n $ i.i.d.~realizations of $ (Z,R) $.

In order to apply this framework to our problem, note that the correspondence $ R=L $ yields

$ m(\rvp)=\cpmf[1][\rvp]. $

Then, we can use the following estimator for $ \pmflz $

$ \cestimi{n}{\rvl}{\rvp} \equiv | 1 - \rvl - m_n(\rvp) |. $

Note that it is required that $ 0\leq m_n(\rvp) \leq 1 $, since we aim at estimating a conditional probability, and that it is satisfied by the $ k_n $-nearest neighbor regressor.

Using Eq. (eq-estim-pmflz), we finally obtain an estimator for $ \deltaz $ :

$ \deltazestim{n} \equiv \divKL{\cestimi{n}{\cdot}{\rvp}}{\cpmf[\cdot]} $

Persistence analysis

Precisely defining persistence is beyond the scope of this user manual, and the reader is referred to [54] , as well as to [24] for all details.

We comment on three issues, though:

  • Practically, our algorithm operates in a discrete setting wince the input consists of samples. The samples are used to build a nearest neighbor graph (NNG), by connecting each sample to a number of nearest neighbors (for a suitable distance metric), from which local minima and saddles can be identified [33] .

  • Equipped with a NNG, reporting the connected components (c.c.) of a sublevel set is a mere application of the so-called union-find algorithm, whose complexity is essentially linear in the number of samples processed. See [33] for details.

  • For large datasets though, in exploratory data analysis, the previous union-find algorithm may be called for varying values $ -\deltaMax $. In that case, a more suitable strategy consists of working on the critical points of the height function (the estimated discrepancy), once all samples which are regular points have been assigned to their so-called stable manifolds. The persistence can then be run on the graph of critical points only (the Morse-Smale-Witten complex), as explained in [22] .

Programmer's Workflow

Method: workflow.
In blue: the parameters.

External dependencies

The following two external dependencies are used in this applications: