update, not to use lfc by default

R/doGSEA.R

@@ -1,4 +1,5 @@
 #' doGSEA
+#' @return GSEA results from clusterProfiler
 #' @export
 doGSEA <- function(stana, candSp=NULL, cl=NULL, eps=1e-2, how=mean,
     zeroPerc=0) {
@@ -27,11 +28,8 @@ doGSEA <- function(stana, candSp=NULL, cl=NULL, eps=1e-2, how=mean,
     ko_df_filt <- data.frame(ko_df_filt[!rowSums(ko_df_filt==0)>dim(ko_df_filt)[2] * zeroPerc, ]) |>
         `colnames<-`(colnames(ko_df_filt))
 
-
     ## eps values
-    ko_df_filt <- ko_df_filt + 1e-2
-    ko_sum <- log2(apply(ko_df_filt[,intersect(inSample, aa)],1,mean) /
-                   apply(ko_df_filt[,intersect(inSample, bb)],1,mean))
+    ko_sum <- L2FC(ko_df_filt, aa, bb, method="t", eps=eps)
     ## Perform GSEA (it will take time)?
     ko_sum <- ko_sum[order(ko_sum, decreasing=TRUE)]
 
@@ -46,6 +44,29 @@ doGSEA <- function(stana, candSp=NULL, cl=NULL, eps=1e-2, how=mean,
     return(enr)
 }
 
+#' @param mat row corresponds to gene, and column samples
+#' @noRd
+L2FC <- function(mat, l1, l2, method="t", eps=0) {
+    if (method == "gmean") {
+        l1_mean <- apply(log2(mat[, intersect(colnames(mat), l1)] + eps), 1, mean)
+        l2_mean <- apply(log2(mat[, intersect(colnames(mat), l2)] + eps), 1, mean)
+        return(l1_mean - l2_mean)
+    } else if (method == "t") {
+        res <- lapply(row.names(mat), function(m) {
+            tres <- t.test(mat[m, intersect(colnames(mat), l1)],
+                mat[m, intersect(colnames(mat), l2)])
+            as.numeric(tres$statistic)
+        }) %>% unlist()
+        names(res) <- row.names(mat)
+        return(res)
+    } else {
+        l1_mean <- apply(mat[, intersect(colnames(mat), l1)], 1, mean)
+        l2_mean <- apply(mat[, intersect(colnames(mat), l2)], 1, mean)
+        return(log2((l1_mean+eps) / (l2_mean+eps)))
+    }
+}
+
+
 #' calcKO
 #' @export
 calcKO <- function(stana, candSp=NULL, how=mean) {

R/plotKEGGPathway.R

@@ -76,17 +76,13 @@ plotKEGGPathway <- function(stana, species, pathway_id,
     if (summarize) {
         commons <- Reduce(intersect, lapply(stana@kos[species], function(x) row.names(x)))
         commoncol <- Reduce(intersect, lapply(stana@kos[species], function(x) colnames(x)))
-        ko_tbl <- Reduce("+", lapply(stana@kos[c("101346","102438")], function(x) x[commons,commoncol]))
+        ko_tbl <- Reduce("+", lapply(stana@kos[species], function(x) x[commons,commoncol]))
 
         if (sum_flag) {
           lfcs[["Sum"]] <- apply(ko_tbl, 1, sum)
         } else {
             qqcat("@{sp}: @{names(cl)[1]} / @{names(cl)[2]}\n")
-            ko_tbl <- ko_tbl + eps
-            inc1 <- intersect(colnames(ko_tbl), cl[[1]])
-            inc2 <- intersect(colnames(ko_tbl), cl[[2]])
-            lfcs[["Sum"]] <- log2((apply(ko_tbl[,inc1], 1, how)) / 
-                                 (apply(ko_tbl[,inc2], 1, how)))
+            lfc[[sp]] <- L2FC(ko_tbl, cl[[1]], cl[[2]])
         }
     } else {
         for (sp in species) {
@@ -95,11 +91,7 @@ plotKEGGPathway <- function(stana, species, pathway_id,
                 lfcs[[sp]] <- apply(ko_tbl, 1, sum)
             } else {
                 qqcat("@{sp}: @{names(cl)[1]} / @{names(cl)[2]}\n")
-                ko_tbl <- ko_tbl + eps
-                inc1 <- intersect(colnames(ko_tbl), cl[[1]])
-                inc2 <- intersect(colnames(ko_tbl), cl[[2]])
-                lfcs[[sp]] <- log2((apply(ko_tbl[,inc1], 1, how)) / 
-                                (apply(ko_tbl[,inc2], 1, how)))
+                lfcs[[sp]] <- lL2FC(ko_tbl, cl[[1]], cl[[2]])
             }
         }        
     }

R/plotTree.R

@@ -59,7 +59,8 @@ plotTree <- function(stana, species=NULL, cl=NULL,
 			tree_args[["x"]] <- t(mat)
 			dm <- do.call(dist, tree_args)			
 		}
-		tre <- NJ(dm)
+        tre <- NJ(dm)
+
 		stana@treeList[[sp]] <- tre			
 
 		## Plot tree		
@@ -90,10 +91,8 @@ plotTree <- function(stana, species=NULL, cl=NULL,
                     meta[[tmp_show_cv]] <- as.factor(change_cv)
                 }
             }
-
             meta <- meta[tre$tip.label,]
 
-
             if (layout=="rectangular") {
                 qqcat("Force setting layout to 'circular'\n")
                 flyt <- "circular"
@@ -128,14 +127,14 @@ plotTree <- function(stana, species=NULL, cl=NULL,
 	                        data = meta,
 	                        geom = geom_point,
 	                        size = point_size, shape=21, 
-	                        mapping = aes(y=id, x=0.5, fill=.data[[tmp_show_cv]])
+	                        mapping = aes(y=id, fill=.data[[tmp_show_cv]])
 	                    )                	
 	                } else {
 	                    p <- p + geom_fruit(
 	                        data = meta,
 	                        geom = geom_star,
 	                        size = point_size,
-	                        mapping = aes(y=id, x=0.5, fill=.data[[tmp_show_cv]]),
+	                        mapping = aes(y=id, fill=.data[[tmp_show_cv]]),
 	                        starshape = starsh[tmp_show_cv]
 	                    )
 	                }