diff --git a/R/plotLoadings.R b/R/plotLoadings.R
index 3f151c59..b32f51c9 100644
--- a/R/plotLoadings.R
+++ b/R/plotLoadings.R
@@ -33,6 +33,14 @@
 #'   \item \code{absolute.scale}: ("barplot", "lollipop") \code{Logical scalar}.
 #'   Specifies whether a barplot or a lollipop plot should be visualized in
 #'   absolute scale. (Default: \code{TRUE})
+#'
+#'   \item \code{sort.method}: ("heatmap") \code{Character scalar} or
+#'   \code{function}. Specifies sorting method. If it specifies column, i.e.,
+#'   principal coordinate, features are ordered in dcreasing order based on the
+#'   specified coordinate. If \code{"neatsort"}, \code{getNeatOrder()} is
+#'   utilized. If \code{function}, the argument is passed to
+#'   \code{\link[bluster:clusterRows]{bluster:clusterRows()}} \code{BLUSPARAM}
+#'   parameter. (Default: \code{NULL})
 #' }
 #'
 #' @details
@@ -65,7 +73,7 @@
 #' plotLoadings(tse, dimred = "PCA", layout = "heatmap", add.tree = TRUE) |>
 #'     # Remove this line to see messages
 #'     suppressMessages()
-#' 
+#'
 #'
 #' # Plotting matrix as a barplot
 #' loadings_matrix <- attr(reducedDim(tse, "PCA"), "rotation")
@@ -78,7 +86,10 @@
 #' plotLoadings(loadings_matrix, layout = "heatmap")
 #'
 #' # Plot with less components
-#' plotLoadings(tse, "PCA", layout = "heatmap", ncomponents = 2)
+#' plotLoadings(
+#'     tse, "PCA", layout = "heatmap",
+#'     ncomponents = 2, sort.method = "neatsort"
+#' )
 #'
 NULL
 
@@ -230,8 +241,6 @@ setMethod("plotLoadings", signature = c(x = "matrix"),
 
 # This function manipulates the loadings data into correct format. The output
 # is data.frame in long format directly usable for ggplot.
-#' @importFrom tibble rownames_to_column
-#' @importFrom tidyr pivot_longer
 .get_loadings_plot_data <- function(df, layout, ncomponents, n = 10, ...) {
     # Transform into a dataframe
     df <- as.data.frame(df)
@@ -244,13 +253,7 @@ setMethod("plotLoadings", signature = c(x = "matrix"),
         res <- do.call(rbind, res)
     } else{
         # For heatmap, the whole data.frame is just converted into long format.
-        components <- colnames(df)
-        res <- df %>%
-            rownames_to_column(var = "Feature") %>%
-            pivot_longer(
-                cols = components,
-                names_to = "PC",
-                values_to = "Value")
+        res <- .process_heatmap_component(df, ...)
     }
     # Convert into data.frame
     res <- as.data.frame(res)
@@ -283,6 +286,64 @@ setMethod("plotLoadings", signature = c(x = "matrix"),
     return(df)
 }
 
+# This function modifies the data so that it is ready to be plotted as
+# heatmap.
+#' @importFrom tibble rownames_to_column
+#' @importFrom tidyr pivot_longer
+.process_heatmap_component <- function(
+        df, sort.method = BLUSPARAM, BLUSPARAM = NULL, ...){
+    components <- colnames(df)
+    supported_methods <- c("neatsort", components)
+    if( !(is.null(sort.method) ||
+          (.is_a_string(sort.method) && sort.method %in% supported_methods) ||
+          is.function(sort.method) ) ){
+        stop("'sort.method' must be NULL, a single character value from the ",
+            "following options '", paste0(supported_methods, collapse = "', '"),
+            "' or a BLUSPARAM function (please see bluster::clusterRows() for ",
+            "details.", call. = FALSE)
+    }
+    #
+    # Get feature order
+    feat_order <- NULL
+    if( !is.null(sort.method) ){
+        feat_order <- .get_feature_order(
+            df = df, sort.method = sort.method, ...)
+    }
+    # Put data into long format
+    res <- df |>
+        rownames_to_column(var = "Feature") |>
+        pivot_longer(
+            cols = components,
+            names_to = "PC",
+            values_to = "Value")
+    # Sort the data
+    if( !is.null(feat_order) ){
+        res[["Feature"]] <- factor(res[["Feature"]], levels = feat_order)
+    }
+    return(res)
+}
+
+# This function gets feature order based on clustering, principal component
+# values or neatsort.
+.get_feature_order <- function(df, sort.method, ...){
+    df <- as.matrix(df)
+    # Get order based on single column
+    if( sort.method %in% colnames(df) ){
+        feat_order <- order(df[[sort.method]], decreasing = TRUE)
+    } else if( sort.method %in% c("neatsort") ){
+        # Get order based on neatsort. It can take only 2 columns as input.
+        feat_order <- getNeatOrder(
+            df[, seq_len(min(2, ncol(df))), drop = FALSE], ...)
+    } else{
+        # Get order based on clustering
+        .require_package("bluster")
+        feat_order <- clusterRows(df, BLUSPARAM = sort.method) |> order()
+    }
+    # Now we have order as indices. Get the feature names in correct order.
+    feat_order <- rownames(df)[ feat_order ]
+    return(feat_order)
+}
+
 # This function calculates place for +/- sign in barplot/lollipop plot
 #' @importFrom dplyr %>% group_by mutate case_when ungroup
 .calculate_max_and_min_for_loadings <- function(df){
diff --git a/man/plotLoadings.Rd b/man/plotLoadings.Rd
index fc20edfd..11e559fd 100644
--- a/man/plotLoadings.Rd
+++ b/man/plotLoadings.Rd
@@ -37,6 +37,14 @@ Applicable when \code{layout="barplot"}. (Default: \code{10}))
 \item \code{absolute.scale}: ("barplot", "lollipop") \code{Logical scalar}.
 Specifies whether a barplot or a lollipop plot should be visualized in
 absolute scale. (Default: \code{TRUE})
+
+\item \code{sort.method}: ("heatmap") \code{Character scalar} or
+\code{function}. Specifies sorting method. If it specifies column, i.e.,
+principal coordinate, features are ordered in dcreasing order based on the
+specified coordinate. If \code{"neatsort"}, \code{getNeatOrder()} is
+utilized. If \code{function}, the argument is passed to
+\code{\link[bluster:clusterRows]{bluster:clusterRows()}} \code{BLUSPARAM}
+parameter. (Default: \code{NULL})
 }}
 
 \item{dimred}{\code{Character scalar}.  Determines the reduced dimension to
@@ -104,6 +112,9 @@ plotLoadings(tse, dimred = "PCA", ncomponents = 2, n = 12)
 plotLoadings(loadings_matrix, layout = "heatmap")
 
 # Plot with less components
-plotLoadings(tse, "PCA", layout = "heatmap", ncomponents = 2)
+plotLoadings(
+    tse, "PCA", layout = "heatmap",
+    ncomponents = 2, sort.method = "neatsort"
+)
 
 }