update NMF

DESCRIPTION

@@ -6,7 +6,7 @@ License: GPL-3
 Description: Analysis toolkit for intra-species diversity from metagenomics.
 Authors@R: person("Noriaki", "Sato", email = "nori@hgc.jp", role = c("cre", "aut"))
 Depends: ggplot2, ggstar, ggraph, igraph
-Imports: GetoptLong, BiocFileCache, RCurl, igraph, vegan, methods, data.table, phangorn, RColorBrewer, ggtree, circlize, ComplexHeatmap, ggkegg, ape, dplyr, exactRankTests, ggblend, ggh4x, scales, tidygraph, ggplotify, ggtreeExtra, ggnewscale, scico, MKmisc
+Imports: GetoptLong, BiocFileCache, RCurl, igraph, vegan, methods, data.table, phangorn, RColorBrewer, ggtree, circlize, ComplexHeatmap, ggkegg, ape, dplyr, exactRankTests, ggblend, ggh4x, scales, tidygraph, ggplotify, ggtreeExtra, ggnewscale, scico, MKmisc, NMF
 Suggests: simplifyEnrichment, stringr, tidyr, Boruta, knitr, ggrepel
 RoxygenNote: 7.3.1
 VignetteBuilder: knitr

NAMESPACE

@@ -1,7 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
+export(NMF)
 export(addGeneAbundance)
 export(alleleStat)
+export(alphaDiversityWithinSpecies)
 export(anno_PATRIC_keywords)
 export(anno_eggNOG_keywords)
 export(calcKO)
@@ -25,6 +27,7 @@ export(doL0)
 export(drawEGGNOG)
 export(drawPATRIC)
 export(exportInteractive)
+export(gene)
 export(genomeHeatmap)
 export(getAdonis)
 export(getCl)
@@ -57,10 +60,12 @@ export(setGroup)
 export(setMap)
 export(setMetadata)
 export(setTree)
+export(snv)
 export(strainClusterHeatmap)
 export(summariseAbundance)
 import(BiocFileCache)
 import(GetoptLong)
+import(NMF)
 import(RColorBrewer)
 import(RCurl)
 import(ggplot2)

R/NMF.R

@@ -0,0 +1,39 @@
+
+#' NMF
+#' @param stana stana object
+#' @param species candidate species ID
+#' @param rank rank of NMF
+#' @param target KO, gene or snv, default to KO
+#' @param seed random seed
+#' @param method NMF method, default to snmf/r
+#' @import NMF
+#' @export
+NMF <- function(stana, species, rank=10, target="KO", seed=53, method="snmf/r",
+    deleteZeroDepth=TRUE, beta=0.01) {
+    if (target=="KO") {
+        mat <- stana@kos[[species]]
+    } else if (target=="genes") {
+        mat <- stana@genes[[species]]
+    } else if (target=="snps") {
+        mat <- stana@snps[[species]]
+        if (deleteZeroDepth) {
+           mat <- mat[rowSums(mat == -1)==0,]
+           qqcat("After filtering `-1`, position numbers: @{dim(mat)[1]}\n")
+        }
+    }
+    res <- NMF::nmf(mat, rank = rank, seed = seed, method=method, beta=beta)
+    stana@NMF[[species]] <- res
+    return(stana)
+}
+
+#' alphaDiversityWithinSpecies
+#' @export
+alphaDiversityWithinSpecies <- function(stana, species) {
+    if (is.null(stana@NMF[[species]])) {
+        stana <- NMF(stana, species)
+    }
+    res <- stana@NMF[[species]]
+    W <- coef(res)
+    div <- vegan::diversity(t(W))
+    return(div)
+}

R/stana.R

@@ -21,6 +21,7 @@ setClass("stana", slots=list(
                             strainClusters="data.table",
                             snps="list",
                             eggNOG="list",
+                            NMF = "list",
                             kos="list",
                             snpsInfo="list",
                             snpsDepth="list",
@@ -238,6 +239,39 @@ changeColors <- function(stana, colors) {
 	return(stana)
 }
 
+#' snv
+#' @param dfs named list of dfs
+#' @export
+snv <- function(dfs) {
+    stana <- new("stana")
+    stana@type <- "manual"
+
+    if (is.null(names(dfs))) {
+        ids <- paste0("species", seq_len(length(dfs)))
+    } else {
+        ids <- names(dfs)
+    }
+    stana@snps <- dfs
+    stana@ids <- ids
+    return(stana)
+}
+
+#' gene
+#' @param dfs named list of dfs
+#' @export
+gene <- function(dfs) {
+    stana <- new("stana")
+    stana@type <- "manual"
+    if (is.null(names(dfs))) {
+        ids <- paste0("species", seq_len(length(dfs)))
+    } else {
+        ids <- names(dfs)
+    }
+    stana@genes <- dfs
+    stana@ids <- ids
+    return(stana)
+}
+
 #' getGenes (concatenated to checkProfile)
 #' 
 #' Obtain gene matrix from midas merge directory from MIDAS.