diff --git a/NAMESPACE b/NAMESPACE
index ee9b190ce48213f52a426cd3a8d906283608bb8d..388d55687355e4c2197493164e35254e05452eb0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -18,6 +18,7 @@ export(export_to_stamp_fun)
 export(fill_tax_fun)
 export(generate_phyloseq_fun)
 export(generate_tree_fun)
+export(heatmapDeseq2_fun)
 export(heatmap_fun)
 export(idTaxa_assign)
 export(idtaxa_assign_fasta_fun)
@@ -49,12 +50,14 @@ import(ggplot2)
 import(here)
 import(metagenomeSeq)
 import(mixOmics)
+import(pheatmap)
 import(phyloseq)
 import(psadd)
 import(readr)
 import(readxl)
 import(reshape2)
 import(scales)
+import(stringr)
 import(tools)
 import(vegan)
 import(vroom)
diff --git a/R/heatmapDeseq2_fun.R b/R/heatmapDeseq2_fun.R
new file mode 100644
index 0000000000000000000000000000000000000000..bebd381516734e8dd263efa8221b4190c8c44935
--- /dev/null
+++ b/R/heatmapDeseq2_fun.R
@@ -0,0 +1,84 @@
+#' heatmap for deseq2_fun output
+#'
+#' @param desq output of deseq2_fun
+#' @param phys the phyloseq object used for deseq2 analyse
+#' @param var Factor to test.
+#' @param workingRank Taxonomy rank to merge features that have same taxonomy at a certain taxonomic rank (among rank_names(data), or 'ASV' for no glom)
+#' @return Return a list object with heatmap plots.
+#' @import futile.logger
+#' @import stringr
+#' @import phyloseq
+#' @import pheatmap
+#' @export
+heatmapDeseq2_fun <- function(desq, phys, var, workingRank) {
+  pheat = list()
+  for (i in names(desq)) {
+    tableDESeqSig <- desq[[i]][["table"]] |> subset(padj < 0.05) |> na.omit()
+
+    # remove unused taxa
+    psTaxaRelSig <- prune_taxa(row.names(tableDESeqSig),
+                               transform_sample_counts(phys, function(x) x/sum(x)))
+
+    # remove unused sample
+    usedName <<- str_split_fixed(names(desq[[i]]$plot$data)[1], "_", 2)
+    keepSamp <<- str_split(usedName[2], "_vs_", simplify = TRUE)
+
+    psTaxaRelSig2 <- paste("psTaxaRelSigSubset <<- subset_samples(psTaxaRelSig, ",
+                           var, "%in%
+                                                   keepSamp)", sep = "")
+    eval(parse(text = psTaxaRelSig2))
+    psTaxaRelSigSubset <- subset_taxa(psTaxaRelSigSubset,
+                                      taxa_sums(psTaxaRelSigSubset) != 0)
+    matrix <- otu_table(psTaxaRelSigSubset) |>
+      data.frame() |>
+      as.matrix()
+    if (!taxa_are_rows(psTaxaRelSigSubset)) {
+      colnames(matrix) <- tax_table(psTaxaRelSigSubset)[, workingRank] |>
+        as.character()
+    } else {
+      rownames(matrix) <- tax_table(psTaxaRelSigSubset)[, workingRank] |>
+        as.character()
+    }
+
+    metadataSub <- sample_data(psTaxaRelSigSubset) |>
+      data.frame()
+    # print(metadataSub)
+
+    if (length(var) == 1) {
+      annotationCol <- data.frame(
+        var1 = metadataSub[, var[1]] |>
+          as.vector() |>
+          as.factor(),
+        check.names = FALSE)
+    } else {
+      annotationCol <- data.frame(
+        var1 = metadataSub[, var[1]] |>
+          as.vector() |>
+          as.factor(),
+        var2 = metadataSub[, var[2]] |>
+          as.vector() |>
+          as.factor(),
+        check.names = FALSE)
+    }
+
+    rownames(annotationCol) <- rownames(metadataSub)
+    colnames(annotationCol) <- var
+
+    annotationRow <- data.frame(
+      Phylum = tax_table(psTaxaRelSigSubset)[, "Phylum"] |>
+        as.factor()
+    )
+    if (taxa_are_rows(psTaxaRelSigSubset)) {
+      rownames(annotationRow) <- rownames(matrix)
+    } else {
+      rownames(annotationRow) <- colnames(matrix)
+    }
+    title <- usedName[,2]
+
+    pheat[[i]] <- pheatmap(matrix, scale = "row",
+                           annotationCol = annotationCol,
+                           annotationRow = annotationRow,
+                           main = title)
+  }
+  return(pheat)
+}
diff --git a/R/plsda_fun.r b/R/plsda_fun.r
index a5599ce41c3f6f55f8ab212bcd906429c67679d9..d9bd4a1e020ef537265c3915420cd2f021007948 100644
--- a/R/plsda_fun.r
+++ b/R/plsda_fun.r
@@ -139,6 +139,9 @@ plsda_fun <- function (data = data, output = "./plsda/", column1 = "",
   }
   r_list <- tune_splsda()
   ncomp <- r_list$ncomp
+  if(r_list$ncomp < 3) {
+    warning("sPLS-DA (regression mode) has less than 3 sPLS-DA components. ")
+  }
   select.keepX <- r_list$selectkeepX
   flog.info(paste("keepX: ", select.keepX, sep = ""))
   flog.info("SPLSDA...")
@@ -199,7 +202,7 @@ plsda_fun <- function (data = data, output = "./plsda/", column1 = "",
   }
 
   outF[["var"]] = list()
-  for (comp in 2:3) {
+  for (comp in 2:splsda.res$ncomp) {
     plotVar(splsda.res, comp = c(1,comp), title = paste("Loadings on comp 1 and ",
                                                         comp, sep = ""))
     outF$var[[glue::glue("comp{comp}")]] <- recordPlot()
diff --git a/man/heatmapDeseq2_fun.Rd b/man/heatmapDeseq2_fun.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..c93748b174b6f1b69adc5f6eeafda34437dd655e
--- /dev/null
+++ b/man/heatmapDeseq2_fun.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/heatmapDeseq2_fun.R
+\name{heatmapDeseq2_fun}
+\alias{heatmapDeseq2_fun}
+\title{heatmap for deseq2_fun output}
+\usage{
+heatmapDeseq2_fun(desq, phys, var, workingRank)
+}
+\arguments{
+\item{desq}{output of deseq2_fun}
+
+\item{phys}{the phyloseq object used for deseq2 analyse}
+
+\item{var}{Factor to test.}
+
+\item{workingRank}{Taxonomy rank to merge features that have same taxonomy at a certain taxonomic rank (among rank_names(data), or 'ASV' for no glom)}
+}
+\value{
+Return a list object with heatmap plots.
+}
+\description{
+heatmap for deseq2_fun output
+}