various updates

NAMESPACE

@@ -74,6 +74,7 @@ export(setAnnotation)
 export(setGroup)
 export(setMap)
 export(setMetadata)
+export(setSNVID)
 export(setSlot)
 export(setTree)
 export(siteFilter)
@@ -93,6 +94,7 @@ exportMethods(getSNVID)
 exportMethods(getSlot)
 exportMethods(getTree)
 exportMethods(getTreePlot)
+exportMethods(setSNVID)
 exportMethods(setSlot)
 exportMethods(siteFilter)
 exportMethods(summary)

R/doAdonis.R

@@ -22,6 +22,7 @@
 #' @param deleteZeroDepth delete zero depth snvs (denoted as `-1`)
 #' Otherwise the cell is treated as NA
 #' @param distArg passed to `dist()`
+#' @param useMeta use metadata slot, not cl slot
 #' @param pcoa if TRUE, performs ape::pcoa on dist and plot them
 #' @importFrom vegan adonis2
 #' @importFrom stats as.formula dist
@@ -32,6 +33,7 @@ doAdonis <- function(stana, specs, cl=NULL,
     distMethod="manhattan", pcoa=FALSE,
     AAfunc=dist.ml, AAargs=list(),
     maj=FALSE, deleteZeroDepth=FALSE,
+    useMeta=FALSE,
     argList=list(), distArg=list()) {
       if (is.null(cl)) {cl <- stana@cl}
       if (pcoa) {pcoaList<-list()}
@@ -133,7 +135,15 @@ doAdonis <- function(stana, specs, cl=NULL,
             pr <- FALSE
         }
         argList[["formula"]] <- formulaPass
-        argList[["data"]] <- met
+        if (useMeta) {
+          ## Assuming all the samples are in the metadata
+          argList[["data"]] <- stana@meta[sn,]
+          if (is.null(formula)) {
+            stop("Please specify formula when using metadata slot")
+          }
+        } else {
+          argList[["data"]] <- met
+        }
 
         adores <- do.call("adonis2", argList)
 
@@ -142,6 +152,9 @@ doAdonis <- function(stana, specs, cl=NULL,
             pr <- pr[!is.na(pr)]
             r2 <- adores$R2[1]
             cat_subtle("#  F: ",adores$F[1],", R2: ",r2,", Pr: ",pr,"\n",sep="")
+        } else {
+          cat_subtle("# Printing raw adonis results ...\n")
+          print(adores)
         }
         stana@adonisList[[sp]] <- adores
 

R/doGSEA.R

@@ -174,15 +174,18 @@ plotGSEA <- function(stana, dataset_names=NULL, padjThreshold=0.05,
 #' @param discNumeric convert discrete value to numeric
 #' @param disc discretize the abundance by the threshold. function for calculating
 #' threshold, like {median}
+#' @param convert conversion such as log10
 #' @export
 #' @return stana
 addGeneAbundance <- function(stana, candSp, IDs,
     target="KO", how=sum, newCol="gene",
-    disc=NULL, discNumeric=TRUE) {
+    disc=NULL, discNumeric=TRUE, convert=NULL) {
     if (target=="KO") {
-        subMat <- stana@kos[[candSp]][IDs, ]
+        ints <- intersect(row.names(stana@kos[[candSp]]), IDs)
+        subMat <- stana@kos[[candSp]][ints, ]
     } else {
-        subMat <- stana@genes[[candSp]][IDs, ]
+        ints <- intersect(row.names(stana@genes[[candSp]]), IDs)
+        subMat <- stana@genes[[candSp]][ints, ]
     }
     if (length(IDs)>1) {
         adda <- apply(subMat, 2, how)    
@@ -200,6 +203,9 @@ addGeneAbundance <- function(stana, candSp, IDs,
         names(adda) <- nm
     }
     meta <- stana@meta
+    if (!is.null(convert)) {
+        adda <- do.call(convert, list(x=adda))
+    }
     meta[[newCol]] <- adda[row.names(stana@meta)]
     stana@meta <- meta
     return(stana)

R/inferAndPlotTree.R

@@ -40,7 +40,7 @@ inferAndPlotTree <- function(stana, species=NULL, cl=NULL,
 	target="fasta", IDs=NULL, use_point=FALSE, branch_col="black",
 	tree_args=list(), branch.length="none", point_size=2,
     subset_samples=NULL,
-    deleteZeroDepth=TRUE, treeFun="upgma", tree_only=FALSE) {
+    deleteZeroDepth=FALSE, treeFun="upgma", tree_only=FALSE) {
 	if (is.null(cl)) {cl <- stana@cl}
 	if (!is.null(meta)) {
 		meta <- checkMeta(stana, meta)
@@ -107,24 +107,27 @@ inferAndPlotTree <- function(stana, species=NULL, cl=NULL,
         	fa <- stana@fastaList[[sp]]
         	nam <- paste0(sp, "_consensus_MSA_stana_tmp.fa")
         	if (file.exists(nam)) {
-        		cat("File already exists!\n")
+        		cat_subtle("# File already exists! Removing the file ...\n")
+        		unlink(nam)
+	        	write.phyDat(fa, nam, format="fasta")
         	} else {
 	        	write.phyDat(fa, nam, format="fasta")    		
         	}
         	trenam <- paste0(sp, "_consensus_tree_stana_tmp.tree")
         	if (file.exists(trenam)) {
-        		cat("Tree file already exists!\n")
+        		cat_subtle("# Tree file already exists! Overwriting ...\n")
         	}
             system2("FastTree", args=c("-out", trenam, "-nt", nam),
                 stdout=TRUE, stderr=TRUE)
             tre <- read.tree(trenam)
         	## Call FastTree
         	## Read tree and save
-        }
-        if (treeFun=="upgma") {
+        } else if (treeFun=="upgma") {
             tre <- upgma(dm)   
-        } else {
+        } else if (treeFun=="NJ") {
             tre <- NJ(dm)
+        } else {
+        	stop("Please specify NJ, upgma, or FastTree")
         }
         if (tree_only) {
             return(tre)

R/plotDist.R

@@ -13,9 +13,10 @@
 #' @param AAargs provided to `AAfunc`
 #' @param distMethod distance method passed to dist() (default, manhattan)
 #' @param distArg passed to `dist()
+#' @param candidate list of IDs included for the calculation (like SNV IDs)
 #' @export
 plotDist <- function(stana, sp, cl=NULL, AAfunc=dist.ml, AAargs=list(),
-    target="snps", distMethod="manhattan", distArg=list()) {
+    target="snps", distMethod="manhattan", distArg=list(), candidate=NULL) {
       if (is.null(cl)) {cl <- stana@cl}
       if (length(sp)>1) {stop("Please specify one species")}
         cat_subtle("# Performing dist in ", sp, " target is ", target, "\n", sep="")
@@ -45,6 +46,12 @@ plotDist <- function(stana, sp, cl=NULL, AAfunc=dist.ml, AAargs=list(),
         } else {
             snps <- stana@genes[[sp]]
         }
+        if (!is.null(candidate)) {
+            cat_subtle("# Subset for candidate IDs\n")
+            snps <- snps[intersect(row.names(snps),candidate), ]
+        }
+        ## Make sure positions are not all NA
+        snps <- snps[,colSums(is.na(snps))!=nrow(snps)]
 
         if (!(target %in% c("tree","fasta"))) {
           distArg[["x"]] <- t(snps)

R/stana.R

@@ -190,6 +190,25 @@ setMethod("siteFilter", "stana",
 })
 
 
+#' setSNVID
+#' @param x stana object
+#' @param sp species ID
+#' @param IDs SNV ID
+#' @export
+setGeneric("setSNVID", function(x, sp, IDs) standardGeneric("setSNVID"))
+
+#' setSNVID
+#' @param x stana object
+#' @param sp species ID
+#' @param IDs SNV ID
+#' @export
+setMethod("setSNVID", "stana",
+    function(x, sp, IDs) {
+    x@includeSNVID[[sp]] <- IDs
+    x
+})
+
+
 
 #' check
 #' check and output statistics based on conditional formulas for summary