One essential functionality in the analysis of large scale datasets is the ability to perform dimensionality reduction techniques and plot the outputs, such as Principal Component Analysis (PCA), t-distributed Stochastic Neighbor Embedding (t-SNE) or Uniform Manifold Approximation and Projection (UMAPs). The mesa package currently supports functionality for PCAs and UMAPs.
Basic use
The function getPCA calculates a PCA on a qseaSet:
pca0 <- getPCA(exampleTumourNormal)
#> ------------------------------
#> Initial number of windows = 819.
#> -----------
#> Generating table of beta values for 819 regions across 10 samples.
#> -----------
#> Filtered out 220 windows with at least one missing value: 599 windows remaining.
#> ------------------------------
#> The following topVarNum values are larger than the number of remaining windows (= 599): 1000
#> These values are not used; PCA will be done with all remaining windows instead.
#> ------------------------------
#> No filtering of windows based on window standard deviation.
#> Performing PCA with 10 samples and 599 windows
#> -> pca1.
#> ------------------------------
By default, this will perform the PCA using beta values, and consider the 1000 most variable windows across all samples (or all windows if there are less than 1000).
We then can use the plotPCA() function to plot the result.
plotPCA(pca0)
#> $pca1
#> $pca1$PC1vsPC2
#>
#> $pca1$PC2vsPC3
The output of getPCA() is the PCA result object (technically a mesaDimRed object). This includes a copy of the sampleTable, detached from the original qseaSet.
As shown above, plotPCA() returns two separate plots by default; one for PC1 vs PC2 and another for PC2 vs PC3.
PlotPCA function arguments
Anotating samples by colour and/or shape
The colour and shape arguments to plotPCA() specify variables to use to set the colour and/or shape of the points (samples) on the PCA plot. These variables need to be in the qseaSet sample table.
Multiple variables can be specified for colour and a PCA plot will be produced for each variable.
Only one variable can be specified for shape.
plotPCA(pca1,
colour = c("patient", "age", "transformed_age"),
shape = "tumour")
#> $pca1
#> $pca1$patient
#> $pca1$patient$PC1vsPC2
#>
#> $pca1$patient$PC2vsPC3
#>
#>
#> $pca1$age
#> $pca1$age$PC1vsPC2
#>
#> $pca1$age$PC2vsPC3
#>
#>
#> $pca1$transformed_age
#> $pca1$transformed_age$PC1vsPC2
#>
#> $pca1$transformed_age$PC2vsPC3
Here, many PCA plots have been produced and displayed (PC1 vs PC2 and PC2 vs PC2 for three different colour variable annotations. To have more control over which plots are displayed you can save the plots to a variable and then choose which ones to display.
PCAplots1 <- plotPCA(pca1,
colour = c("patient", "age", "transformed_age"),
shape = "tumour")
PCAplots1$pca1$patient$PC1vsPC2
PCAplots1$pca1$age$PC1vsPC2
PCAplots1$pca1$transformed_age$PC1vsPC2
Colour and shape palettes
The function aims to select sensible colour and shape palettes based on the properties of the variable being used.
For the transformed_age variable (that was generated in the mutate example above), a diverging colour palette is used, but it is not symmetric around zero. That is, the darkest blue maps to approx -0.5 whereas the darkest red maps to approx 2.
Depending on the variable, it may be appropriate that the darkest blue maps to -2 (i.e. equivalent to the darkest red mapping to 2), and for this the symDivColourScale argument can be set to TRUE, as seen below (and we also replot the above version for comparison).
PCAplots2 <- plotPCA(pca1, colour = "transformed_age", shape = "tumour", symDivColourScale = TRUE)
PCAplots1$pca1$transformed_age$PC1vsPC2
PCAplots2$pca1$transformed_age$PC1vsPC2
By default, filled shapes (21-25) with a black outline are used for points unless a categorical variable with more than 5 categories (including NA) is set to the shape, in which case the unfilled shapes are used (0-20)*. Filled shapes are generally preferred, as they remain more distinguishable when many points overlap.
*Note: Shapes 19 and 20 are excluded from default palettes as they are difficult to distinguish from shape 16.
As a rule, we do not allow mixing of filled and unfilled shapes on the same plot. This is because the variable provided to the colour argument is mapped to the fill aesthetic for filled shapes and the colour aesthetic for unfilled shapes in the underlying ggplot2 geom.
plotPCA(pca1,
colour = "tumour",
shape = "tissue")
#> $pca1
#> $pca1$tumour
#> $pca1$tumour$PC1vsPC2
#>
#> $pca1$tumour$PC2vsPC3
plotPCA(pca1,
colour = "tumour",
shape = "stage")
#> $pca1
#> $pca1$tumour
#> $pca1$tumour$PC1vsPC2
#>
#> $pca1$tumour$PC2vsPC3
It is possible to specify your own colour and shape palettes using the colourPalette and shapePalette arguments.
plotPCA(pca1,
colour = "tumour",
colourPalette = c("firebrick","blue"),
shape = "tumour",
shapePalette = c(15, 8))
#> $pca1
#> $pca1$tumour
#> $pca1$tumour$PC1vsPC2
#>
#> $pca1$tumour$PC2vsPC3
These can be a named vector, in order to specify which value of the annotation gets which colour/shape:
plotPCA(pca1,
colour = "tumour",
colourPalette = c("Tumour" = "firebrick", "Normal" = "blue"),
shape = "tumour",
shapePalette = c("Tumour" = 15, "Normal" = 8))
#> $pca1
#> $pca1$tumour
#> $pca1$tumour$PC1vsPC2
#>
#> $pca1$tumour$PC2vsPC3
NA values
If a variable used for colour annotation contains NA values then these are displayed with a dark grey colour by default (this can be changed using the NAcolour argument).
If a variable used for shape annotation contains NA values then these are displayed explicitly with shape 25 by default for a filled shape palette and shape 7 by default for an unfilled shape palette (this can be changed using the NAshape argument).
plotPCA(pca1, components = c(1,2), colour = "stage")
#> $pca1
#> $pca1$stage
#> $pca1$stage$PC1vsPC2
plotPCA(pca1, components = c(1,2), shape = "stage", NAshape = 5)
#> $pca1
#> $pca1$PC1vsPC2
Components
As can be seen in the last example above, the components argument to plotPCA() can be used to specify which components to plot. Above we only ask for PC1 vs PC2.
If you would like more than one plot, then use a list. For example, the default behaviour is to set components = list(c(1, 2), c(2, 3)) which plots both PC1 vs PC2 and PC2 vs PC3.
Point size
This is controlled using the pointSize argument, which has a default value of 2.
plotPCA(pca1,
components = c(1, 2),
colour = "patient",
shape = "tumour")
#> $pca1
#> $pca1$patient
#> $pca1$patient$PC1vsPC2
plotPCA(pca1,
components = c(1, 2),
colour = "patient",
shape = "tumour",
pointSize = 4)
#> $pca1
#> $pca1$patient
#> $pca1$patient$PC1vsPC2
Text annotations
Sample names can be displayed for each point by setting showSampleNames = TRUE.
plotPCA(pca1,
components = c(1, 2),
colour = "patient",
shape = "tumour",
showSampleNames = TRUE)
#> $pca1
#> $pca1$patient
#> $pca1$patient$PC1vsPC2
Alternatively, plotly::ggplotly() can be used to generate an interactive version of the plot which gives information for each point when you rollover the point with the mouse. The plotlyAnnotations argument can be used to specify which variables to display in the pop-up boxes. The sample name will always be displayed, as will any variables being used colour and shape annotations.
PCAplots3 <- plotPCA(pca1,
components = c(1, 2),
colour = "patient",
shape = "tumour",
plotlyAnnotations = c("stage", "age", "gender"))
Plotly plots are great for interactive use in html reports.
The downside of using them is that the legend often gets displayed in a complicated, unhelpful way and the formatting of the plot is often not as good and hard to control.
Sometimes it makes sense to plot a non-plotly and plotly version of the same plot.
getPCA() function arguments
By default, the getPCA function will only consider a subset of the most variable windows across the qseaSet. This has two benefits, the first being computational efficiency, but the second is to prevent the effects of minor batch effects from appearing when they are summed over many windows (although that is also worth investigating separately). This is seen particularly when using beta values, as how close a beta value can get to 0 or 1 is determined by the enrichment and number of fragments sequenced; a sample with 5 million fragments might be assigned a beta value of 0.95 in a window, but 0.99 in the same window if it had 10 million fragments.
Setting number of top variable windows to use
This selection of the number of windows can be set using the topVarNum argument of getPCA(). If not explicitly set, the default value is the top 1000 most variable windows.
More than one value can be specified for topVarNum and PCA results will be produced for each. This is computationally more efficient than running getPCA() multiple times with a single topVarNum value each time.
Specifying a value of NA will result in all windows being used (i.e. no filtering based on most variable windows).
pca2 <- getPCA(exampleTumourNormal, topVarNum = c(100, 250, 500, NA))
#> ------------------------------
#> Initial number of windows = 819.
#> -----------
#> Generating table of beta values for 819 regions across 10 samples.
#> -----------
#> Filtered out 220 windows with at least one missing value: 599 windows remaining.
#> ------------------------------
#> Calculating standard deviation for each window across all 10 samples:
#> Colon1_N, Colon1_T, Colon2_N, Colon2_T, Lung1_N, Lung1_T, Lung2_N, Lung2_T, Lung3_N, Lung3_T.
#> -> column name windowSd1.
#> ------------------------------
#> Filtering windows based on standard deviation across all 10 samples (windowSd1).
#> Standard deviation threshold = 0.19 resulting in 100 windows.
#> Performing PCA with 10 samples and 100 windows
#> -> pca1.
#> ------------------------------
#> Filtering windows based on standard deviation across all 10 samples (windowSd1).
#> Standard deviation threshold = 0.13 resulting in 250 windows.
#> Performing PCA with 10 samples and 250 windows
#> -> pca2.
#> ------------------------------
#> Filtering windows based on standard deviation across all 10 samples (windowSd1).
#> Standard deviation threshold = 0.0485 resulting in 500 windows.
#> Performing PCA with 10 samples and 500 windows
#> -> pca3.
#> ------------------------------
#> No filtering of windows based on window standard deviation.
#> Performing PCA with 10 samples and 599 windows
#> -> pca4.
#> ------------------------------
Here, four different values have been used and the messages output tell you the name assigned to each PCA result.
When plugging this output into plotPCA, plots are produced for all four PCA versions, with all options for the colour parameter:
PCAplots4 <- plotPCA(pca2, colour = c("patient", "age"), shape = "tumour")
PCAplots4$pca1$patient$PC1vsPC2
PCAplots4$pca2$age$PC1vsPC2
PCAplots4$pca3$patient$PC2vsPC3
PCAplots4$pca4$age$PC1vsPC2
If you specify a value larger than the number of available windows, then PCA will be performed with all available windows and redundant values will be ignored:
pca <- getPCA(exampleTumourNormal, topVarNum = c(250, 500, 750, 1000, NA))
#> ------------------------------
#> Initial number of windows = 819.
#> -----------
#> Generating table of beta values for 819 regions across 10 samples.
#> -----------
#> Filtered out 220 windows with at least one missing value: 599 windows remaining.
#> ------------------------------
#> The following topVarNum values are larger than the number of remaining windows (= 599): 750, 1000
#> These values are not used; PCA is already being done with all remaining windows.
#> Calculating standard deviation for each window across all 10 samples:
#> Colon1_N, Colon1_T, Colon2_N, Colon2_T, Lung1_N, Lung1_T, Lung2_N, Lung2_T, Lung3_N, Lung3_T.
#> -> column name windowSd1.
#> ------------------------------
#> Filtering windows based on standard deviation across all 10 samples (windowSd1).
#> Standard deviation threshold = 0.13 resulting in 250 windows.
#> Performing PCA with 10 samples and 250 windows
#> -> pca1.
#> ------------------------------
#> Filtering windows based on standard deviation across all 10 samples (windowSd1).
#> Standard deviation threshold = 0.0485 resulting in 500 windows.
#> Performing PCA with 10 samples and 500 windows
#> -> pca2.
#> ------------------------------
#> No filtering of windows based on window standard deviation.
#> Performing PCA with 10 samples and 599 windows
#> -> pca3.
#> ------------------------------
Samples used to calculate variability
It is possible to select which samples are used to calculate the standard deviation for each window using the topVarSamples argument. This may be either an explicit vector of sample names,
tumourNames <- exampleTumourNormal |> filter(tumour == "Tumour") |> getSampleNames()
getPCA(exampleTumourNormal, topVarNum = c(50), topVarSamples = tumourNames)
#> ------------------------------
#> Initial number of windows = 819.
#> -----------
#> Generating table of beta values for 819 regions across 10 samples.
#> -----------
#> Filtered out 220 windows with at least one missing value: 599 windows remaining.
#> ------------------------------
#> Calculating standard deviation for each window across 5 of 10 samples:
#> Colon1_T, Colon2_T, Lung1_T, Lung2_T, Lung3_T.
#> -> column name windowSd1.
#> ------------------------------
#> Filtering windows based on standard deviation across 5 samples (windowSd1).
#> Standard deviation threshold = 0.254 resulting in 50 windows.
#> Performing PCA with 10 samples and 50 windows
#> -> pca1.
#> ------------------------------
#> Object containing 1 dimensionality reduction objects for 10 samples
Or a string to apply (as a regular expression):
getPCA(exampleTumourNormal, topVarNum = c(50), topVarSamples = "_T")
#> ------------------------------
#> Initial number of windows = 819.
#> -----------
#> Generating table of beta values for 819 regions across 10 samples.
#> -----------
#> Filtered out 220 windows with at least one missing value: 599 windows remaining.
#> ------------------------------
#> Calculating standard deviation for each window across 5 of 10 samples:
#> Colon1_T, Colon2_T, Lung1_T, Lung2_T, Lung3_T.
#> -> column name windowSd1.
#> ------------------------------
#> Filtering windows based on standard deviation across 5 samples (windowSd1).
#> Standard deviation threshold = 0.254 resulting in 50 windows.
#> Performing PCA with 10 samples and 50 windows
#> -> pca1.
#> ------------------------------
#> Object containing 1 dimensionality reduction objects for 10 samples
The default is to use all samples.
Inputting a pre-calculated dataTable
The slowest part of the PCA is often generating the table of values (e.g. beta values), which we often refer to as the dataTable. For large numbers of windows and/or samples this can take a considerable amount of time and could become a bottleneck if you need to run getPCA() multiple times.
The getPCA() function therefore allows for the dataTable to be generated externally and then input into getPCA(). This means getPCA() can be run multiple times with the dataTable generated only once.
pca1 <- getPCA(exampleTumourNormal) # dataTable generated within getPCA
#> ------------------------------
#> Initial number of windows = 819.
#> -----------
#> Generating table of beta values for 819 regions across 10 samples.
#> -----------
#> Filtered out 220 windows with at least one missing value: 599 windows remaining.
#> ------------------------------
#> The following topVarNum values are larger than the number of remaining windows (= 599): 1000
#> These values are not used; PCA will be done with all remaining windows instead.
#> ------------------------------
#> No filtering of windows based on window standard deviation.
#> Performing PCA with 10 samples and 599 windows
#> -> pca1.
#> ------------------------------
For example, we can use getBetaTable to generate the table to use:
dataTable <- getBetaTable(exampleTumourNormal)
#> Generating table of beta values for 819 regions across 10 samples.
pca2 <- getPCA(exampleTumourNormal, dataTable = dataTable) # dataTable generated externally and passed into getPCA
#> Generating table of beta values for 50 regions across 5 samples.
#> ------------------------------
#> Initial number of windows = 819.
#> Filtered out 220 windows with at least one missing value: 599 windows remaining.
#> ------------------------------
#> The following topVarNum values are larger than the number of remaining windows (= 599): 1000
#> These values are not used; PCA will be done with all remaining windows instead.
#> ------------------------------
#> No filtering of windows based on window standard deviation.
#> Performing PCA with 10 samples and 599 windows
#> -> pca1.
#> ------------------------------
And this gives the same result:
identical(pca1, pca2)
#> [1] TRUE
Another option is to run getPCA() and set the returnDataTable argument to TRUE. The dataTable will then be part of the function output and can be passed into any further calls to getPCA().
The dataTable is not given as part of the output by default since the dataTable can often be large and storing it every time getPCA() is run would be memory-inefficient.
pca3 <- getPCA(exampleTumourNormal, returnDataTable = TRUE) # dataTable generated within getPCA and returned as part of output
#> ------------------------------
#> Initial number of windows = 819.
#> -----------
#> Generating table of beta values for 819 regions across 10 samples.
#> -----------
#> Filtered out 220 windows with at least one missing value: 599 windows remaining.
#> ------------------------------
#> The following topVarNum values are larger than the number of remaining windows (= 599): 1000
#> These values are not used; PCA will be done with all remaining windows instead.
#> ------------------------------
#> No filtering of windows based on window standard deviation.
#> Performing PCA with 10 samples and 599 windows
#> -> pca1.
#> ------------------------------
methods::slot(pca1,"dataTable") # NULL
#> data frame with 0 columns and 0 rows
methods::slot(pca3, "dataTable") # dataTable has been stored in output
#> # A tibble: 599 × 14
#> seqnames start end CpG_density Colon1_N Colon1_T Colon2_N Colon2_T
#> <fct> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 7 25005001 25005300 13.4 0.0565 0.0171 0.0275 0.0165
#> 2 7 25017601 25017900 8.38 0.122 0.0695 0.918 0.818
#> 3 7 25017901 25018200 17.9 0.0240 0.0660 0.265 0.250
#> 4 7 25018201 25018500 16.5 0.0472 0.0143 0.112 0.0669
#> 5 7 25018501 25018800 17.5 0.0241 0.0140 0.0426 0.0133
#> 6 7 25018801 25019100 15.3 0.0493 0.0149 0.0692 0.0559
#> 7 7 25019101 25019400 18.2 0.0932 0.0380 0.220 0.0646
#> 8 7 25019401 25019700 17.1 0.0466 0.0528 0.0881 0.0258
#> 9 7 25035001 25035300 9.29 0.0855 0.327 0.308 0.389
#> 10 7 25047301 25047600 9.54 0.225 0.0787 0.0722 0.0447
#> # ℹ 589 more rows
#> # ℹ 6 more variables: Lung1_N <dbl>, Lung1_T <dbl>, Lung2_N <dbl>,
#> # Lung2_T <dbl>, Lung3_N <dbl>, Lung3_T <dbl>
We can then use the dataTable calculated in a previous getPCA call, but ask getPCA to return a different number of top windows:
pca4 <- getPCA(exampleTumourNormal,
dataTable = methods::slot(pca3, "dataTable"),
topVarNum = 250)
#> Generating table of beta values for 50 regions across 5 samples.
#> ------------------------------
#> Initial number of windows = 599.
#> No windows have missing values.
#> ------------------------------
#> Calculating standard deviation for each window across all 10 samples:
#> Colon1_N, Colon1_T, Colon2_N, Colon2_T, Lung1_N, Lung1_T, Lung2_N, Lung2_T, Lung3_N, Lung3_T.
#> -> column name windowSd1.
#> ------------------------------
#> Filtering windows based on standard deviation across all 10 samples (windowSd1).
#> Standard deviation threshold = 0.13 resulting in 250 windows.
#> Performing PCA with 10 samples and 250 windows
#> -> pca1.
#> ------------------------------
Normalisation Method
By default, getPCA will use beta values to generate the PCA. This may be changed via the normMethod option, for instance to nrpm as shown below, although this has a tendency to show more enrichment-based variability than the use of the beta values, which correct this to some extent.
pca_beta <- getPCA(exampleTumourNormal, normMethod = "beta")
#> ------------------------------
#> Initial number of windows = 819.
#> -----------
#> Generating table of beta values for 819 regions across 10 samples.
#> -----------
#> Filtered out 220 windows with at least one missing value: 599 windows remaining.
#> ------------------------------
#> The following topVarNum values are larger than the number of remaining windows (= 599): 1000
#> These values are not used; PCA will be done with all remaining windows instead.
#> ------------------------------
#> No filtering of windows based on window standard deviation.
#> Performing PCA with 10 samples and 599 windows
#> -> pca1.
#> ------------------------------
pca_nrpm <- getPCA(exampleTumourNormal, normMethod = "nrpm")
#> ------------------------------
#> Initial number of windows = 819.
#> -----------
#> Generating table of nrpm values for 819 regions across 10 samples.
#> -----------
#> No windows have missing values.
#> ------------------------------
#> The following topVarNum values are larger than the number of remaining windows (= 819): 1000
#> These values are not used; PCA will be done with all remaining windows instead.
#> ------------------------------
#> No filtering of windows based on window standard deviation.
#> Performing PCA with 10 samples and 819 windows
#> -> pca1.
#> ------------------------------
pca_beta |> plotPCA(colour = "tissue", shape = "tumour", components = c(1,2))
#> $pca1
#> $pca1$tissue
#> $pca1$tissue$PC1vsPC2
pca_nrpm |> plotPCA(colour = "tissue", shape = "tumour", components = c(1,2))
#> $pca1
#> $pca1$tissue
#> $pca1$tissue$PC1vsPC2
Filtering regions
Whilst PCA is an unsupervised method for dimensionality reduction, sometimes there is a set of windows which are of particular interest, and you wish to see how the samples cluster when considering these. Ideally these have been selected in an orthogonal manner rather than generated using the samples that are being plotted in the PCA. For instance, differentially methylated regions may have been found using a comparison between methylation arrays (a different technology), and you wish to see how these same windows split the samples.
This may be done by pre-filtering the qseaSet object via filterByOverlaps:
exampleTumourNormal |>
filterByOverlaps(regionsToOverlap = tibble(seqnames = 7, start = 25002001, end = 25117900)) |>
getPCA()
#> ------------------------------
#> Initial number of windows = 38.
#> -----------
#> Generating table of beta values for 38 regions across 10 samples.
#> -----------
#> Filtered out 11 windows with at least one missing value: 27 windows remaining.
#> ------------------------------
#> The following topVarNum values are larger than the number of remaining windows (= 27): 1000
#> These values are not used; PCA will be done with all remaining windows instead.
#> ------------------------------
#> No filtering of windows based on window standard deviation.
#> Performing PCA with 10 samples and 27 windows
#> -> pca1.
#> ------------------------------
#> Object containing 1 dimensionality reduction objects for 10 samples
For a full list of arguments, see the documentation (?getPCA).
UMAPs
Uniform Manifold and Projection (UMAP) is another method for performing dimensionality reduction. This is a nonlinear technique, attempting to reduce the N-dimensional space onto a 2 or 3 dimensional projection, attempting to preserve both local distances between points as well as the global topology. This is more appropriate for large datasets so we use the [qsea::getExampleQseaSet()] function to create an example dataset with two types of samples, before calculating and plotting a UMAP using the getUMAP and plotUMAP functions:
set.seed(1)
qseaSetRandom <- qsea::getExampleQseaSet(repl = 20)
#> selecting windows with low CpG density for background read estimation
#> 1.48% of the windows have enrichment pattern density of at most 0.5 per fragment
#> and are used for background reads estimation
umap1 <- qseaSetRandom |>
getUMAP()
#> ------------------------------
#> Initial number of windows = 10000.
#> -----------
#> Generating table of beta values for 9800 regions across 2 samples.
#> -----------
#> Filtered out 3239 windows with at least one missing value: 6561 windows remaining.
#> ------------------------------
#> Calculating standard deviation for each window across all 40 samples:
#> Sim1T, Sim2T, Sim3T, Sim4T, Sim5T, Sim6T, Sim7T, Sim8T, Sim9T, Sim10T, Sim11T, Sim12T, Sim13T, Sim14T, Sim15T, Sim16T, Sim17T, Sim18T, Sim19T, Sim20T, Sim1N, Sim2N, Sim3N, Sim4N, Sim5N, Sim6N, Sim7N, Sim8N, Sim9N, Sim10N, Sim11N, Sim12N, Sim13N, Sim14N, Sim15N, Sim16N, Sim17N, Sim18N, Sim19N, Sim20N.
#> -> column name windowSd1.
#> ------------------------------
#> Filtering windows based on standard deviation across all 40 samples (windowSd1).
#> Standard deviation threshold = 0.203 resulting in 1000 windows.
#> Performing UMAP with 40 samples and 1000 windows
#> -> umap1.
#> ------------------------------
plotUMAP(umap1, colour = "group")
#> $umap1
#> $umap1$group
#> $umap1$group$UMAP1vsUMAP2
The vast majority of the options for getUMAP are the same as for getPCA, but while the PCA algorithm has no further options, the UMAP algorithm has a number of options which can affect the output. The underlying implementation is via the [uwot::umap()] function, and any additional parameters provided to getUMAP are passed to this function. Of particular note are the two options, n_neighbors and min_dist. n_neighbors controls how many neighbouring samples are considered to build the local structure, while min_dist controls how far apart the resulting embedded points must be. There is an excellent interactive example available to see how these two parameters influence the output, for a range of different datatypes.
Here we show two UMAPs built off the same dataset with different values of min_dist and n_neighbors:
qseaSetRandom |>
getUMAP(min_dist = 0.01, n_neighbors = 2) |>
plotUMAP(colour = "group")
#> ------------------------------
#> Initial number of windows = 10000.
#> -----------
#> Generating table of beta values for 9800 regions across 2 samples.
#> -----------
#> Filtered out 3239 windows with at least one missing value: 6561 windows remaining.
#> ------------------------------
#> Calculating standard deviation for each window across all 40 samples:
#> Sim1T, Sim2T, Sim3T, Sim4T, Sim5T, Sim6T, Sim7T, Sim8T, Sim9T, Sim10T, Sim11T, Sim12T, Sim13T, Sim14T, Sim15T, Sim16T, Sim17T, Sim18T, Sim19T, Sim20T, Sim1N, Sim2N, Sim3N, Sim4N, Sim5N, Sim6N, Sim7N, Sim8N, Sim9N, Sim10N, Sim11N, Sim12N, Sim13N, Sim14N, Sim15N, Sim16N, Sim17N, Sim18N, Sim19N, Sim20N.
#> -> column name windowSd1.
#> ------------------------------
#> Filtering windows based on standard deviation across all 40 samples (windowSd1).
#> Standard deviation threshold = 0.203 resulting in 1000 windows.
#> Performing UMAP with 40 samples and 1000 windows
#> -> umap1.
#> ------------------------------
#> $umap1
#> $umap1$group
#> $umap1$group$UMAP1vsUMAP2
qseaSetRandom |>
getUMAP(min_dist = 1, n_neighbors = 40) |>
plotUMAP(colour = "group")
#> ------------------------------
#> Initial number of windows = 10000.
#> -----------
#> Generating table of beta values for 9800 regions across 2 samples.
#> -----------
#> Filtered out 3239 windows with at least one missing value: 6561 windows remaining.
#> ------------------------------
#> Calculating standard deviation for each window across all 40 samples:
#> Sim1T, Sim2T, Sim3T, Sim4T, Sim5T, Sim6T, Sim7T, Sim8T, Sim9T, Sim10T, Sim11T, Sim12T, Sim13T, Sim14T, Sim15T, Sim16T, Sim17T, Sim18T, Sim19T, Sim20T, Sim1N, Sim2N, Sim3N, Sim4N, Sim5N, Sim6N, Sim7N, Sim8N, Sim9N, Sim10N, Sim11N, Sim12N, Sim13N, Sim14N, Sim15N, Sim16N, Sim17N, Sim18N, Sim19N, Sim20N.
#> -> column name windowSd1.
#> ------------------------------
#> Filtering windows based on standard deviation across all 40 samples (windowSd1).
#> Standard deviation threshold = 0.203 resulting in 1000 windows.
#> Performing UMAP with 40 samples and 1000 windows
#> -> umap1.
#> ------------------------------
#> $umap1
#> $umap1$group
#> $umap1$group$UMAP1vsUMAP2
As a final note, we state that UMAPs should be used with care, as they are easy to misinterpret.