Introduction

The ClusterSignificance package provides tools to assess if “class clusters” in dimensionality reduced data representations have a separation different from permuted data. The term “class clusters” here refers to, clusters of points representing known classes in the data. This is particularly useful to determine if a subset of the variables, e.g. genes in a specific pathway, alone can separate samples into these established classes. Evaluation of this is accomplished in a three stage process Projection, Separation classification, and Score permutation. To be able to compare class cluster separations, we give them a score based on this separation. First, all data points are projected onto a line (Projection), after which the best separation for two groups at a time is identified and scored (Separation classification). Finally, to get a p-value for the separation, we compare the separation score for our real data to the separation score for permuted data (Score permutation).

The package includes 2 different methods for accomplishing the Projection step outlined above, Mean line projection (Mlp) and Principal curve projection (Pcp). Here we will first underline the assumptions made by the ClusterSignificance method, followed by a description of the similarities of the three stages of ClusterSignificance independent of which Projection method is used. This is followed by an example of the Mlp and Pcp methods and a description of the unique features of each. Furthermore, we provide an example where ClusterSignificance is used to examine the separation between 2 class clusters downsteam of a principal component analysis (PCA).

Assumptions

The ClusterSignificance package operates based on three simple assumptions shown below:

1. The dimensionality reduction method was sufficiently capable of detecting the dissimilarities, characterised by features in high-dimensional space, that correspond to the class separations under investigation.

2. The obtained principal curve accurately depicts the data with respect to the class separations under consideration.

3. Under the null hypothesis, the joint distribution remains invariant under all rearrangements of the subscripts.

Methods

Choosing a projection method

As previously mentioned, there are two projection methods provided with the ClusterSignificance package; Principal curve projection (Pcp) and Mean line projection (Mlp) . These are outlined below together with the situations where each can be used.

Examples

library(ClusterSignificance)

data(mlpMatrix)

data(pcpMatrix)

## Create the group variable.
classes <- rownames(pcpMatrix)

## Run Pcp and plot.
prj <- pcp(pcpMatrix, classes)
plot(prj)

cl <- classify(prj)
plot(cl)

getData(cl, "AUC")

## class1 vs. class3 class2 vs. class3 class1 vs. class2 
##           0.61250           0.95750           0.96875

#permute matrix
iterations <- 100 ## Set the number of iterations.

## Permute and plot.
pe <- permute(
    mat = pcpMatrix, 
    iter = iterations, 
    classes = classes, 
    projmethod = "pcp"
)

plot(pe)

pvalue(pe)

## class1 vs class2 class1 vs class3 class2 vs class3 
##             0.01             0.11             0.01

conf.int(pe, conf.level = .99)

##            class1 vs class2 class1 vs class3 class2 vs class3
## CI99%-low        0.00000000       0.04450171       0.00000000
## CI99%-high       0.05160403       0.21454670       0.05160403

## Create the group variable.
classes <- rownames(mlpMatrix)

## Run Mlp and plot.
prj <- mlp(mlpMatrix, classes)
plot(prj)

## Classify and plot.
cl <- classify(prj)
plot(cl)

## Set the seed and number of iterations.
iterations <- 500 

## Permute and plot.
pe <- permute(
    mat = mlpMatrix, 
    iter = iterations, 
    classes = classes, 
    projmethod = "mlp"
)

plot(pe)

## class1 vs class2 
##            0.002

##            class1 vs class2
## CI99%-low        0.00000000
## CI99%-high       0.01054069

Principal Component Analysis

One probable application for ClusterSignificance would be the analysis of cluster separation post-PCA. Here we use the leukemia dataset from the plsgenomics package to demonstrate how this may be done. The dataset includes “gene expression data (3051 genes and 38 leukemia samples) from the leukemia microarray study of Golub et al. (1999)”. The tumor samples are from both acute lymphoblastic leukemia (ALL) and acute myeloid leukemia (AML) patients. Due to the inherent differences in the transcriptome of these two cancer types, we hypothesise that they would form separate clusters post-PCA and that we can use ClusterSignificance to examine if this separation is significant.

## The leukemia dataset is provided as a list with 3 elements. 
## An explaination of these can be found using the command:
## help(leukemia)

library(plsgenomics, quietly = TRUE)
data(leukemia)

## Extract the expression matrix.
expression <- leukemia$X 

## Run PCA and extract the principal components.
pca <- prcomp(expression, scale = TRUE)
mat <- pca$x

## Extract the grouping variable (coded as 1 (ALL) and 2 (AML)).
classes <- ifelse(leukemia$Y == 1, "ALL", "AML")

ob <- pcp(mat, classes, normalize = TRUE)
plot(ob)

cl <- classify(ob)
plot(cl)

iterations <- 100
pe <- permute(mat = mat, iter = iterations, classes = classes, projmethod = "pcp")
pvalue(pe)

## ALL vs AML 
##       0.01

conf.int(pe)

##            ALL vs AML
## CI99%-low  0.00000000
## CI99%-high 0.05160403

As we can see from the results, at alpha 0.05, we would reject the null hypothesis concluding that there is a separation between the AML and ALL class clusters. Despite this, by examining the confidence interval, we see that the upper limit is greater than 0.05 which should give us pause before rejecting the null at alpha = 0.05. Although not run here, 10000 iterations results in p = < 1e-04 with the confidence interval of CI99%-low = 0.00000 and CI99%-high = 0.00052 giving us greater confidence to reject the null.

Common questions

Do I have to use dimensionality reduced data as input to ClusterSignificance.

Well, technically no, although ClusterSignificance was designed to function downstream of dimensionality reduction algorithms and we see this being the most frequent use-case for ClusterSignificance. One of the reasons for this is, it is often the case that datasets containing a large number of features, such as those originating from high throughput experiments, will often contain features that cannot discriminate between the class separations of interest. This is often accompanied by “background noise” arising during the sample preparation Dimensionality reduction thus serves here as a feature selection method, identifying sources of large variation within the data. Despite this, there are some potential pitfalls with dimensionality reduction which are highlighted in the next question.

What does the assumption “The obtained principal curve accurately depicts the data with respect to the class separations under consideration” mean?

We believe the principal curve method to be a robust method that will perform well in the majority of situations we predict ClusterSignificance to be used in. Despite this, the principal curve method is an unsupervised method and therefore, it is not necessarily the case that it will be drawn in an optimal way to access the class separations of interest. We currently have some functionality allowing the user to regulate the smoothness of the principal curve. This can be accomplished using the df argument to the pcp function BUT you must supply the same argument to the permute function or your results will not be valid. We recommend altering the df argument as a last resort and to only be used under situations where no other alternative is available.

What does the assumption “The dimensionality reduction method was sufficiently capable of detecting the dissimilarities, characterised by features in high-dimensional space, that correspond to the class separations under investigation” mean?

Dimensionality reduction can be said to detect dissimilarities, characterised by features in high-dimensional space although, there is no guarantee that the dissimilarities it detects will be those discriminating the class separations you are interested in. In addition, linear dimensionality reduction methods, such as PCA, may not provide a useful representation of the data when the discriminatory features between classes under investigation vary in a non-linear fashion.

What does the assumption “The joint distribution remains invariant under all rearrangements of the subscripts.” mean?

This is a standard assumption for permutation tests which means that the data must exhibit “exchangeability”. A basic description of exchangeability can be found here. Extensions Of The Concept Of Exchangeability And Their Applications Good (2012) provides a detailed discussion on this topic.

Does ClusterSignificance actually perform a permutation test?

Despite the name of one of ClusterSignificance’s stages being permutation, strictly speaking, ClusterSignificance is not designed to perform a permutation test. Due to the fact that a real permutation test tests all possible combinations without replacement, which is often infeasible, most users will actually be performing a monte carlo test. Instead of permuting all possible combinations to form the null distribution (or permutation distribution), a monte carlo test samples the permutation distribution and uses this sample as the null distribution. For this reason, a p-value from a monte carlo test is actually a estimate of the true p-value. For this reason we provide the conf.int function that will calculate the confidence interval for the p-value estimate. It is theoretically, possible to perform a permutation test with ClusterSignificance by providing your own permutation matrix instead of letting ClusterSignificance calculate one for you which we discuss in the Advanced usage section.

How many permutations should I run?

The complete answer to this question is too in-depth to answer here. We would recommend 10^4 permutations as a starting point. Ideally, the permutation test returns a exact p-value (i.e. not p-value < X), which means that at least one permuted max score was greater than the max score for the real data, although this is not always possible. Furthermore, it is unnecessary to run more permutations than possible class combinations, which provides one concrete way to calculate the number of ideal permutations. Despite this, it will often challenge your computers memory to calculate this number if many combinations are possible. Under such circumstances, it is common that the permutation distribution is randomly sampled and the more times this sampling occurs, the nearer your estimate of the p-value will come to the true p-value (see the question above). A confidence interval can be calculated for the reported p-value using the conf.int function.

Why does the principal curve change shape between steps 3 and 4 of the projection phase?

The line is actually the same shape, it is the scale of the X-axis that changes. This is necessary due to the fact that as the points move to one dimension more space is needed on the x-axis to show all the points.

Why are the points projected diagonally onto the line in step 4 of the Mlp projection?

Step 4 shows the orthogonal projection onto the regression line and appears visually to not be perpendicular to the line. This is due to the automatic scaling of the y-axis which is necessary when large differences in proportions between the x and y dimensions would otherwise cause data to disappear outside of the margins.

Advanced usage

mat <- mlpMatrix
classes <- rownames(mlpMatrix)
iterations <- 10

pe1 <- permute(mat=mat, iter=iterations, classes=classes, projmethod="pcp")
pe2 <- permute(mat=mat, iter=iterations, classes=classes, projmethod="pcp")
pe3 <- permute(mat=mat, iter=iterations, classes=classes, projmethod="pcp")
pvalue(pe1)

## class1 vs class2 
##              0.1

pe <- c(pe1, pe2, pe3)
pvalue(pe)

## class1 vs class2 
##       0.03333333

##define permMatrix function
.pca <- function(y, classes, uniq.groups, mat) {
    pca <- prcomp(
        sapply(1:ncol(mat[classes %in% uniq.groups[, y], ]), 
            function(i)
                mat[classes %in% uniq.groups[, y], i] <- 
                    sample(
                        mat[classes %in% uniq.groups[, y], i], 
                        replace=FALSE
                    )
        ), scale=TRUE
    )
    return(pca$x)
}

permMatrix <- function(iterations, classes, mat) {
    uniq.groups <- combn(unique(classes), 2)
    permats <- lapply(1:ncol(uniq.groups), 
        function(y)
            lapply(1:iterations, 
                function(x)
                    .pca(y, classes, uniq.groups, mat)
            )
    )
    return(permats)
}

set.seed(3)
mat <- pcpMatrix
classes <- rownames(pcpMatrix)
iterations = 100
permats <- permMatrix(iterations, classes, mat)

pe <- permute(
    mat = mat, 
    iter=100, 
    classes=classes, 
    projmethod="pcp",
    userPermats=permats
)

Conclusion

The ClusterSignificance package can estimate a p-value for a separation projected on a line. Typically, the package might be used post-PCA-analysis, but clusters coming from other sources can also be tested.

Session info

sessionInfo()

## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] plsgenomics_1.5-3          ClusterSignificance_1.35.0
## [3] ggplot2_3.5.1              BiocStyle_2.35.0          
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.9           utf8_1.2.4           stringi_1.8.4       
##  [4] pracma_2.4.4         digest_0.6.37        magrittr_2.0.3      
##  [7] evaluate_1.0.1       grid_4.4.1           RColorBrewer_1.1-3  
## [10] fastmap_1.2.0        maps_3.4.2           plyr_1.8.9          
## [13] jsonlite_1.8.9       princurve_2.1.6      BiocManager_1.30.25 
## [16] fansi_1.0.6          spam_2.11-0          viridisLite_0.4.2   
## [19] scales_1.3.0         RhpcBLASctl_0.23-42  jquerylib_0.1.4     
## [22] cli_3.6.3            rlang_1.1.4          scatterplot3d_0.3-44
## [25] munsell_0.5.1        withr_3.0.2          cachem_1.1.0        
## [28] yaml_2.3.10          tools_4.4.1          parallel_4.4.1      
## [31] reshape2_1.4.4       colorspace_2.1-1     boot_1.3-31         
## [34] buildtools_1.0.0     vctrs_0.6.5          R6_2.5.1            
## [37] lifecycle_1.0.4      stringr_1.5.1        MASS_7.3-61         
## [40] pkgconfig_2.0.3      pillar_1.9.0         bslib_0.8.0         
## [43] gtable_0.3.6         glue_1.8.0           Rcpp_1.0.13         
## [46] fields_16.3          xfun_0.48            tibble_3.2.1        
## [49] highr_0.11           sys_3.4.3            knitr_1.48          
## [52] htmltools_0.5.8.1    rmarkdown_2.28       maketools_1.3.1     
## [55] dotCall64_1.2        compiler_4.4.1