From 4077cc307664a9541a3d37691d9b2513fbd1237f Mon Sep 17 00:00:00 2001
From: Ghanbari <mghanba@ohler34.mdc-berlin.net>
Date: Thu, 5 Oct 2017 17:42:28 +0200
Subject: [PATCH] update R package

---
 NAMESPACE           |  2 +-
 R/fisher_test.R     | 51 ---------------------------------------------
 R/kmeans_links.R    | 45 +++++++++++++++++++++++++++++++++++++++
 README.md           | 32 ++++++++++++++++++++--------
 man/fisher_test.Rd  | 41 ------------------------------------
 man/kmeans_links.Rd | 35 +++++++++++++++++++++++++++++++
 6 files changed, 104 insertions(+), 102 deletions(-)
 delete mode 100644 R/fisher_test.R
 create mode 100644 R/kmeans_links.R
 delete mode 100644 man/fisher_test.Rd
 create mode 100644 man/kmeans_links.Rd

diff --git a/NAMESPACE b/NAMESPACE
index f59fae7..5e740e3 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,8 +2,8 @@
 
 export(Dcenter)
 export(dpm)
-export(fisher_test)
 export(get_link)
+export(kmeans_links)
 export(nldata)
 export(rdata)
 export(reg.dpm)
diff --git a/R/fisher_test.R b/R/fisher_test.R
deleted file mode 100644
index d5c72fc..0000000
--- a/R/fisher_test.R
+++ /dev/null
@@ -1,51 +0,0 @@
-#' Function to compute fisher z-transformation pvalues for the links
-#'
-#' @param score_matrix score matrix obtained from DPM or reg.DPM.
-#' @param alpha a significance level
-#' @param nsamp is the number of samples (number of rows of input matrix for (reg-)DPM ). Note that we need nsamp such that: nsamp*(nsamp+1)> 2*(nvar-5)
-#' @param nvar is the number of variables (number of columns of input matrix for (reg-)DPM )
-#'
-#' @return a list containing
-#' - the list of significant links as a data frame. Each line of the data frame corresponds to a link.
-#' The first and second columns are the corresponding nodes of the link, the third column is the score of the link and the forth column is the corresponding pvalue.
-#' - the threshold that all the links with scores more than it are significant at defined significance level.
-#'
-#' @examples
-#' ## Load DREAM4 data (insilico_size10_1_knockouts)
-#' data(dream4)
-#'
-#'## Run regularized DPM
-#' res <- reg.dpm(data)
-#'
-#' ## Get the list of significant links at significance level alpha=0.05
-#' fisher_test(res,alpha=0.05,10,10)$LinkList
-#'
-#' ##Get the threshold at significance level alpha=0.05
-#' fisher_test(score_matrix=res,alpha=0.05,nsamp=10,nvar=10)$threshold
-#'
-#'
-#' @export
-
-###########################################################
-### Function to compute fisher z transformation pvalues ###
-###########################################################
-
-fisher_test <- function(score_matrix, alpha, nsamp, nvar){
-  nsamp <- (nsamp*(nsamp+1))/2
-  if (nsamp-(nvar-2)-3 <= 0 ) stop("the test requires nsamp such that: nsamp*(nsamp+1)> 2*(nvar-5)")
-  diag(score_matrix) <- NA
-  score_matrix[upper.tri(score_matrix)] <- NA
-  linkList <- reshape2::melt(score_matrix, na.rm=TRUE)
-  colnames(linkList) <- c("Node1", "Node2", "Score")
-  p <- linkList["Score"]
-  z <- sqrt(nsamp-(nvar-2)-3) * abs(atan(p))
-  pval <- 2 * pnorm(-as.matrix(z))
-  colnames(pval) <- "Pvalue"
-  linkList <- cbind(linkList,pval)
-  linkList <- linkList[linkList$Pvalue<alpha,]
-  temp <- exp( 2 * pval_tr / (sqrt(nsamp-(nvar-2)-3)) )
-  threshold <- (temp -1)/(temp + 1)
-  return(list(LinkList=linkList,threshold=threshold))
-}
-
-
diff --git a/R/kmeans_links.R b/R/kmeans_links.R
new file mode 100644
index 0000000..c7ba6e2
--- /dev/null
+++ b/R/kmeans_links.R
@@ -0,0 +1,45 @@
+#' Function to apply kmeans clustering on the score of links
+#'
+#' @param score_matrix score matrix obtained from DPM or reg.DPM.
+#'
+#' @return a list containing
+#' - the list of the links as a data frame. Each line of the data frame corresponds to a link.
+#' The first and second columns are the corresponding nodes of the link and the third column is the score of the link.
+#' - the threshold that all the links with scores more than it are considered as an edge.
+#'
+#' @examples
+#' ## Load DREAM4 data (insilico_size10_1_knockouts)
+#' data(dream4)
+#'
+#'## Run regularized DPM
+#' res <- reg.dpm(data)
+#'
+#' ## Get the list of the links based on kmeans clustering
+#' linklist <- kmeans_links(res)$LinkList
+#'
+#' ##Get the threshold at significance level alpha=0.05
+#' kmeans_links(res)$threshold
+#'
+#'
+#' @export
+
+###########################################################
+### Function to compute fisher z transformation pvalues ###
+###########################################################
+kmeans_links <- function(score_matrix){
+  score_matrix <- abs(score_matrix)
+  diag(score_matrix) <- NA
+  score_matrix[upper.tri(score_matrix)] <- NA
+  linkList <- reshape2::melt(score_matrix, na.rm=TRUE)
+  colnames(linkList) <- c("Node1", "Node2", "Score")
+
+  # do kmeans on scores
+  res = kmeans(linkList$Score, centers=2)
+  linkList <- linkList[res$cluster==which.max(res$centers),]
+
+  # compute threshold
+  threshold = mean(a$centers)
+
+return(list(LinkList=linkList,threshold=threshold))
+}
+
diff --git a/README.md b/README.md
index f5c287a..0b81a77 100644
--- a/README.md
+++ b/README.md
@@ -15,6 +15,7 @@ Then install DPM using the `install_github` function in the
 ```r
 library(devtools)
 devtools::install_github("ghanbari/DPM",host="github.molgen.mpg.de/api/v3")
+library(DPM)
 ```
 
 ### Format of input data
@@ -29,6 +30,7 @@ head(DataMat)
 In the context of gene networks variables correspond to genes. The input is gene expression data as a matrix in which columns correspond to genes and rows correspond to samples. As an example, we have one data from DREAM4 challenge (insilico_size10_1_knockouts):
 
 ```r
+library(DPM)
 data(dream4)
 head(dream4)
 ```
@@ -90,22 +92,34 @@ ggsave(“ScoreDistribution.png”, plot = last_plot(), device = “png”)
 When you decide for the threshold you can visualize the resulted network (see section [Visualization of the network](#visualization-of-the-network) ) and get all the links with a score higher than the threshold (see subsection [Get only the links with a score higher than a threshold](#get-only-the-links-with-a-score-higher-than-a-threshold) )
 
 
-However, we provide a threshold selection method based on the [fisher z-transformation](https://en.wikipedia.org/wiki/Fisher_transformation). You can compute the threshold for the DREAM4 data with the following command where:
+However, we provide a threshold selection method based on k-means clustering. You can compute the threshold for the DREAM4 data with the following command where:
 
 * `score_matrix` is the score matrix obtained from DPM or reg.DPM.
-* `alpha` is a significance value
-* `nsamp` is the number of samples (number of rows of input matrix for (reg-)DPM )
-* `nvar` is the number of variables (number of columns of input matrix for (reg-)DPM )
 
 ```r
-tr <- fisher_test(dpmres, alpha=0.05, nsamp=nrow(dream4), nnodes=ncol(dream4))$threshold
-tr <- fisher_test(regdpmres, alpha=0.05, nsamp=nrow(dream4), nnodes=ncol(dream4))$threshold
+dpmtr <- kmeans_links(dpmres)$threshold
+regdpmtr <- kmeans_links(regdpmres)$threshold
 ```
 
-You can also get all the significant links based on fisher z-transformation at significance level alpha:
+You can visualize the position of the threshold on the distribution of the scores:
+
+```r
+library("ggplot2")
+dpmdat <- data.frame(score = abs(dpmres[lower.tri(dpmres)]),name ='DPM')
+regdpmdat <-data.frame(score = abs( regdpmres[lower.tri(regdpmres)]), name='reg-DPM')
+dat <- rbind(dpmdat,regdpmdat)
+
+ggplot(dat, aes(x=score,color=name))+ theme_bw()  + geom_density(size=1.5)+scale_color_manual(values=c("DPM"="SteelBlue3","reg-DPM"="dodgerblue4"))+ geom_vline(xintercept=dpmtr, size=1, color="SteelBlue3") + geom_vline(xintercept=regdpmtr, size=1,color="dodgerblue4")
+
+ggsave(“ScoreDistribution.png”, plot = last_plot(), device = “png”)
+```
+
+
+You can also get all the links considered as true links based on k-means clustering:
+
 ```r
-dpm_sig_links <- fisher_test(dpmres, alpha=0.05, nsamp=nrow(dream4), nnodes=ncol(dream4))$Listlink
-regdpm_sig_links <- fisher_test(regdpmres, alpha=0.05, nsamp=nrow(dream4), nnodes=ncol(dream4))$Listlink
+dpm_links <- kmeans_links(dpmres)$Listlink
+regdpm_links <- kmeans_links(regdpmres)$Listlink
 ```
 
 As an alternative, you can select the top-ranked links based on (reg.)DPM score which are most likely to be true links.
diff --git a/man/fisher_test.Rd b/man/fisher_test.Rd
deleted file mode 100644
index 6d4a87f..0000000
--- a/man/fisher_test.Rd
+++ /dev/null
@@ -1,41 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/fisher_test.R
-\name{fisher_test}
-\alias{fisher_test}
-\title{Function to compute fisher z-transformation pvalues for the links}
-\usage{
-fisher_test(score_matrix, alpha, nsamp, nvar)
-}
-\arguments{
-\item{score_matrix}{score matrix obtained from DPM or reg.DPM.}
-
-\item{alpha}{a significance level}
-
-\item{nsamp}{is the number of samples (number of rows of input matrix for (reg-)DPM ). Note that we need nsamp such that: nsamp*(nsamp+1)> 2*(nvar-5)}
-
-\item{nvar}{is the number of variables (number of columns of input matrix for (reg-)DPM )}
-}
-\value{
-a list containing
-- the list of significant links as a data frame. Each line of the data frame corresponds to a link.
-The first and second columns are the corresponding nodes of the link, the third column is the score of the link and the forth column is the corresponding pvalue.
-- the threshold that all the links with scores more than it are significant at defined significance level.
-}
-\description{
-Function to compute fisher z-transformation pvalues for the links
-}
-\examples{
-## Load DREAM4 data (insilico_size10_1_knockouts)
-data(dream4)
-
-## Run regularized DPM
-res <- reg.dpm(data)
-
-## Get the list of significant links at significance level alpha=0.05
-fisher_test(res,alpha=0.05,10,10)$LinkList
-
-##Get the threshold at significance level alpha=0.05
-fisher_test(score_matrix=res,alpha=0.05,nsamp=10,nvar=10)$threshold
-
-
-}
diff --git a/man/kmeans_links.Rd b/man/kmeans_links.Rd
new file mode 100644
index 0000000..06e8e92
--- /dev/null
+++ b/man/kmeans_links.Rd
@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/kmeans_links.R
+\name{kmeans_links}
+\alias{kmeans_links}
+\title{Function to apply kmeans clustering on the score of links}
+\usage{
+kmeans_links(score_matrix)
+}
+\arguments{
+\item{score_matrix}{score matrix obtained from DPM or reg.DPM.}
+}
+\value{
+a list containing
+- the list of the links as a data frame. Each line of the data frame corresponds to a link.
+The first and second columns are the corresponding nodes of the link and the third column is the score of the link.
+- the threshold that all the links with scores more than it are considered as an edge.
+}
+\description{
+Function to apply kmeans clustering on the score of links
+}
+\examples{
+## Load DREAM4 data (insilico_size10_1_knockouts)
+data(dream4)
+
+## Run regularized DPM
+res <- reg.dpm(data)
+
+## Get the list of the links based on kmeans clustering
+linklist <- kmeans_links(res)$LinkList
+
+##Get the threshold at significance level alpha=0.05
+kmeans_links(res)$threshold
+
+
+}