From 74e675cb88c96f9d3f1083e698c2ee96a61b8990 Mon Sep 17 00:00:00 2001
From: philipperuiz <philippe.ruiz@inrae.fr>
Date: Wed, 13 Dec 2023 15:37:20 +0100
Subject: [PATCH 1/3] a function to createheatmap from deseq2 output
---
R/heatmapDeseq2_fun.R | 84 +++++++++++++++++++++++++++++++++++++++++++
1 file changed, 84 insertions(+)
create mode 100644 R/heatmapDeseq2_fun.R
diff --git a/R/heatmapDeseq2_fun.R b/R/heatmapDeseq2_fun.R
new file mode 100644
index 0000000..bebd381
--- /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)
+}
--
GitLab
From 41efa5db5639520cb8c455bc144f04c3ac8af5c1 Mon Sep 17 00:00:00 2001
From: philipperuiz <philippe.ruiz@inrae.fr>
Date: Fri, 15 Dec 2023 11:26:34 +0100
Subject: [PATCH 2/3] update namespace and plsda
---
NAMESPACE | 3 +++
man/heatmapDeseq2_fun.Rd | 23 +++++++++++++++++++++++
man/plsda_fun.Rd | 22 +++++++++++++++++++++-
3 files changed, 47 insertions(+), 1 deletion(-)
create mode 100644 man/heatmapDeseq2_fun.Rd
diff --git a/NAMESPACE b/NAMESPACE
index ee9b190..388d556 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/man/heatmapDeseq2_fun.Rd b/man/heatmapDeseq2_fun.Rd
new file mode 100644
index 0000000..c93748b
--- /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
+}
diff --git a/man/plsda_fun.Rd b/man/plsda_fun.Rd
index 475a46f..c0be210 100644
--- a/man/plsda_fun.Rd
+++ b/man/plsda_fun.Rd
@@ -4,7 +4,17 @@
\alias{plsda_fun}
\title{PLSDA from MixOmics package}
\usage{
-plsda_fun(data = data, output = "./plsda/", column1 = "", rank = "ASV")
+plsda_fun(
+ data = data,
+ output = "./plsda/",
+ column1 = "",
+ rank = "ASV",
+ axis = c(1, 2),
+ multilevel = NULL,
+ ind.names = FALSE,
+ ellipse = FALSE,
+ progressBar = FALSE
+)
}
\arguments{
\item{data}{A phyloseq object (output from decontam or generate_phyloseq)}
@@ -14,6 +24,16 @@ plsda_fun(data = data, output = "./plsda/", column1 = "", rank = "ASV")
\item{column1}{Factor to test.}
\item{rank}{Taxonomy rank to merge features that have same taxonomy at a certain taxonomic rank (among rank_names(data), or 'ASV' for no glom)}
+
+\item{axis}{Select the axis to plot}
+
+\item{multilevel}{Use Within matrix decomposition for repeated measurements}
+
+\item{ind.names}{either a character vector of names for the individuals to be plotted, or FALSE for no names. If TRUE, the row names of the first (or second) data matrix is used as names (see mixOmics Details).}
+
+\item{ellipse}{Logical indicating if ellipse plots should be plotted (see mixOmics Details).}
+
+\item{progressBar}{Silence the progress bar.}
}
\value{
Return a list object with plots, and loadings table.
--
GitLab
From a50440461ffbe2245c99fa717c587b238fb22547 Mon Sep 17 00:00:00 2001
From: philipperuiz <philippe.ruiz@inrae.fr>
Date: Wed, 3 Apr 2024 15:23:33 +0200
Subject: [PATCH 3/3] =?UTF-8?q?plsda=20tune=20return=20a=20warning=20if=20?=
=?UTF-8?q?only=202=C2=A0components=20are=20selected=20and=20don't=20stop?=
=?UTF-8?q?=20during=20plotVar?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
R/plsda_fun.r | 5 ++++-
1 file changed, 4 insertions(+), 1 deletion(-)
diff --git a/R/plsda_fun.r b/R/plsda_fun.r
index a5599ce..d9bd4a1 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()
--
GitLab