updated the readme file

.Rbuildignore

@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,3 @@
+.Rproj.user
+.Rhistory
+.RData

DESCRIPTION

@@ -0,0 +1,25 @@
+Package: directPA
+Type: Package
+Title: Direction Analysis for Pathways and Kinases
+Version: 1.5.1
+Date: 2023-11-14
+Author: Pengyi Yang & Ellis Patrick
+Maintainer: Pengyi Yang <[email protected]>
+Description: Direction analysis is a set of tools designed to identify
+    combinatorial effects of multiple treatments/conditions on pathways
+    and kinases profiled by microarray, RNA-seq, proteomics, or phosphoproteomics
+    data. See Yang P et al (2014) <doi:10.1093/bioinformatics/btt616>; and Yang P et al. (2016) <doi:10.1002/pmic.201600068>.
+License: GPL-3
+Depends:
+    R (>= 3.10.0),
+Imports:
+    grDevices,
+    graphics,
+    stats,
+    plotly,
+    calibrate
+Packaged: 2016-03-22 13:08:32 UTC; Pengyi
+NeedsCompilation: no
+Repository: CRAN
+Date/Publication: 2016-03-22 18:00:22
+RoxygenNote: 7.1.1

MD5

@@ -0,0 +1,30 @@
+e05721202d78707adb25b5dadb9a905c *DESCRIPTION
+32f7448fafd8f272e98f8f3053cc3f51 *NAMESPACE
+61911534dedf36247467a549fbf57ec1 *R/directEnrichPlot2d.R
+76cad031899e402b3ab2c8ee3bad588f *R/directExplorer2d.R
+c193001c14dbc241fd64be4f11202254 *R/directPA.R
+016bae5eb003dc005c4dd37ccc50f58e *R/geneStats.R
+c94800bbd93527de53b09fc6bbce19a3 *R/pathwayStats.R
+0fb8763a7e2602a1004522f451bbea85 *R/rotate2d.R
+989f12754c30a855e238dd4a160a9806 *R/rotate3d.R
+83648a7a500619337f5f19913ad14f71 *data/HEK.rda
+079cc541e6cfd73ded8ab5b20a8687e6 *data/PM.rda
+d8b015f9de0c68ca738fa61b8604add3 *data/Pathways.rda
+a4829f119c5d601a1663ff68b6f033c9 *data/PhosphoELM.rda
+48bf3aaf03017c2519d62984ea9587a2 *data/PhosphoSite.rda
+a8abd01dded8fe74481dfc5d6ee6284a *man/HEK.Rd
+c02262b58e37b1fefd695441707b8809 *man/PM.Rd
+19d4e0276d10cd41f6d577ed5245ea12 *man/Pathways.KEGG.Rd
+ffd89f58f9af47d0ff556011a5c3521c *man/Pathways.reactome.Rd
+50f5b7321b62385a4af7e5c1ee1ee901 *man/PhosphoELM.human.Rd
+8e9a9b800cbdb9e5e9764143c54ef6fb *man/PhosphoELM.mouse.Rd
+9224f661163899dbf13b3b3a8775a167 *man/PhosphoSite.human.Rd
+02df716ef54fb50fdecd8648c6883235 *man/PhosphoSite.mouse.Rd
+3d1ba4de7361b2947991bf8b6f1de90b *man/directEnrichPlot2d.Rd
+29b9f292ce686fa8368ab99a2a527a74 *man/directExplorer2d.Rd
+8809fdbeedb03dd0c8672be14ef11383 *man/directPA-package.Rd
+02069d7456d6905cc361baab683de7f8 *man/directPA.Rd
+fcc77f5e75cdbce83340ec66f2777f47 *man/geneStats.Rd
+aefc1a7e9de81ead1d8f8d1911002812 *man/pathwayStats.Rd
+3e13753f1c19b7fb239f048db2376b20 *man/rotate2d.Rd
+574ba8d8e806479f12b078b52b6de344 *man/rotate3d.Rd

NAMESPACE

@@ -0,0 +1,15 @@
+# Generated by roxygen2: do not edit by hand
+
+export(directExplorer2d)
+export(directPA)
+export(geneStats)
+export(kinasePA)
+export(pathwayStats)
+export(perturbPlot2d)
+export(perturbPlot3d)
+export(rotate2d)
+export(rotate3d)
+import(calibrate)
+import(grDevices)
+import(graphics)
+import(stats)

source/directPA_1.1/R/directExplorer2d.R → R/directExplorer2d.R

@@ -1,7 +1,8 @@
 #' Batch Direction Analysis in 2-dimentional space
 #'
 #' Rotate to the direction of interest in polar coordinates by degree (e.g. pi/4).
-#' @usage directExplorer2d(Tc, annotation=NULL, gene.method="OSP", path.method="Stouffer", top=10, ...)
+#' @usage directExplorer2d(Tc, annotation=NULL, gene.method="OSP", 
+#' path.method="Stouffer", top=10, nd=8, ...)
 #' @param Tc a numeric matrix with 2 columns. The rows are genes or phosphorylation sites and the columns 
 #' are treatments vs control statistics.
 #' @param annotation a list with names correspond to pathways or kinases and elements correspond to 
@@ -13,6 +14,7 @@
 #' sites that belongs to a pathway or kinase. Available methods are Stouffer, OSP, Fisher, and maxP. 
 #' Default method is Stouffer. 
 #' @param top the number of entries to be highlighted in the plot.
+#' @param nd the number of directions to plot (4 or 8)
 #' @param ... parameters for controlling the plot.
 #' @return The the list of enrichment analysis in tables. 
 #' @export
@@ -41,26 +43,44 @@
 #' bda$gene.tab[order(bda$gene.tab[,"--"]),][1:10,]
 #' bda$path.tab[order(bda$path.tab[,"--"]),][1:10,]
 #'
-directExplorer2d <- function(Tc, annotation=NULL, gene.method="OSP", path.method="Stouffer", top=10, ...) {
+directExplorer2d <- function(Tc, annotation=NULL, gene.method="OSP", path.method="Stouffer", top=10, nd=8, ...) {
 
   directionCode <- c('++','+*','+-','*-','--','-*','-+','*+')
-  gene.tab <- matrix(NA, nrow(Tc), 8) 
-  colnames(gene.tab) <- directionCode
+  gene.tab <- matrix(NA, nrow(Tc), nd)
   rownames(gene.tab) <- rownames(Tc)
-
+  if(nd != 8) {
+    colnames(gene.tab) <- directionCode[c(1,3,5,7)]
+  } else {
+    colnames(gene.tab) <- directionCode
+  }
+
   path.tab <- NULL
   if (!is.null(annotation)) {
-    path.tab <- matrix(NA, length(annotation), 9)
-    colnames(path.tab) <- c("size", directionCode)
-    rownames(path.tab) <- names(annotation)   
+    path.tab <- matrix(NA, length(annotation), nd+1)
+    rownames(path.tab) <- names(annotation)
+    if(nd != 8) {
+      colnames(path.tab) <- c("size", directionCode[c(1,3,5,7)])
+    } else {
+      colnames(path.tab) <- c("size", directionCode)
+    }
+  }
+
+  ds <- c()
+  if (nd != 8) {
+    ds <- c(0, 2, 4, 6)
+  } else {
+    ds <- 0:7
   }
+  count <- 0
 
-  for(i in 0:7){
+  for(i in ds){
+    count <- count + 1
+
     # gene or phosphorylation site level
     Tc.zscores <- apply(Tc, 2, function(x){qnorm(rank(x)/(nrow(Tc)+1))})
     Tc.rotated <- rotate2d(Tc.zscores, pi/4*i)
     gene.pvalues <- apply(Tc.rotated, 1, geneStats, gene.method)
-    gene.tab[,(i+1)] <- gene.pvalues
+    gene.tab[,count] <- gene.pvalues
 
 	  # pathway or kinase level
     if (!is.null(annotation)) {
@@ -69,19 +89,21 @@ directExplorer2d <- function(Tc, annotation=NULL, gene.method="OSP", path.method
 	    if (i == 0) {
 		    path.tab[,1] <- unlist(gst[,"size"])
 	    }
-	    path.tab[,(i+2)] <- unlist(gst[,"pvalue"])
+	    path.tab[,(count+1)] <- unlist(gst[,"pvalue"])
     }
   }
 
   # highlight candidates
   plot(Tc, col="gray", pch=16, ...)
   abline(h=0, v=0, col="gold", lty=2)
   abline(a=0, b=1, lty = 2, col="darkgreen")
-  color <- rainbow(8)
-  for (i in 1:8) {
-     ids <- names(sort(gene.tab[,i])[1:top])
-     points(Tc[ids,], col=color[i], pch=16)
-     textxy(Tc[ids,1], Tc[ids,2], ids, col = color[i])
+  color <- c("red", "orange3", "green4", "#169B48", "blue4", "#0F8AB5", "purple", "gray50")
+  count <- 0
+  for (i in ds) {
+    count <- count + 1
+     ids <- names(sort(gene.tab[,count])[1:top])
+     points(Tc[ids,], col=color[i+1], pch=16)
+     textxy(Tc[ids,1], Tc[ids,2], ids, col = color[i+1])
   }
 
   results <- list()

R/directPA.R

@@ -0,0 +1,129 @@
+#' Direction Analysis for Pathways
+#'
+#' The main function of direction Analysis. This function takes in a matrix of test statistics with 
+#' two (2-dimensional space) or three (3-dimensional space) columns, the direction of interests, and 
+#' the annotation list such as pathway annotation, and test for enrichment of pathways on the specified 
+#' direction.
+#' 
+#' @usage directPA(Tc, direction, annotation, minSize=5, gene.method="OSP", 
+#' path.method="Stouffer", visualize=TRUE, ...)
+#' 
+#' @param Tc a numeric matrix. Rows are genes and columns are treatments vs control statistics.
+#' @param direction the direction to be tested for enrichment. Either specified as a degree for 
+#' two-dimensional analysis or as contrast (in a triplet) for three-dimensional analysis.
+#' @param annotation a list with names correspond to pathways and elements correspond to genes belong to 
+#' each pathway, respectively.
+#' @param minSize the size of annotation groups to be considered for calculating enrichment. Groups 
+#' that are smaller than the minSize will be removed from the analysis.
+#' @param gene.method the method to be used for integrating statistics across treatments for each gene. 
+#' Available methods are Stouffer, OSP, Fisher, and maxP. 
+#' Default method is OSP.
+#' @param path.method the method to be used for integrating statistics of all genes that belongs to a 
+#' pathway. Available methods are Stouffer, OSP, Fisher, and maxP. 
+#' Default method is Stouffer.
+#' @param visualize whether to visualize the plot.
+#' @param ... other visualization parameters to pass on.
+#' @return a list that contains directional p-values for each gene and directional enrichment for each pathway.
+#' @export
+#' @examples
+#' 
+#' # load the proteomics dataset
+#' data(PM)
+#' 
+#' # load pathway annotations
+#' data(Pathways)
+#' 
+#' # display reactome pathways. Could be replaced by any other pathway databases
+#' Pathways.reactome[1:5]
+#' 
+#' # direction pathway analysis in 3-dimensional space. Implemnted as rotating by contrast
+#' # (1) test combined effect of all 3 treatments (stimulation and inhibitions) vs control (basal) 
+#' # on the original direction.
+#' dPA <- directPA(Tc=PM, direction=c(1,1,1), annotation=Pathways.reactome)
+#' dPA$gst[order(unlist(dPA$gst[,1])),][1:20,]
+#' # rank substrates on the direciton of interest
+#' sort(dPA$gene.pvalues)[1:20]
+#' 
+#' # (2) test combined effect of all 3 treatments vs controls on direction c(1,-1, 0)
+#' # this rotates Ins by 0 degree, Wmn by 90 degree, and MK by 45 degree.
+#' dPA <- directPA(Tc=PM, direction=c(1,-1,0), annotation=Pathways.reactome)
+#' dPA$gst[order(unlist(dPA$gst[,1])),][1:20,]
+#' 
+#' # (3) test combined effect of all 3 perturbations vs controls on direction c(1,-1, 1)
+#' # this rotates Ins by 0 degree, Wmn by 90 degree, and MK by 0 degree.
+#' dPA <- directPA(Tc=PM, direction=c(1,-1,1), annotation=Pathways.reactome)
+#' dPA$gst[order(unlist(dPA$gst[,1])),][1:20,]
+#'
+directPA <- function(Tc, direction, annotation, minSize=5, gene.method="OSP", path.method="Stouffer", visualize=TRUE, ...){
+  ## spherical coordinates for three-dimensional rotation
+  if (length(direction) == 3) {
+
+    # step 1. convert statistics into z-scores
+    Tc.zscores <- apply(Tc, 2, function(x){qnorm(rank(x)/(nrow(Tc)+1))})
+
+    # step 2. rotate z-scores
+    Tc.rotated <- rotate3d(Tc.zscores, direction)
+
+    # step 3. integrate statistics across treatments
+    gene.pvalues <- apply(Tc.rotated, 1, geneStats, gene.method)
+
+    if (visualize == TRUE) {
+      #HC = rainbow(length(gene.pvalues)*1.2)
+      #plot3d(Tc, col=HC[rank(gene.pvalues)], size=5, ...)
+      #abclines3d(x=0, y=0, z=0, a=diag(3), col="black", lwd=3)
+      #abclines3d(x=0, a=direction, col="pink", lwd=5)
+
+      df <- data.frame(Tc)
+      colnames(df) <- c("x", "y", "z")
+      df$pvalue <- gene.pvalues[rownames(df)]
+
+      my_col = colorRampPalette(rainbow(12))(100)
+
+      p <- plotly::plot_ly(df, x=~x, y=~y, z=~z, color=~pvalue, colors=my_col, size=5)
+      p <- plotly::add_markers(p) 
+      p <- plotly::layout(p, scene = list(xaxis = list(title = colnames(Tc)[[1]]),
+                                     yaxis = list(title = colnames(Tc)[[2]]),
+                                     zaxis = list(title = colnames(Tc)[[3]])))
+      print(p)
+    }
+
+    # step 4. integrate statistics for pathways
+    gene.zscores <- qnorm(gene.pvalues, lower.tail = FALSE)
+    gst <- t(sapply(annotation, pathwayStats, gene.zscores, minSize=5, path.method))
+
+    result <- list()
+    result$gene.pvalues <- gene.pvalues
+    result$gst <- gst
+    return(result)
+  }
+
+  ## polar coordinates for two-dimensional rotation
+  if (length(direction) == 1) {
+
+    # step 1. convert statistics into z-scores
+    Tc.zscores <- apply(Tc, 2, function(x){qnorm(rank(x)/(nrow(Tc)+1))})
+
+    # step 2. rotate z-scores
+    Tc.rotated <- rotate2d(Tc.zscores, direction)
+
+    # step 3. integrate statistics across treatments
+    gene.pvalues <- apply(Tc.rotated, 1, geneStats, gene.method)
+
+    if (visualize == TRUE) {
+      HC = rainbow(length(gene.pvalues)*1.2)
+      plot(Tc, col=HC[rank(gene.pvalues)], pch=16, ...)
+      abline(v = 0,h = 0, lty=2, col="gold")
+      abline(a=0, b=1, col="darkgreen", lty=2)
+      abline(a=0, b=-1, col="darkgreen", lty=2)
+    }
+
+    # step 4. integrate statistics for pathways
+    gene.zscores <- qnorm(gene.pvalues, lower.tail = FALSE)
+    gst <- t(sapply(annotation, pathwayStats, gene.zscores, minSize=5, path.method))
+
+    result <- list()
+    result$gene.pvalues <- sort(gene.pvalues)
+    result$pathways <- gst[order(as.numeric(gst[,1])),]
+    return(result)
+  }
+}

source/directPA_1.1/R/geneStats.R → R/geneStats.R

@@ -1,4 +1,4 @@
-#' Gene Level Statistics
+#' Molecule Level Statistics
 #' 
 #' Takes a vector of statistics with each element corresponds to a treatment vs control comparison, 
 #' and calculates a combined statistics accross multiple treatments.

R/kinasePA.R

@@ -0,0 +1,59 @@
+#' Direction Analysis for Kinases
+#'
+#' This is a wrapper for runing directPA for kinase perturbation analysis (kinasePA)
+#' 
+#' @usage kinasePA(Tc, direction, annotation, minSize=5, substrate.method="OSP", 
+#' kinase.method="Stouffer", visualize=TRUE, ...)
+#' 
+#' @param Tc a numeric matrix. The columns are phosphorylation sites and the columns are treatments vs 
+#' control statistics.
+#' @param direction the direction to be tested for enrichment. Either specified as a degree for 
+#' two-dimensional analysis or as contrast (in a triplet) for three-dimensional analysis.
+#' @param annotation a list with names correspond to kinases and elements correspond to substrates belong 
+#' to each kinase, respectively.
+#' @param minSize the size of annotation groups to be considered for calculating enrichment. Groups 
+#' that are smaller than the minSize will be removed from the analysis.
+#' @param substrate.method the method to be used for integrating statistics across treatments for each 
+#' substrate (phosphorylation site). Available methods are Stouffer, OSP, Fisher, and maxP. 
+#' Default method is OSP.
+#' @param kinase.method the method to be used for integrating statistics of all phosphorylation
+#' sites that belongs to a kinase. Available methods are Stouffer, OSP, Fisher, and maxP. 
+#' Default method is Stouffer. 
+#' @param visualize whether to visualize the plot.
+#' @param ... other visualization parameters to pass on. 
+#' @return a list that contains directional p-values for each substrate and directional enrichment for 
+#' each kinase.
+#' @export
+#' @examples
+#' 
+#' # load the phosphoproteomics dataset
+#' data(HEK)
+#' 
+#' # load the kinase-substrate annoations
+#' data(PhosphoSite)
+#' 
+#' # direction pathway analysis in 2-dimensional space. Implemented as rotating by degree 
+#' # (1) test combined effect of Torin1 and Rapamycin vs insul both on "down-regulation"
+#' # (180 degree to original direction)
+#' kPA <- kinasePA(Tc=HEK, direction=pi, annotation=PhosphoSite.mouse)
+#' kPA$kinase[order(unlist(kPA$kinase[,1])),][1:20,]
+#' # rank substrates on the direciton of interest
+#' sort(kPA$substrate.pvalues)[1:20]
+#' 
+#' # (2) test combined effect of Torin1 and Rapamycin vs insul on "no change and down-regulation"
+#' # (135 degree to the original direction) 
+#' kPA <- kinasePA(Tc=HEK, direction=pi*3/4, annotation=PhosphoSite.mouse)
+#' kPA$kinase[order(unlist(kPA$kinase[,1])),][1:20,]
+#' 
+#' # (3) test combined effect of Torin1 and Rapamycin vs insul on "down-regulation and no change"
+#' # (225 degree to the original direction) 
+#' kPA <- kinasePA(Tc=HEK, direction=pi*5/4, annotation=PhosphoSite.mouse)
+#' kPA$kinase[order(unlist(kPA$kinase[,1])),][1:20,]
+#' 
+kinasePA <- function(Tc, direction, annotation, minSize=5, substrate.method="OSP", kinase.method="Stouffer", visualize=TRUE, ...){
+   dPA <- directPA(Tc, direction, annotation, minSize, gene.method=substrate.method, path.method=kinase.method, visualize, ...)
+   kPA <- list()
+   kPA$substrate.pvalues <- dPA$gene.pvalues 
+   kPA$kinase <- dPA$pathways
+   return(kPA)
+}

source/directPA_1.1/R/pathwayStats.R → R/pathwayStats.R

@@ -28,7 +28,7 @@
 #' 
 #' # Rotate the matrix by contrast 1, -1, -1 (i.e. up-regulation, down-regulation, dow-regulation).
 #' PM.rotated <- rotate3d(PM.zscores, contrast = c(1, -1, -1))
-#'  
+#'
 #' # combine rotated statistics across treatments
 #' gene.pvalues <- apply(PM.rotated, 1, geneStats)
 #'