version 1.0

cran · Nov 22, 2017 · de5d2b8 · de5d2b8
1 parent f6b8ec3
commit de5d2b8
Show file tree

Hide file tree

Showing 23 changed files with 481 additions and 393 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,25 +1,25 @@
 Type: Package
 Title: Multivariate Information Inductive Causation
-Version: 0.1
-Date: 2017-10-09
+Version: 1.0
+Date: 2017-11-11
 Package: miic
 Description: We report an information-theoretic method which learns a large class of causal or non-causal graphical models from purely observational data, while including the effects of unobserved latent variables, commonly found in many datasets. Starting from a complete graph, the method iteratively removes dispensable edges, by uncovering significant information contributions from indirect paths, and assesses edge-specific confidences from randomization of available data. The remaining edges are then oriented based on the signature of causality in observational data. This approach can be applied on a wide range of datasets and provide new biological insights on regulatory networks from single cell expression data, genomic alterations during tumor development and co-evolving residues in protein structures. For more information you can refer to: Verny et al. Plos Comput Biol. (2017) <doi:10.1371/journal.pcbi.1005662>.
 Authors@R: c(person("Nadir", "Sella", role = c("aut","cre"), email = "nadir.sella@curie.fr"),
            person("Louis", "Verny",role = "aut"),
            person("Severine", "Affeldt", role = "aut"),
            person("Hervé", "Isambert", role = c("aut"), email = "Herve.Isambert@curie.fr"))
 Maintainer: Nadir Sella <nadir.sella@curie.fr>
-Imports: MASS,igraph, methods, plotrix, bnlearn, Rcpp, ppcor
+Imports: MASS, igraph, bnlearn, ppcor, stats, Rcpp
 License: GPL (>= 2)
 NeedsCompilation: yes
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 6.0.1
 LinkingTo: Rcpp
-Packaged: 2017-10-09 14:35:54 UTC; nadir
+Packaged: 2017-11-22 16:47:43 UTC; nadir
 Author: Nadir Sella [aut, cre],
   Louis Verny [aut],
   Severine Affeldt [aut],
   Hervé Isambert [aut]
 Repository: CRAN
-Date/Publication: 2017-10-09 15:54:34 UTC
+Date/Publication: 2017-11-22 16:57:01 UTC
diff --git a/MD5 b/MD5
@@ -1,17 +1,17 @@
-a6fdbc256cb2630ced8bbd4c7c218f6b *DESCRIPTION
+025686c756633c295cdadeca670e5189 *DESCRIPTION
 b50205a65608e6b31f7140a07576af20 *NAMESPACE
 8a83329490f3bc92cbb5a3aa459ab69c *R/data.R
-b87cae829aa82bad8fd7a2a25309885b *R/evaluate_efn.R
+15d73959b48006dd12eb39f61f0c6c16 *R/evaluate_efn.R
 3a37815bb41d603c82883dbcbcd99626 *R/gmPlot.lib.R
 710b1febf091e82b3f33f6f7f8b7eae6 *R/gmStatistics.orient.lib.R
 23a27aa69c6e7fdd30d8ead240e655a1 *R/gmStatistics.skeleton.lib.R
 0a6be7d93839500f11bc1991d97c78d1 *R/gmSummary.R
 c19120d1c18e25e86db5767dcfcc8baf *R/gmSummary.lib.R
-a394543dff5a97048508c9181e922912 *R/miic.R
-f338965e7fa825c25f4f266d335dfd55 *R/miic.orient.R
+1f8e49561cf3be399f6078676bf13cd1 *R/miic.R
+f56a076b1cd577c7eca658d3ce3dae5f *R/miic.orient.R
 f6cd4f12a2bed4020d51f6390ee3c7f5 *R/miic.plot.R
-5697c881c566cab959170e4c18c067e3 *R/miic.skeleton.R
-1db813647f8d46f894e6aa7c8640155e *R/miic.utils.R
+9e249cf685f2d2c3e3e2405e76c9b429 *R/miic.skeleton.R
+e67bc298c655e8835df7dd757287b284 *R/miic.utils.R
 95897d5733d3844b947af2488d8d2ed4 *R/shared.utils.lib.R
 6e20ed7b7261c755396ccd6a2fb5b13e *R/write.cytoscape.R
 56f5fb69772a96249d878df39f9bd695 *R/write.style.R
@@ -24,37 +24,36 @@ f27116874bf4e6c19e8d178eb33cac92 *data/ohno_stateOrder.rda
 e9304c35e64bb0d689e9337ee62f3885 *man/cosmicCancer.Rd
 ee82115615692931b2533d7f9b01c67c *man/cosmicCancer_stateOrder.Rd
 6936b8176557a4dfc94e9fd1e3df4551 *man/hematoData.Rd
-7f43f6ac4c416548d6329f5b251470a7 *man/miic.Rd
+b1af86a80fc0dbb3580bc1f19f821e39 *man/miic.Rd
 641f80420a52817ce17ece6671e66a15 *man/miic.evaluate.effn.Rd
 843ed9f299b77a3102519ab320a7d22c *man/miic.plot.Rd
 dd2b3833e4079d797dfb278008da3639 *man/miic.write.network.cytoscape.Rd
 d0a7e312ede8df92a1240855c6853240 *man/miic.write.style.cytoscape.Rd
 6edc3db98f08bb12a0d416417e6d344e *man/ohno.Rd
 7817c89c2c9177e6f5cbdeab9742ee56 *man/ohno_stateOrder.Rd
 9dfda68c998b80c1539be1ea141f1c6f *src/Makevars
-b0e0a870b963f94fb3d0e33504c91b65 *src/computeEnsInformation.cpp
+3455c4ca65015e052e24fc7b818722d9 *src/computeEnsInformation.cpp
 fd6cb900f21d3eade279ebcd8e7dbba8 *src/computeEnsInformation.h
-a6d7806d77368b2f2765c98c61c856e0 *src/computeInfo.cpp
+e274e55bfec8e6738152044be2faa9f1 *src/computeInfo.cpp
 b92099fde178adad20c48b3d1b826928 *src/computeInfo.h
-9d3d5c5dd59800f6cc38a46581a875d1 *src/computeInfo_interface.cpp
+94ca49ed6a869bef8509a56a32657564 *src/computeInfo_interface.cpp
 5cce2bdc080024287e6f2dcadd9f9fdc *src/computeInfo_interface.h
-e0825921467980d53160f447963cee2d *src/confidenceCut.cpp
+3ef290b00de4a93037ee99824300ab93 *src/confidenceCut.cpp
 254a0cc0d3a35cfc7d8751ec67021a78 *src/confidenceCut.h
-ea8c2e36af98501e6c6a4d554efd2fe1 *src/crosscorrelation.cpp
+610c9a669fbc5b06ebcc50833517bc4b *src/crosscorrelation.cpp
 6f386fbfa1601b35da8b8f9921c8a600 *src/memory.h
-0d7b6ebdba5c16c9354e95da5b7466f1 *src/miic_init.c
-53af7dedc80e724265ec834e5e6dff96 *src/orientationProbability.cpp
+548accdd415f9466471ac55f36cc4b37 *src/miic_init.c
+9c832029fb04eddcc0e7d8a37edc14c9 *src/orientationProbability.cpp
 df81b15d577e2284c88218618f86c029 *src/orientationProbability.h
-15774d652f8b40cbda603e776f757c16 *src/probaOrientation.cpp
+a3b3a987a3b248022a60647b547ebafb *src/probaOrientation.cpp
 bf4e298c2e147cfa6e17d239c9d3a381 *src/probaOrientation.h
 2a3d543c6b5132ac56366927a7b38f68 *src/probaOrientation_interface.cpp
 0631912416f071d2dff6da62f0c1fc77 *src/probaOrientation_interface.h
-3c958989b52b7268b3b4e137251ac817 *src/skeleton.cpp
-f10078f15d83c1321742cfad6aa75f6e *src/skeleton.h
-b174ec1055fb1cdcc13742af1041f9e3 *src/skeletonInitialization.cpp
+cd8d370b069ccf5952f531bf3c1f909d *src/skeleton.cpp
+5b0686d955df95fadb5ce553d5685c4d *src/skeletonInitialization.cpp
 aad43dd4310bfbc15eca725c08f4db17 *src/skeletonInitialization.h
-99e7531c85e5e01872e2c86f91b277bb *src/skeletonIteration.cpp
+dbd4e526309abd82b43065a0bdbe6c10 *src/skeletonIteration.cpp
 28a37bc6efab65d5b2d965a9b897f80f *src/skeletonIteration.h
-73800446d485674606d7151c9f78fcd4 *src/structure.h
-ffb184f3990522d7948fd3632d3c29e1 *src/utilities.cpp
-21e78659556c299a5e80d952126585f9 *src/utilities.h
+7efacf8d0693eeb4c50b05b244127575 *src/structure.h
+9232cefaed632a26a0d5e28d83813446 *src/utilities.cpp
+23a6dbcf0782dc004c4280e8160b87f7 *src/utilities.h
diff --git a/R/evaluate_efn.R b/R/evaluate_efn.R
@@ -13,29 +13,32 @@
 
 miic.evaluate.effn <- function(inputData = NULL, plot=T)
 {
+  result = list()
   #### Check the input arguments
   if( is.null( inputData ) )
   { stop("The input data file is required") }
   inData <- c(colnames(inputData), as.vector(as.character(t(as.matrix(inputData)))))
-  if (requireNamespace("Rcpp", quietly = TRUE)) {
-    res <- .Call('evaluateEffn', inData, ncol(inputData), nrow(inputData),PACKAGE = "miic")
+  if (base::requireNamespace("Rcpp", quietly = TRUE)) {
+    result <- .Call('evaluateEffn', inData, ncol(inputData), nrow(inputData),PACKAGE = "miic")
   }
-  if(length(which(res$correlation > 0)) == length(res$correlation)){
+  if(length(which(result$correlation > 0)) == length(result$correlation)){
 
-    fit1 <- MASS::fitdistr(res$correlation, "exponential")
-    pval = stats::ks.test(res$correlation, "pexp", fit1$estimate)$p.value
+    fit1 <- MASS::fitdistr(result$correlation, "exponential")
+    pval = stats::ks.test(result$correlation, "pexp", fit1$estimate)$p.value
     if(pval < 0.05){
-      res$exponential_decay= FALSE
+      result$exponential_decay= FALSE
     } else {
-      res$exponential_decay= TRUE
+      result$exponential_decay= TRUE
     }
 
     if(plot){
-      graphics::plot(res$correlation, type="l", log="y", ylab="Autocorrelation with lag", xlab="n")
-      graphics::title("Autocorrelation between n distant samples")
+      graphics::par(mar=rep(1.5, 4), oma=c(3,3,3,1), las=1)
+
+      graphics::plot(0:(length(result$correlation)-1),result$correlation, type="l", log="y")
+      graphics::title("Autocorrelation between n distant samples",ylab="Autocorrelation with lag", xlab="n")
     }
   }else {
-    res$exponential_decay= FALSE
+    result$exponential_decay= FALSE
   }
-  res
+   result
 }
diff --git a/R/miic.R b/R/miic.R
@@ -187,11 +187,13 @@
 #'
 #' # plot graph
 #' miic.plot(miic.res)
+#'\dontrun{
 #'
 #' # write graph to graphml format. Note that to correctly visualize
 #' # the network we created the miic style for Cytoscape (http://www.cytoscape.org/).
-#' miic.write.network.cytoscape(g = miic.res, file=paste(tempdir(),"/temp",sep=""))
-#'\dontrun{
+#'
+#' miic.write.network.cytoscape(g = miic.res, file.path(tempdir(),"/temp"))
+#'
 #' # EXAMPLE CANCER
 #' data(cosmicCancer)
 #' data(cosmicCancer_stateOrder)
@@ -204,7 +206,7 @@
 #'
 #' # write graph to graphml format. Note that to correctly visualize
 #' # the network we created the miic style for Cytoscape (http://www.cytoscape.org/).
-#' miic.write.network.cytoscape(g = miic.res, file = paste(tempdir(),"/temp",sep=""))
+#' miic.write.network.cytoscape(g = miic.res, file = file.path(tempdir(),"/temp"))
 #'
 #' # EXAMPLE OHNOLOGS
 #' data(ohno)
@@ -218,7 +220,7 @@
 #'
 #' # write graph to graphml format. Note that to correctly visualize
 #' # the network we created the miic style for Cytoscape (http://www.cytoscape.org/).
-#' miic.write.network.cytoscape(g = miic.res, file = paste(tempdir(),"/temp",sep=""))
+#' miic.write.network.cytoscape(g = miic.res, file = file.path(tempdir(),"/temp"))
 #'}
 
 miic <- function(inputData = NULL, categoryOrder= NULL, trueEdges = NULL, blackBox = NULL, nThreads=1,
@@ -242,6 +244,17 @@ miic <- function(inputData = NULL, categoryOrder= NULL, trueEdges = NULL, blackB
   if( !is.data.frame(inputData))
   { stop("The input data is not a dataframe") }
 
+  effnAnalysis = miic.evaluate.effn(inputData, plot=F)
+  if(effnAnalysis$neff < 0.5 * nrow(inputData) ){
+    if(effnAnalysis$exponential_decay){
+        print(paste("Warning! Your samples in the datasets seem to be correlated! We suggest to re run the method specifying",
+          effnAnalysis$neff, " in the neff parameter. See the autocorrelation plot for more details."))
+      } else {
+        print("Warning! Your samples in the datasets seem to be correlated but the correlation decay is not exponential.
+          Are your samples correlated in some way? See the autocorrelation plot for more details.")
+      }
+  }
+
   # if(!is.null(trueEdges)){
   #   if(length(which(!c(as.vector(trueEdges[,1]), as.vector(trueEdges[,2])) %in% colnames(inputData))) > 0){
   #     stop("True edges file does not correspond to the input data matrix. Please check node names.")
@@ -288,10 +301,13 @@ miic <- function(inputData = NULL, categoryOrder= NULL, trueEdges = NULL, blackB
     if(skeleton){
       if(verbose)
        cat("\t# -> START skeleton...\n")
-      res <- miic.skeleton(inputData = inputData, stateOrder= categoryOrder, nThreads= nThreads, cplx = cplx, latent = latent, effN = neff, blackBox = blackBox, verbose= verbose)
+      res <- miic.skeleton(inputData = inputData, stateOrder= categoryOrder, nThreads= nThreads, cplx = cplx,
+                          latent = latent, effN = neff, blackBox = blackBox,
+                          confidenceShuffle = confidenceShuffle, confidenceThreshold= confidenceThreshold, verbose= verbose)
       # print(res)
       edges = res$edges
-
+      confData = res$confData
+      time = res$time
       if(verbose)
         cat("\t# -> END skeleton...\n\t# --------\n")
     }
@@ -302,16 +318,21 @@ miic <- function(inputData = NULL, categoryOrder= NULL, trueEdges = NULL, blackB
     }
 
 
+    timeOrt=0
     if(orientation){
 
       if(verbose)
         cat("\tSTART orientation...")
-
+      ptm <- proc.time()
       res = miic.orient(inputData= inputData, method = orientationMethod, stateOrder = categoryOrder,
                           edges = edges, effN = neff,  cplx = cplx,  latent = latent, propagation = propagation,
-                          verbose = FALSE, confidenceShuffle = confidenceShuffle, confidenceThreshold= confidenceThreshold)
-      res$edges <- edges
+                          verbose = FALSE)
 
+      timeOrt=(proc.time() - ptm)["elapsed"]
+      timeInitIterOrt = time["initIter"]+timeOrt
+
+      res$edges <- edges
+      res$confData <- confData
       if(verbose)
         cat("\tEND orientation...")
     }
@@ -341,9 +362,18 @@ miic <- function(inputData = NULL, categoryOrder= NULL, trueEdges = NULL, blackB
 
     # Call the function
 
+    ptm <- proc.time()
     resGmSummary <- gmSummary(inputData = inputData, edges = res$edges,
                               adjMatrix= res$adjMatrix, trueEdges = trueEdges, stateOrder = categoryOrder, verbose = verbose)
 
+    timeInitIterOrt = timeOrt + time[4]
+    timeSum=(proc.time() - ptm)["elapsed"]
+    timeTotal = timeInitIterOrt+timeSum
+    timeVec = c(time, timeOrt, timeInitIterOrt, timeSum, timeTotal)
+    timeVec[which(timeVec==0)]=NA
+    res$time = stats::setNames(timeVec,
+                c("initailization", "iteration", "confidenceCut", "skeleton", "orientation", "skeleton+Orientation", "summary", "total"))
+
     res$all.edges.summary <- resGmSummary$all.edges.summary
     res$retained.edges.summary <- resGmSummary$retained.edges.summary
     res$statistics = resGmSummary$statistics

diff --git a/R/miic.orient.R b/R/miic.orient.R
@@ -1,9 +1,8 @@
 
 miic.orient <- function(inputData = NULL, method = c("probabilistic", "logic"), stateOrder = NULL, edges = NULL,
                              effN = -1, cplx = c("nml", "mdl"), eta = 1,
-                             latent = FALSE, propagation = TRUE, hvs = FALSE, continuous = 0, verbose = FALSE, 
-                             confidenceShuffle = 0, confidenceThreshold= 0
-){
+                             latent = FALSE, propagation = TRUE, hvs = FALSE, continuous = 0, verbose = FALSE)
+{
   isK23 = TRUE
   isDegeneracy = FALSE
 
@@ -50,7 +49,7 @@ miic.orient <- function(inputData = NULL, method = c("probabilistic", "logic"),
   if(method=="logic" ){
 
     res <- .Call('orientation',inData, numNodes, edges, effN, cplx, eta, latent,
-                isK23, isDegeneracy, propagation, continuous, verbose,PACKAGE = "miic")
+                isK23, isDegeneracy, propagation, continuous, verbose, PACKAGE = "miic")
 
     # create the data frame of the structures before orientation
     tmp <- unlist(res$tableOfOrientationsBeforePropagation)[1:length(res$tableOfOrientationsBeforePropagation[[1]])]
@@ -84,7 +83,7 @@ miic.orient <- function(inputData = NULL, method = c("probabilistic", "logic"),
   } else if( method=="probabilistic" ){
 
     res <- .Call('orientationProbability',inData, numNodes, edges, effN, cplx, eta, hvs, latent,
-                 isK23, isDegeneracy, propagation, continuous, verbose, confidenceShuffle, confidenceThreshold)
+                 isK23, isDegeneracy, propagation, continuous, verbose)
 
     #create the data frame of the structures after orientation
     df = res$orientations.prob
@@ -117,21 +116,5 @@ miic.orient <- function(inputData = NULL, method = c("probabilistic", "logic"),
   # update the returned adj matrix
   res$adjMatrix <- df
 
-
-
-  if(confidenceShuffle > 0){
-    # create the data frame for the confidence file
-    tmp <- unlist(res$confData)[1:length(res$confData[[1]])]
-    res1 <- unlist(res$confData)[(length(res$confData[[1]])+1):length(unlist(res$confData))]
-    df <- data.frame(matrix(res1, nrow=length(res$confData)-1, byrow=T),stringsAsFactors=FALSE)
-    colnames(df) <-tmp
-    df[,"confidence_ratio"] = as.numeric(df[,"confidence_ratio"])
-
-    df <-df[order(df[,"confidence_ratio"]),]
-
-
-    res$confData <- df
-
-  }
   res
 }