R/omics_network.R

In Patterns: Deciphering Biological Networks with Patterned Heterogeneous Measurements

#' @rdname show-methods
setMethod(f= "show", signature = "omics_network", definition = function(object) {
cat(
  paste(
    "This is a S4 class with : \n - (@omics_network) a matrix of dimension ",
    dim(object@omics_network)[1],
    "*",
    dim(object@omics_network)[2],
    " .... [the omics_network] \n - (@name) a vector of length ",
    length(object@name),
    " .... [omics names] \n",
    "- (@F) a array of dimension ",
    dim(object@F)[1],
    "*",
    dim(object@F)[2],
    "*",
    dim(object@F)[3],
    " .... [F matrices] \n",
    "- (@convF) a matrix of dimension ",
    dim(object@convF)[1],
    "*",
    dim(object@convF)[2],
    " .... [convergence (L1 norm) of array F] \n",
    "- (@convO)a vector of length ",
    length(object@convO),
    " .... [convergence (L1 norm) of matrix Omega]\n",
    "- (@time_pt) an vector of length",
    length(object@time_pt),
    "  .... [time points]"
  )
)
invisible(object)
})

#' Analysing the network
#' 
#' Calculates some indicators for each node in the network.
#' 
#' 
#' @aliases analyze_network analyze_network-methods
#' analyze_network,omics_network-method
#' @param Omega a omics_network object
#' @param nv the level of cutoff at which the analysis should be done
#' @param label_v (optionnal) the name of the genes
#' @return A matrix containing, for each node, its betweenness,its degree, its
#' output, its closeness.
#' @author Bertrand Frederic, Myriam Maumy-Bertrand.
#' @keywords methods
#' @examples
#' 
#' data(network)
#' analyze_network(network,nv=0)
#' 
setMethod("analyze_network", "omics_network", function(Omega, nv, label_v = NULL) {
  require(tnet)
  if (is.null(label_v)) {
    label_v <- 1:dim(Omega@omics_network)[1]
  }
  O <- Omega@omics_network
  Omega <- Omega@omics_network
  O[abs(O) <= nv] <- 0
  G <- graph.adjacency(O, weighted = TRUE)
  get.edgelist(G) -> Q
  weight <- rep(0, dim(Q)[1])
  for (i in 1:dim(Q)[1]) {
    weight[i] <- Omega[Q[i, 1], Q[i, 2]]
  }
  R <- cbind(Q, abs(weight))
  G <- as.tnet(R,type="weighted one-mode tnet")
  C <-
    data.frame(label_v,
               betweenness_w(G)[, 2],
               degree_w(G)[, 2:3],
               closeness_w(G, gconly = FALSE)[, 2])
  colnames(C) <- c("node", "betweenness", "degree", "output", "closeness")
  return(C)
})


#' See the evolution of the network with change of cutoff
#' 
#' See the evolution of the network with change of cutoff
#' 
#' Several types of outputs are available using the type.ani option. \itemize{
#' \item html \item latex (requires latex) \item swf (requires swftools)
#' \item video (requires ffmpeg) \item gif \item manual_gif }
#' 
#' @aliases evolution evolution-methods evolution,omics_network-method
#' @param net a omics_network object
#' @param list_nv a vector of cutoff at which the network should be shown
#' @param gr a vector giving the group of each genee. Defaults to NULL 
#' @param color.vertex a vector giving the color of each nodee. Defaults to NULL 
#' @param color.edge a vector giving the color of each edge. Defaults to NULL 
#' @param fix logical, should the position of the node in the network be calculated once at the beginning ? Defaults to TRUE. 
#' @param size vector giving the size of the plot. Defaults to c(2000,1000) 
#' @param label_v vector giving the labels of each vertex. Defaults to 1:dim(net@omics_network)[1] 
#' @param legend string giving the position of the legend. Defaults to "topleft" 
#' @param legend.position string giving the position of the legend. Defaults to "topleft" 
#' @param frame.color string giving the color of the frame of the plot. Defaults to "black" 
#' @param label.hub label hubs. Defaults to FALSE 
#' @param outdir Directory to save the animation. No default value since it must be specified by the user. 
#' @param type.ani Type of animation. Defaults to "html" 
#' @return A HTML page with the evolution of the network.
#' @author Bertrand Frederic, Myriam Maumy-Bertrand.
#' @keywords methods
#' @examples
#' 
#' \donttest{
#' data(network)
#' sequence<-seq(0,0.2,length.out=20)
#' 
#' #Change the destdir to have the animation created where you want. 
#' destdir = tempdir()
#' 
#' #Example of use of the evolution method with an html output.
#' evolution(network,sequence,type.ani = "html", outdir=destdir)
#' }
#' 
#' \dontrun{
#' #Example of use of the evolution method with an animated gif output.
#' evolution(network,sequence,type.ani = "gif", outdir=destdir)
#' }
#' 
setMethod("evolution", "omics_network", 
          function(net
                   ,
                   list_nv
                   ,
                   gr = NULL
                   ,
                   color.vertex = NULL
                   ,
                   color.edge = NULL
                   ,
                   fix = TRUE
                   ,
                   size = c(2000, 1000)
                   ,
                   label_v = 1:dim(net@omics_network)[1]
                   ,
                   legend.position = "topleft"
                   ,
                   frame.color = "black"
                   ,
                   label.hub = FALSE
                   ,
                   outdir
                   ,
                   type.ani = "html"
                   #,edge.arrow.size=0.6*(1+size.ed)
                   #,edge.thickness=1
                   )
                   {
                    if(missing(outdir)){stop("An output directory is required")}
                    #We save the user working directory.
                    orig_working_directory <- getwd()
                    
                    #We restore the user working directory on function exit.
                    on.exit(setwd(orig_working_directory))
                    
                    #Set the working directory to the outdir directory
                    setwd(outdir)
            
                      Omega <- net
                     
                     if (length(list_nv) == 1) {
                       list_nv = seq(0, max(net@omics_network) - 0.05, length.out = list_nv)
                     }
                     
                     if (type.ani == "html") {
                       require(animation)
                       animation::ani.options(ani.height = size[2],
                                              ani.width = size[1])
                       animation::saveHTML({
                         POS <- position(net, nv = list_nv[1])
                         
                         for (i in list_nv) {
                           if (fix == TRUE) {
                             plot(
                               Omega
                               ,
                               nv = i
                               ,
                               gr = gr
                               ,
                               ini = POS
                               ,
                               color.vertex = color.vertex
                               ,
                               color.edge = color.edge
                               ,
                               label_v = label_v
                               ,
                               legend.position = legend.position
                               ,
                               frame.color = frame.color
                               ,
                               label.hub = label.hub
                               #,edge.arrow.size=edge.arrow.size
                               #,edge.thickness=edge.thickness
                             )
                           }
                           else{
                             plot(
                               Omega
                               ,
                               nv = i
                               ,
                               gr = gr
                               ,
                               color.vertex = color.vertex
                               ,
                               color.edge = color.edge
                               ,
                               label_v = label_v
                               ,
                               legend.position = legend.position
                               ,
                               frame.color = frame.color
                               ,
                               label.hub = label.hub
                               #,edge.arrow.size=edge.arrow.size
                               #,edge.thickness=edge.thickness
                             )
                             
                           }
                           text(-1, 1, round(i, 3))
                         }
                       })
                     }
                     if (type.ani == "latex") {
                       require(animation)
                       animation::ani.options(ani.height = size[2],
                                              ani.width = size[1])#, qpdf = '/opt/local/bin/qpdf'
                       animation::saveLatex({
                         POS <- position(net, nv = list_nv[1])
                         
                         for (i in list_nv) {
                           if (fix == TRUE) {
                             plot(
                               Omega
                               ,
                               nv = i
                               ,
                               gr = gr
                               ,
                               ini = POS
                               ,
                               color.vertex = color.vertex
                               ,
                               color.edge = color.edge
                               ,
                               label_v = label_v
                               ,
                               legend.position = legend.position
                               ,
                               frame.color = frame.color
                               ,
                               label.hub = label.hub
                               #,edge.arrow.size=edge.arrow.size
                               #,edge.thickness=edge.thickness
                             )
                           }
                           else{
                             plot(
                               Omega
                               ,
                               nv = i
                               ,
                               gr = gr
                               ,
                               color.vertex = color.vertex
                               ,
                               color.edge = color.edge
                               ,
                               label_v = label_v
                               ,
                               legend.position = legend.position
                               ,
                               frame.color = frame.color
                               ,
                               label.hub = label.hub
                               #,edge.arrow.size=edge.arrow.size
                               #,edge.thickness=edge.thickness
                             )
                             
                           }
                           text(-1, 1, round(i, 3))
                         }
                       })
                     }
                     if (type.ani == "swf") {
                       require(animation)
                       animation::ani.options(ani.height = size[2],
                                              ani.width = size[1])#, qpdf = '/opt/local/bin/qpdf'
                       animation::saveSWF({
                         POS <- position(net, nv = list_nv[1])
                         
                         for (i in list_nv) {
                           if (fix == TRUE) {
                             plot(
                               Omega
                               ,
                               nv = i
                               ,
                               gr = gr
                               ,
                               ini = POS
                               ,
                               color.vertex = color.vertex
                               ,
                               color.edge = color.edge
                               ,
                               label_v = label_v
                               ,
                               legend.position = legend.position
                               ,
                               frame.color = frame.color
                               ,
                               label.hub = label.hub
                               #,edge.arrow.size=edge.arrow.size
                               #,edge.thickness=edge.thickness
                             )
                           }
                           else{
                             plot(
                               Omega
                               ,
                               nv = i
                               ,
                               gr = gr
                               ,
                               color.vertex = color.vertex
                               ,
                               color.edge = color.edge
                               ,
                               label_v = label_v
                               ,
                               legend.position = legend.position
                               ,
                               frame.color = frame.color
                               ,
                               label.hub = label.hub
                               #,edge.arrow.size=edge.arrow.size
                               #,edge.thickness=edge.thickness
                             )
                             
                           }
                           text(-1, 1, round(i, 3))
                         }
                       })
                     }
                     if (type.ani == "video") {
                       require(animation)
                       animation::ani.options(ani.height = size[2],
                                              ani.width = size[1])#, qpdf = '/opt/local/bin/qpdf'
                       animation::saveVideo({
                         POS <- position(net, nv = list_nv[1])
                         
                         for (i in list_nv) {
                           if (fix == TRUE) {
                             plot(
                               Omega
                               ,
                               nv = i
                               ,
                               gr = gr
                               ,
                               ini = POS
                               ,
                               color.vertex = color.vertex
                               ,
                               color.edge = color.edge
                               ,
                               label_v = label_v
                               ,
                               legend.position = legend.position
                               ,
                               frame.color = frame.color
                               ,
                               label.hub = label.hub
                               #,edge.arrow.size=edge.arrow.size
                               #,edge.thickness=edge.thickness
                             )
                           }
                           else{
                             plot(
                               Omega
                               ,
                               nv = i
                               ,
                               gr = gr
                               ,
                               color.vertex = color.vertex
                               ,
                               color.edge = color.edge
                               ,
                               label_v = label_v
                               ,
                               legend.position = legend.position
                               ,
                               frame.color = frame.color
                               ,
                               label.hub = label.hub
                               #,edge.arrow.size=edge.arrow.size
                               #,edge.thickness=edge.thickness
                             )
                             
                           }
                           text(-1, 1, round(i, 3))
                         }
                       })
                     }
                     if (type.ani == "gif") {
                       setwd(orig_working_directory)
                       require(animation)
                       animation::ani.options(ani.height = size[2],
                                              ani.width = size[1])#, qpdf = '/opt/local/bin/qpdf'
                       animation::saveGIF({
                         POS <- position(net, nv = list_nv[1])
                         
                         for (i in list_nv) {
                           if (fix == TRUE) {
                             plot(
                               Omega
                               ,
                               nv = i
                               ,
                               gr = gr
                               ,
                               ini = POS
                               ,
                               color.vertex = color.vertex
                               ,
                               color.edge = color.edge
                               ,
                               label_v = label_v
                               ,
                               legend.position = legend.position
                               ,
                               frame.color = frame.color
                               ,
                               label.hub = label.hub,
                               movie.name=paste(outdir,'animation.gif',sep='/')
                               #,edge.arrow.size=edge.arrow.size
                               #,edge.thickness=edge.thickness
                             )
                           }
                           else{
                             plot(
                               Omega
                               ,
                               nv = i
                               ,
                               gr = gr
                               ,
                               color.vertex = color.vertex
                               ,
                               color.edge = color.edge
                               ,
                               label_v = label_v
                               ,
                               legend.position = legend.position
                               ,
                               frame.color = frame.color
                               ,
                               label.hub = label.hub
                               #,edge.arrow.size=edge.arrow.size
                               #,edge.thickness=edge.thickness
                             )
                             
                           }
                           text(-1, 1, round(i, 3))
                         }
                       })
                     }
                     if (type.ani == "manual_gif") {
                       par(ask = TRUE)
                       POS <-
                         position(net, nv = list_nv[1])
                       if (is.null(gr)) {
                         gr <- rep(1, dim(Omega@omics_network)[1])
                       }
                       for (i in list_nv) {
                         plot(
                           Omega
                           ,
                           nv = i
                           ,
                           gr = gr
                           ,
                           ini = POS
                           ,
                           color.vertex = color.vertex
                           ,
                           color.edge = color.edge
                           ,
                           label_v = label_v
                           ,
                           legend.position = legend.position
                           ,
                           frame.color = frame.color
                           ,
                           label.hub = label.hub
                           #,edge.arrow.size=edge.arrow.size
                           #,edge.thickness=edge.thickness
                         )
                       }
                       par(ask = FALSE)
                     }
                   }
          )
          

#' Returns the position of edges in the network
#' 
#' Returns the position of edges in the network
#' Retrieve network position for consistent plotting.
#' Utility function to plot networks.
#' 
#' @name position-methods
#' @aliases position position-methods position,omics_network-method
#' @docType methods
#' @section Methods: \describe{
#' 
#' \item{list("signature(net = \"omics_network\")")}{ Returns a matrix with the
#' position of the node. This matrix can then be used as an argument in the
#' plot function. } }
#' 
#' @param net a omics_network object
#' @param nv the level of cutoff at which the analysis should be done
#' @return Matrix with as many rows as the number of edges of the network and three
#' columns (name, xcoord, ycoord).
#' @author Bertrand Frederic, Myriam Maumy-Bertrand.
#' @keywords methods dplots
#' @examples
#' 
#' data(network)
#' position(network)
#' 
setMethod("position", "omics_network", function(net, nv = 0) {
  require(igraph, quietly = TRUE, warn.conflicts = FALSE);on.exit(unloadNamespace("package:igraph"))
  O <- net@omics_network
  Omega <- net@omics_network
  O[abs(O) <= nv] <- 0
  O[abs(O) > nv] <- 1
  
  nom <- 1:dim(O)[1]
  enle <- which(apply(O, 1, sum) + apply(O, 2, sum) == 0)
  if (length(enle) != 0) {
    nom <- nom[-enle]
    O <- O[-enle, -enle]
  }
  
  G <- graph.adjacency(O, weighted = TRUE)
  #    	if(.Platform$OS.type=="unix"){
  #    get.edgelist(G)+1->Q
  #    }else{
  get.edgelist(G) -> Q
  #    	}
  Q[, 1] <- nom[Q[, 1]]
  Q[, 2] <- nom[Q[, 2]]
  
  # dev.new(width=150,heigth=300,xlim=c(min(L[,1]),max(L[,1])))
  #par(mar=c(0,0,0,0), oma=c(0,0,0,0),mai=c(0,0,0,0))
  L <-
    layout.fruchterman.reingold(
      G,
      minx = rep(0, vcount(G)),
      maxx = rep(150, vcount(G)),
      miny = rep(0, vcount(G)),
      maxy = rep(300, vcount(G))
    )
  return(cbind(nom, L))
})
          
#######################################
#######################################
#######################################

#' Find the neighborhood of a set of nodes.
#' 
#' Find the neighborhood of a set of nodes.
#' 
#' 
#' @aliases geneNeighborhood geneNeighborhood-methods
#' geneNeighborhood,omics_network-method
#' @param net a omics_network object
#' @param targets a vector containing the set of nodes
#' @param nv the level of cutoff. Defaut to 0.
#' @param order of the neighborhood. Defaut to `length(net@time_pt)-1`.
#' @param label_v vector defining the vertex labels.
#' @param ini using the ``position'' function, you can
#' fix the position of the nodes.
#' @param frame.color color of the frames.
#' @param label.hub logical ; if TRUE only the hubs are labeled.
#' @param graph plot graph of the network. Defaults to `TRUE`.
#' @param names return names of the neighbors. Defaults to `FALSE`.
#' @return The neighborhood of the targeted genes.
#' @author Bertrand Frederic, Myriam Maumy-Bertrand.
#' @keywords methods
#' @examples
#' 
#' data(Selection)
#' data(infos)
#' #Find probesets for EGR1
#' pbst_EGR1 = infos[infos$hgnc_symbol=="EGR1", "affy_hg_u133_plus_2"]
#' 
#' gene_IDs = infos[match(Selection@name, infos$affy_hg_u133_plus_), "hgnc_symbol"]
#' 
#' data(network)
#' #A nv value can chosen using the cutoff function
#' nv=.11 
#' EGR1<-which(is.element(Selection@name,pbst_EGR1))
#' P<-position(network,nv=nv)
#' 
#' geneNeighborhood(network,targets=EGR1,nv=nv,ini=P,
#' label_v=gene_IDs)
#'           
setMethod("geneNeighborhood", "omics_network"
          , function(net
                     ,
                     targets
                     ,
                     nv = 0
                     ,
                     order = length(net@time_pt) - 1
                     ,
                     label_v = NULL
                     ,
                     ini = NULL
                     ,
                     frame.color = "white"
                     ,
                     label.hub = FALSE
                     #,edge.arrow.size=0.7
                     #,edge.thickness=1
                     ,
                     graph = TRUE
                     ,
                     names = F) {
            require(igraph, quietly = TRUE, warn.conflicts = FALSE);on.exit(unloadNamespace("package:igraph"))
            O <- net@omics_network
            Omega <- net@omics_network
            O[abs(O) <= nv] <- 0
            O[abs(O) > nv] <- 1
            G <- graph.adjacency(O, weighted = TRUE)
            coloring <- colorRamp(c("red", "blue"))
            color <- rep(grey(0.92), length(V(G)))
            if (is.null(label_v)) {
              label_v <- 1:dim(O)[1]
            }
            K <- vector("list", order)
            for (j in order:1) {
              #		if(.Platform$OS.type=="unix"){
              #    N<-neighborhood(G, order=j, nodes=targets-1, mode=c( "out"))
              #    }else{
              N <-
                neighborhood(G,
                             order = j,
                             nodes = targets,
                             mode = c("out"))
              #    	}
              
              
              K[[j]] <- N
              gg <- length(targets)
              
              for (i in 1:gg) {
                #if(.Platform$OS.type=="unix"){
                #color[N[[i]]+1]<-rgb(coloring((j-1)*1/(order-1)),alpha=255,maxColorValue=255)
                #}else{
                color[N[[i]]] <-
                  rgb(coloring((j - 1) * 1 / (order - 1)),
                      alpha = 255,
                      maxColorValue = 255)
                
                #}
                
              }
            }
            if (graph == TRUE) {
              color[targets] <- "green"
              plot(
                net
                ,
                nv = nv
                ,
                ini = ini
                ,
                color.vertex = color
                ,
                frame.color = frame.color
                ,
                label.hub = label.hub
                ,
                label_v = label_v
              )
              legend(
                "topright"
                ,
                legend = paste("order", 1:order)
                ,
                pch = 16
                ,
                col = rgb(coloring((1:order - 1) * 1 / (order - 1)), maxColorValue = 255)
              )
            }
            ind1 <- K[[1]][[1]]
            neighbors <- list(net@name[ind1])
            for (ord in 2:order) {
              neighbors[[ord]] <- list()
              for (nod in 1:gg) {
                neighbPrec <- K[[ord - 1]][[nod]]
                neighbCur <- K[[ord]][[nod]]
                ind <- neighbCur[!(neighbCur %in% neighbPrec)]
                neighbors[[ord]][[nod]] <- net@name[ind]
              }
            }
            if (names) {
              return(neighbors)
            }
            else{
              return(K)
            }
          })


#' Choose the best cutoff
#' 
#' Allows estimating the best cutoff. For a sequence of cutoff, the p value
#' corresponding to each cutoff value of the sequence. Mainly recommended for
#' single time cascade networks. To achieve more sparsity in other settings,
#' please use a fiiting function based on the stability selection or
#' selectboost algorithms.
#' 
#' 
#' @aliases cutoff cutoff-methods cutoff,omics_network-method
#' @param Omega a omics_network object
#' @param sequence a vector corresponding to the sequence of cutoffs that will be tested.
#' @param x_min an integer ; only values over x_min are further retained for performing the test.
#' 
#' @return A list containing two objects : \item{p.value}{the p values
#' corresponding to the sequence of cutoff} \item{p.value.inter}{the smoothed p
#' value vector, using the loess function}
#' @author Bertrand Frederic, Myriam Maumy-Bertrand.
#' @keywords methods
#' @examples
#' 
#' \donttest{
#' 		data(network)
#' 		cutoff(network)
#' 		#See vignette for more details
#' }
#' 
setMethod("cutoff", "omics_network", function(Omega,
                                        sequence = NULL,
                                        x_min = 0) {
  plfit <- function(x = VGAM::rpareto(1000, 10, 2.5)
                    ,
                    method = "limit"
                    ,
                    value = c()
                    ,
                    finite = FALSE
                    ,
                    nowarn = FALSE
                    ,
                    nosmall = FALSE) {
    #init method value to NULL
    vec <- c()
    sampl <- c()
    limit <- c()
    fixed <- c()
    ###########################################################################################
    #
    #  test and trap for bad input
    #
    switch(
      method
      ,
      range = vec <- value
      ,
      sample = sampl <- value
      ,
      limit = limit <- value
      ,
      fixed = fixed <- value
      ,
      argok <- 0
    )
    
    if (exists("argok")) {
      stop("(plfit) Unrecognized method")
    }
    
    if (!is.null(vec) &&
        (!is.vector(vec) || min(vec) <= 1 || length(vec) <= 1)) {
      print(
        paste(
          "(plfit) Error: ''range'' argument must contain a vector > 1; using default."
        )
      )
      vec <- c()
    }
    if (!is.null(sampl) &&
        (!(sampl == floor(sampl)) ||  length(sampl) > 1 || sampl < 2)) {
      print(
        paste(
          "(plfit) Error: ''sample'' argument must be a positive integer > 2; using default."
        )
      )
      sample <- c()
    }
    if (!is.null(limit) && (length(limit) > 1 || limit < 1)) {
      print(paste(
        "(plfit) Error: ''limit'' argument must be a positive >=1; using default."
      ))
      limit <- c()
    }
    if (!is.null(fixed) && (length(fixed) > 1 || fixed <= 0)) {
      print(paste(
        "(plfit) Error: ''fixed'' argument must be a positive >0; using default."
      ))
      fixed <- c()
    }
    
    #  select method (discrete or continuous) for fitting and test if x is a vector
    fdattype <- "unknow"
    if (is.vector(x, "numeric")) {
      fdattype <- "real"
    }
    if (all(x == floor(x)) && is.vector(x)) {
      fdattype <- "integer"
    }
    if (all(x == floor(x)) &&
        min(x) > 1000 && length(x) > 100) {
      fdattype <- "real"
    }
    if (fdattype == "unknow") {
      stop("(plfit) Error: x must contain only reals or only integers.")
    }
    
    #
    #  end test and trap for bad input
    #
    ###########################################################################################
    
    ###########################################################################################
    #   CONTINUOUS CASE
    #  estimate xmin and alpha in the continuous case
    #
    if (fdattype == "real") {
      xmins <- sort(unique(x))
      xmins <- xmins[-length(xmins)]
      
      if (!is.null(limit)) {
        xmins <- xmins[xmins >= limit]
      }
      if (!is.null(fixed)) {
        xmins <- fixed
      }
      if (!is.null(sampl)) {
        xmins <- xmins[unique(round(seq(1, length(xmins), length.out = sampl)))]
      }
      #Vecteur pour stocker les valeurs de D
      dat <- rep(0, length(xmins))
      #Echantillon trie
      z <- sort(x)
      
      #Boucle sur toutes les valeurs de xmin pour lesquelles on calcule D
      for (xm in 1:length(xmins)) {
        #xmin courant
        xmin <- xmins[xm]
        #Exces
        z    <- z[z >= xmin]
        n    <- length(z)
        # estimate alpha-1 using direct MLE
        a    <- n / sum(log(z / xmin))
        # truncate search if nosmall is selected
        if (nosmall) {
          if ((a - 1) / sqrt(n) > 0.1) {
            dat <- dat[1:(xm - 1)]
            print(paste(
              "(plfit) Warning : xmin search truncated beyond",
              xmins[xm - 1]
            ))
            break
          }
        }
        # compute KS statistic
        cx   <- c(0:(n - 1)) / n
        cf   <- 1 - (xmin / z) ^ a
        dat[xm] <- max(abs(cf - cx))
      }
      #min de D
      D     <- min(dat)
      #argmin de D
      xmin  <- xmins[min(which(dat <= D))]
      #exces et estimations finaux
      z     <- x[x >= xmin]
      n     <- length(z)
      alpha <- 1 + n / sum(log(z / xmin))
      
      if (finite) {
        # finite-size correction
        alpha <- alpha * (n - 1) / n + 1 / n
      }
    }
    #
    #  end continuous case
    #
    ###########################################################################################
    
    ###########################################################################################
    #  DISCRETE CASE
    #  estimate xmin and alpha in the discrete case
    #
    if (fdattype == "integer") {
      if (is.null(vec)) {
        vec <-
          seq(1.5, 3.5, .01)
      } # covers range of most practical scaling parameters
      zvec <- VGAM::zeta(vec)
      #On enleve la plus grande valeur
      xmins <- sort(unique(x))
      xmins <- xmins[-length(xmins)]
      
      if (!is.null(limit)) {
        limit <- round(limit)
        xmins <- xmins[xmins >= limit]
      }
      
      if (!is.null(fixed)) {
        xmins <- fixed
      }
      
      if (!is.null(sampl)) {
        xmins <- xmins[unique(round(seq(1, length(xmins), length.out = sampl)))]
      }
      
      if (is.null(xmins) || length(xmins) < 2) {
        stop("(plfit) error: x must contain at least two unique values.")
      }
      
      if (length(which(xmins == 0) > 0)) {
        stop("(plfit) error: x must not contain the value 0.")
      }
      
      xmax <- max(x)
      dat <- matrix(0, nrow = length(xmins), ncol = 2)
      z <- x
      for (xm in 1:length(xmins)) {
        xmin <- xmins[xm]
        z    <- z[z >= xmin]
        n    <- length(z)
        # estimate alpha via direct maximization of likelihood function
        #  vectorized version of numerical computation
        # matlab: zdiff = sum( repmat((1:xmin-1)',1,length(vec)).^-repmat(vec,xmin-1,1) ,1);
        if (xmin == 1) {
          zdiff <- rep(0, length(vec))
        } else{
          zdiff <-
            apply(rep(t(1:(xmin - 1)), length(vec)) ^ -t(kronecker(t(
              array(1, xmin - 1)
            ), vec)), 2, sum)
        }
        # matlab: L = -vec.*sum(log(z)) - n.*log(zvec - zdiff);
        L <- -vec * sum(log(z)) - n * log(zvec - zdiff)
        
        I <- which.max(L)
        # compute KS statistic
        fit <-
          cumsum((((xmin:xmax) ^ -vec[I])) / (zvec[I] - sum((1:(
            xmin - 1
          )) ^ -vec[I])))
        cdi <-
          cumsum(hist(z, c(
            min(z) - 1, (xmin + .5):xmax, max(z) + 1
          ), plot = FALSE)$counts / n)
        dat[xm, ] <- c(max(abs(fit - cdi)), vec[I])
      }
      D     <- min(dat[, 1])
      I     <- which.min(dat[, 1])
      xmin  <- xmins[I]
      n     <- sum(x >= xmin)
      alpha <- dat[I, 2]
      
      if (finite) {
        alpha <- alpha * (n - 1) / n + 1 / n # finite-size correction
      }
      
    }
    #
    #  end discrete case
    #
    ###########################################################################################
    
    #  return xmin, alpha and D in a list
    return(list(
      xmin = xmin,
      alpha = alpha,
      D = D
    ))
  }
  
  ###########################################################################################
  ###########################################################################################
  ###########################################################################################
  ###########################################################################################
  ###########################################################################################
  ###########################################################################################
  plpva <-
    function(x = VGAM::rpareto(1000, 10, 2.5),
             xmin = 1,
             method = "limit",
             value = c(),
             Bt = 1000,
             quiet = FALSE,
             vec = seq(1.5, 3.5, .01)) {
      #init method value to NULL
      sampl <- c()
      limit <- c()
      ###########################################################################################
      #
      #  test and trap for bad input
      #
      switch(method,
             sample = sampl <- value,
             limit = limit <- value,
             argok <- 0)
      
      if (exists("argok")) {
        stop("(plvar) Unrecognized method")
      }
      
      if (!is.null(vec) &&
          (!is.vector(vec) || min(vec) <= 1 || length(vec) <= 1)) {
        print(
          paste(
            "(plvar) Error: ''range'' argument must contain a vector > 1; using default."
          )
        )
        vec <- c()
      }
      if (!is.null(sampl) &&
          (!(sampl == floor(sampl)) ||  length(sampl) > 1 || sampl < 2)) {
        print(
          paste(
            "(plvar) Error: ''sample'' argument must be a positive integer > 2; using default."
          )
        )
        sample <- c()
      }
      if (!is.null(limit) && (length(limit) > 1 || limit < 1)) {
        print(paste(
          "(plvar) Error: ''limit'' argument must be a positive >=1; using default."
        ))
        limit <- c()
      }
      if (!is.null(Bt) && (!is.vector(Bt) || Bt <= 1 ||
                           length(Bt) > 1)) {
        print(paste(
          "(plvar) Error: ''Bt'' argument must be a positive value > 1; using default."
        ))
        vec <- c()
      }
      
      #  select method (discrete or continuous) for fitting and test if x is a vector
      fdattype <- "unknow"
      if (is.vector(x, "numeric")) {
        fdattype <- "real"
      }
      if (all(x == floor(x)) && is.vector(x)) {
        fdattype <- "integer"
      }
      if (all(x == floor(x)) &&
          min(x) > 1000 && length(x) > 100) {
        fdattype <- "real"
      }
      if (fdattype == "unknow") {
        stop("(plfit) Error: x must contain only reals or only integers.")
      }
      
      N   <- length(x)
      nof <- rep(0, Bt)
      
      if (!quiet) {
        print("Power-law Distribution, parameter error computation")
        print("Warning: This can be a slow computation; please be patient.")
        print(paste(" n =", N, "xmin =", xmin, "- reps =", Bt, fdattype))
      }
      #
      #  end test and trap for bad input
      #
      ###########################################################################################
      
      ###########################################################################################
      #
      #  estimate xmin and alpha in the continuous case
      #
      if (fdattype == "real") {
        # compute D for the empirical distribution
        z     <- x[x >= xmin]
        nz   <- length(z)
        y     <- x[x < xmin]
        ny   <- length(y)
        alpha <- 1 + nz / sum(log(z / xmin))
        cz    <- (0:(nz - 1)) / nz
        cf    <- 1 - (xmin / sort(z)) ^ (alpha - 1)
        gof   <- max(abs(cz - cf))
        pz    <- nz / N
        
        # compute distribution of gofs from semi-parametric bootstrap
        # of entire data set with fit
        for (B in 1:length(nof)) {
          # semi-parametric bootstrap of data
          n1 <- sum(runif(N) > pz)
          q1 <- y[sample(ny, n1, replace = TRUE)]
          n2 <- N - n1
          q2 <- xmin * (1 - runif(n2)) ^ (-1 / (alpha - 1))
          q  <- sort(c(q1, q2))
          
          # estimate xmin and alpha via GoF-method
          qmins <- sort(unique(q))
          qmins <- qmins[-length(qmins)]
          if (!is.null(limit)) {
            qmins <- qmins[qmins <= limit]
          }
          if (!is.null(sampl)) {
            qmins <- qmins[unique(round(seq(1, length(qmins), length.out = sampl)))]
          }
          dat <- rep(0, length(qmins))
          for (qm in 1:length(qmins)) {
            qmin <- qmins[qm]
            zq   <- q[q >= qmin]
            nq   <- length(zq)
            a    <- nq / sum(log(zq / qmin))
            cq   <- (0:(nq - 1)) / nq
            cf   <- 1 - (qmin / zq) ^ a
            dat[qm] <- max(abs(cq - cf))
          }
          if (!quiet) {
            print(paste(B, sum(nof[1:B] >= gof) / B))
          }
          # store distribution of estimated gof values
          nof[B] <- min(dat)
        }
      }
      
      ###########################################################################################
      #
      #  estimate xmin and alpha in the discrete case
      #
      if (fdattype == "integer") {
        if (is.null(vec)) {
          vec <-
            seq(1.5, 3.5, .01)
        } # covers range of most practical scaling parameters
        zvec <- VGAM::zeta(vec)
        
        # compute D for the empirical distribution
        z     <- x[x >= xmin]
        nz   <- length(z)
        xmax <- max(z)
        y     <- x[x < xmin]
        ny   <- length(y)
        
        if (xmin == 1) {
          zdiff <- rep(0, length(vec))
        } else{
          zdiff <-
            apply(rep(t(1:(xmin - 1)), length(vec)) ^ -t(kronecker(t(
              array(1, xmin - 1)
            ), vec)), 2, sum)
        }
        # matlab: L = -vec.*sum(log(z)) - n.*log(zvec - zdiff);
        L <- -vec * sum(log(z)) - nz * log(zvec - zdiff)
        
        I <- which.max(L)
        alpha <- vec[I]
        
        fit <- cumsum((((xmin:xmax) ^ -alpha))
                      / (zvec[I]
                         - (if (xmin == 1)
                           0
                           else
                             sum((
                               1:(xmin - 1)
                             ) ^ -alpha))))
        
        hist = aggregate(z, list(z), length)
        cdi = rep(0, xmax - xmin + 1)
        cdi[hist$Group.1 - xmin + 1] = hist$x / nz
        cdi = cumsum(cdi)
        
        gof <- max(abs(fit - cdi))
        pz  <- nz / N
        
        mmax <- 20 * xmax
        pdf <- c(rep(0, xmin - 1), (((xmin:mmax) ^ -alpha))
                 / (zvec[I]
                    - (if (xmin == 1)
                      0
                      else
                        sum((
                          1:(xmin - 1)
                        ) ^ -alpha))))
        cdf <- cbind(1:(mmax + 1), c(cumsum(pdf), 1))
        
        # compute distribution of gofs from semi-parametric bootstrap
        # of entire data set with fit
        for (B in 1:Bt) {
          # semi-parametric bootstrap of data
          n1 <- sum(runif(N) > pz)
          q1 <- y[sample(ny, n1, replace = TRUE)]
          n2 <- N - n1
          # simple discrete zeta generator
          r2 <- sort(runif(n2))
          d <- 1
          q2 <- rep(0, n2)
          k <- 1
          for (i in xmin:(mmax + 1)) {
            while ((d <= length(r2)) && (r2[d] <= cdf[i, 2])) {
              d <- d + 1
            }
            if (k <= d - 1)
              q2[k:(d - 1)] <- i
            k <- d
            if (k > n2) {
              break
            }
          }
          q <- c(q1, q2)
          ########################################
          # estimate xmin and alpha via GoF-method
          qmins <- sort(unique(q))
          qmins <- qmins[-length(qmins)]
          if (!is.null(limit)) {
            qmins <- qmins[qmins <= limit]
          }
          if (!is.null(sampl)) {
            qmins <-
              qmins[sort(unique(round(
                seq(1, length(qmins), length.out = sampl)
              )))]
          }
          dat   <- rep(0, length(qmins))
          qmax  <- max(q)
          zq    <- q
          
          if (length(qmins) > 0)
            for (qm in 1:length(qmins)) {
              qmin <- qmins[qm]
              zq   <- zq[zq >= qmin]
              nq    <- length(zq)
              if (nq > 1) {
                # vectorized version of numerical computation
                if (qmin == 1) {
                  #WARNING
                  #zdiff <- rep(1,length(vec))
                  #correction added the 6/10/2009 (following Naoki Masuda for plfit)
                  zdiff <- rep(0, length(vec))
                } else{
                  zdiff <-
                    apply(rep(t(1:(
                      qmin - 1
                    )), length(vec)) ^ -t(kronecker(t(
                      array(1, qmin - 1)
                    ), vec)), 2, sum)
                }
                L <- -vec * sum(log(zq)) - nq * log(zvec - zdiff)
                
                I <- which.max(L)
                # compute KS statistic
                fit <- cumsum((((
                  qmin:qmax
                ) ^ -vec[I]))
                / (zvec[I]
                   - (if (qmin == 1)
                     0
                     else
                       sum((1:(qmin - 1)) ^ -vec[I]
                       ))))
                
                hist = aggregate(zq, list(zq), length)
                cdi = rep(0, qmax - qmin + 1)
                cdi[hist$Group.1 - qmin + 1] = hist$x / nq
                cdi = cumsum(cdi)
                
                dat[qm] <- max(abs(fit - cdi))
              } else{
                dat[qm] <- -Inf
              }
            }
          if (!quiet) {
            print(paste(B, "-", sum(nof[1:B] >= gof) / B))
          }
          # -- store distribution of estimated gof values
          nof[B] <- min(c(dat, Inf))
          
          
        }
        
        ####################################
        
      }
      #
      #  end discrete case
      #
      ###########################################################################################
      #  return p and gof in a list
      p <- sum(nof >= gof) / length(nof)
      return(list(p = p, gof = gof))
    }
  
  sequence_test <- sequence
  
  #   if(is.null(x_min)){
  #     x_min<-round(dim(Omega@omics_network)[1]*0.02)
  #   }
  if (is.null(sequence_test)) {
    sequence_test <-
      seq(0, min(max(abs(
        Omega@omics_network * 0.9
      )), 0.4), length.out = 10)
  }
  #
  #install_github('poweRlaw', 'csgillespie')
  #cutoff courant
  cc <- 0
  #compteur
  u <- 0
  #pvaleur
  f <- rep(0, length(sequence_test))
  
  #Boucle sur les valeurs candidates pour le cutoff
  print("This computation may be long")
  for (cc in sequence_test) {
    u <- u + 1
    print(paste(u, length(sequence_test), sep = "/"))
    O <- Omega@omics_network
    #On garde les liens d'intensite superieure au cutoff
    O[abs(O) > cc] <- 1
    O[abs(O) <= cc] <- 0
    #nombre de liens pour chaque gene
    liste <- apply(O, 1, sum) + 1
    if (length(which(liste > round(x_min) + 1)) < 4) {
      break
    }
    #On prend les genes ayant plus d'un lien
    xmin = plpva(liste, quiet = TRUE)$p
    f[u] <- xmin
  }
  #index du premier element superieur a 0.25
  #   p<-which(sequence_test>0.25)[1]
  #   #On ne conserve que ceux qui sont inferieurs a 0.25
  #   f<-f[1:p]
  print(f)
  K <- 1:length(f)
  f2 <- f
  #interpolation des points de f pour lisser
  f <- predict(loess(f ~ K))
  par(mar = c(5, 4, 0.2, 0.2))
  plot(
    sequence_test
    ,
    f
    ,
    xlab = "cutoff sequence"
    ,
    ylab = "scale-freeness test p-value"
    ,
    type = "l"
    ,
    lwd = 2
    ,
    ylim = c(min(f) - 0.02, max(f))
  )
  abline(h = 0.1, col = "red")
  rect(0.10, 0.10, 0.15, 1, col = rgb(0, 1, 0, 0.5))
  rect(0.05, 0.10, 0.10, 1, col = rgb(0.5, 1, 0, 0.4))
  rect(-0.5, 0.10, 0.05, 1, col = rgb(1, 0, 0, 0.4))
  rect(0.15, 0.10, 0.5, 1, col = rgb(0.5, 1, 0, 0.4))
  abline(h = 0.1, col = "red", lwd = 3)
  legend(
    "bottomright"
    ,
    lty = c(1, 0, 0, 0)
    ,
    col = c("red"
            , rgb(0, 1, 0, 0.5)
            , rgb(0.5, 1, 0, 0.4)
            , rgb(1, 0, 0, 0.4))
    ,
    pch = c(NA, 15, 15, 15)
    ,
    legend = c(
      "p-value=0.1:  above this line, scale-freeness may be assumed"
      ,
      "best area of choice (determined by simulation)"
      ,
      "less recommended area of choice (determined by simulation)"
      ,
      "area of choice to be avoided (determined by simulation)"
    )
    ,
    lwd = c(3, 5, 5, 5)
    ,
    cex = 1
  )
  
  return(list(
    p.value = f2,
    p.value.inter = f,
    sequence = sequence_test
  ))
})


#' Some basic criteria of comparison between actual and inferred network.
#' 
#' Allows comparison between actual and inferred network.
#' 
#' 
#' @name compare-methods
#' @aliases compare-methods compare compare,omics_network,omics_network,numeric-method
#' @docType methods
#' @return A vector containing : sensitivity, predictive positive value, the
#' usual F-score (2*ppv*sens/(sppvpe+sens)), the 1/2 ponderated Fscore
#' ((1+0.5^2)*ppv*sens/(ppv/4+sens)) and the 2 ponderated Fscore
#' ((1+2^2)*ppv*sens/(ppv*4+sens)).
#' @section Methods: \describe{
#' 
#' \item{list("signature(Net = \"omics_network\", Net_inf = \"omics_network\", nv =
#' \"numeric\")")}{ \describe{ \item{Net}{ A omics_network object containing the
#' actual network.  } \item{Net_inf}{ A omics_network object containing the inferred
#' network.  } \item{nv}{ A number that indicates at which level of cutoff the
#' comparison should be done.  } } }
#' 
#' }
#' @param Net A omics_network object containing the
#' actual network.
#' @param Net_inf A omics_network object containing the inferred
#' network.
#' @param nv A number that indicates at which level of cutoff the
#' comparison should be done.
#' @author Bertrand Frederic, Myriam Maumy-Bertrand.
#' @examples
#' 
#' data(Net)
#' data(Net_inf_PL)
#' 
#' #Comparing true and inferred networks
#' Crit_values=NULL
#' 
#' #Here are the cutoff level tested
#' test.seq<-seq(0,max(abs(Net_inf_PL@omics_network*0.9)),length.out=200)
#' for(u in test.seq){
#' 	Crit_values<-rbind(Crit_values,Patterns::compare(Net,Net_inf_PL,u))
#' }
#' matplot(test.seq,Crit_values,type="l",ylab="Criterion value",xlab="Cutoff level",lwd=2)
#' legend(x="topleft", legend=colnames(Crit_values), lty=1:5,col=1:5,ncol=2,cex=.9)
#' 
setMethod("compare",c("omics_network","omics_network","numeric"),function(Net,
                                                              Net_inf,
                                                              nv=1){
  N1<-Net@omics_network
  N2<-Net_inf@omics_network
  N1[abs(N1)>0]<-1
  N1[abs(N1)<=0]<-0
  N2[abs(N2)>nv]<-1
  N2[abs(N2)<=nv]<-0
  Nb<-sum(N1)
  sens<-0
  if(sum(N1==1)){
    sens<-sum((N1-2*N2)==-1)/sum(N1==1)
  }
  spe<-0
  if(sum(N2==1)){
    spe<-sum((N1-2*N2)==-1)/sum(N2==1)
  }
  Fscore<-0
  Fscore12<-0
  Fscore2<-0
  if(sens !=0){
    Fscore<-2*spe*sens/(spe+sens)
    Fscore12<-(1+0.5^2)*spe*sens/(spe/4+sens)
    Fscore2<-(1+2^2)*spe*sens/(spe*4+sens)
    crits=c(sens,spe,Fscore,Fscore12,Fscore2)
    names(crits) <- c("Sensitivity","PPV","F-score","F-score 1/2","F-score 2")
    return(crits)
  }
  crits=c(sens,spe,Fscore,Fscore12,Fscore2)
  names(crits) <- c("Sensitivity","PPV","F-score","F-score 1/2","F-score 2")
  return(crits)
}
)



#' @rdname plot-methods
setMethod("plot"
          , "omics_network"
          , function(x
                     ,
                     y
                     ,
                     choice = "omics_network"
                     ,
                     nv = 0
                     ,
                     gr = NULL
                     ,
                     ini = NULL
                     ,
                     color.vertex = NULL
                     ,
                     color.edge = NULL
                     ,
                     video = TRUE
                     ,
                     weight.node = NULL
                     ,
                     ani = FALSE
                     ,
                     size = c(2000, 1000)
                     ,
                     label_v = 1:dim(x@omics_network)[1]
                     ,
                     horiz = TRUE
                     ,
                     legend.position = "topleft"
                     ,
                     frame.color = "black"
                     ,
                     label.hub = FALSE
                     ,
                     nround = 2
                     ,
                     ani.img.name = "Rplot"
                     ,
                     ani.imgdir = "images"
                     ,
                     ani.htmlfile = "index.html"
                     ,
                     outdir
                     ,
                     ani.group.legend = "Cluster"
                     ,
                     layout = ini
                     ,
                     alpha = 1
                     ,
                     pixmap.color = terrain.colors(20)
                     #,edge.arrow.size=0.6*(1+size.ed)
                     #,edge.thickness=1
                     ,
                     ...)
          {
            if (choice == "F") {
              require(plotrix, quietly = TRUE)
              F <- x@F
              sizeF <- dim(F)[1]
              nF <- dim(F)[3]
              ngrp = sqrt(dim(F)[3])
              ymax <- max(F)
              coloring <- rainbow(ngrp)
              #  par(mfrow=c(ngrp,ngrp))
              FF = NULL
              for (i in 1:ngrp) {
                FFa = NULL
                for (j in 1:ngrp) {
                  FFa = cbind(FFa, round(F[, , (i - 1) * ngrp + j], nround))
                  #par(mai=c(0.1,0.1,0.1,0.1))
                  #color2D.matplot(x=round(F[,,(i-1)*ngrp+j],2),cs1=c(0,1,1),cs2=c(1,1,0),cs3=0,show.values=TRUE,axes=FALSE,main="",xlab="",ylab="",show.legend=TRUE)
                }
                FF = rbind(FF, FFa)
              }
              par(
                mar = c(0, 0, 0, 0),
                oma = c(0, 0, 0, 0),
                mai = c(0, 0, 0, 0)
              )
              color2D.matplot(
                x = FF,
                cs1 = c(0, .5, 1),
                cs2 = c(.5, 1, 0),
                cs3 = c(1, 2, 0),
                show.values = nround,
                axes = FALSE,
                main = "",
                xlab = "",
                ylab = "",
                show.legend = FALSE
              )
              abline(h = sizeF * (1:(ngrp - 1)), lwd = 3, col="grey50")
              abline(v = sizeF * (1:(ngrp - 1)), lwd = 3, col="grey50")
            }
            if (choice == "Fpixmap") {
              require(plotrix, quietly = TRUE)
              F <- x@F
              nF <- dim(F)[3]
              ngrp = sqrt(dim(F)[3])
              ymax <- max(F)
              coloring <- rainbow(ngrp)
              #  par(mfrow=c(ngrp,ngrp))
              FF = NULL
              for (i in 1:ngrp) {
                FFa = NULL
                for (j in 1:ngrp) {
                  FFa = cbind(FFa, round(F[, , (i - 1) * ngrp + j], nround))
                  #par(mai=c(0.1,0.1,0.1,0.1))
                  #color2D.matplot(x=round(F[,,(i-1)*ngrp+j],2),cs1=c(0,1,1),cs2=c(1,1,0),cs3=0,show.values=TRUE,axes=FALSE,main="",xlab="",ylab="",show.legend=TRUE)
                }
                FF = rbind(FF, FFa)
              }
              par(
                mar = c(0, 0, 0, 0),
                oma = c(0, 0, 0, 0),
                mai = c(0, 0, 0, 0)
              )
              x.temp <- pixmap::pixmapGrey(data = FF, cellres = c(2, 2))
              plot(x.temp)
              #              abline(h=4*(1:(ngrp-1)),lwd=3)
              #              abline(v=4*(1:(ngrp-1)),lwd=3)
              rm(x.temp)
            }
            if (choice == "omics_network") {
              coloring <- 0
              O <- x@omics_network
              Omega <- x@omics_network
              O[abs(O) <= nv] <- 0
              O[abs(O) > nv] <- 1
              F <- x@F
              nF <- dim(F)[3]
              ngrp = sqrt(dim(F)[3])
              
              
              
              require(igraph, quietly = TRUE, warn.conflicts = FALSE);on.exit(unloadNamespace("package:igraph"))
              if (is.null(gr)) {
                gr <- rep(1, dim(O)[1])
              }
              if (is.null(color.vertex)) {
                color.vertex <- rainbow(length(unique(gr)), alpha = alpha)
              }
              coloring <- 1
              
              nom <- 1:dim(O)[1]
              enle <- which(apply(O, 1, sum) + apply(O, 2, sum) == 0)
              if (length(enle) != 0) {
                nom <- nom[-enle]
                
                O <- O[-enle, -enle]
                
                if (!is.null(weight.node)) {
                  weight.node <- weight.node[-enle]
                }
              }
              
              if (!is.null(ini)) {
                ini <- ini[which(ini[, 1] %in% nom), ]
              }
              #              nom2<-nom
              nom2 <- label_v[nom]
              if (!is.null(weight.node)) {
                size <- weight.node
              } else
              {
                size <- apply(O, 1, sum)
              }
              
              size <- size / max(size) * 10
              size[size < 2] <- 2
              
              if (label.hub == TRUE) {
                nom2[which(size < 3)] <- NA
              }
              
              
              G <- igraph::graph.adjacency(O, weighted = TRUE)
              
              #    	if(.Platform$OS.type=="unix"){
              #    get.edgelist(G)+1->Q
              #    }else{
              igraph::get.edgelist(G) -> Q
              #    	}
              Q[, 1] <- nom[Q[, 1]]
              Q[, 2] <- nom[Q[, 2]]
              size.ed <- rep(0, dim(Q)[1])
              
              for (i in 1:dim(Q)[1]) {
                size.ed[i] <- abs(Omega[Q[i, 1], Q[i, 2]])
              }
              #fpl<-function(x){3*exp(8*x)/(14+exp(8*x))}
              size.ed <- size.ed / max(size.ed)
              #size.ed<-fpl(size.ed)
              #size.ed<-size.ed
              
              if (ani == TRUE) {
                if(missing(outdir)){stop("An output directory is required")}
                #We save the user working directory.
                orig_working_directory <- getwd()
                
                #We restore the user working directory on function exit.
                on.exit(setwd(orig_working_directory))
                
                #Set the working directory to the outdir directory
                setwd(outdir)
                
                if (is.null(gr)) {stop("Need of groups")
                }
                T <- length(unique(gr))#length(x@time_pt)
                require(animation)
                animation::ani.options(ani.height = size[2],
                                       ani.width = size[1])
                
                if (is.null(ini)) {
                  ini <- position(x, nv)[]
                }
                
                
                #ini<-ini[which(ini[,1] %in% nom),]
                ini <- ini[, 2:3]
                gr2 <- gr[nom]
                
                animation::saveHTML({
                  for (i in c(0.999, 0.99, 0.98, 0.97, 0.96, 0.95)) {
                    plot(
                      G
                      ,
                      layout = ini
                      ,
                      vertex.size = size
                      ,
                      edge.arrow.size = 0.6 * (1 + size.ed)
                      ,
                      edge.width = size.ed * 5
                      ,
                      edge.arrow.width = 0.5 * (1 + size.ed)
                      #,edge.width=size.ed*5*edge.thickness
                      #,edge.arrow.size=edge.arrow.size*size.ed*edge.thickness
                      #,edge.arrow.width=edge.arrow.size*size.ed*edge.thickness
                      ,
                      vertex.color = grey(i)
                      ,
                      vertex.label.cex = 1
                      ,
                      edge.color = grey(i)
                      ,
                      asp = 0
                      ,
                      vertex.label = nom2
                      ,
                      vertex.frame.color = frame.color
                    )
                  }
                  for (i in 1:(T)) {
                    Ver.col <- rep(grey(0.95), length(nom))
                    Edge.col <- rep(grey(0.95), dim(Q)[1])
                    mms <- 20
                    PP <- 1:mms
                    for (p in PP) {
                      coul <- col2rgb(color.vertex[i])
                      indi <- which(gr2 == i)
                      Ver.col[indi] <-
                        rgb(
                          coul[1],
                          coul[2],
                          coul[3],
                          alpha = (p - 1) * 255 / (mms),
                          maxColorValue = 255
                        )
                      
                      
                      for (k in 1:i) {
                        coul <- col2rgb(color.vertex[k])
                        indi3 <- which(gr[Q[, 1]] == k &
                                         gr[Q[, 2]] == i)
                        
                        Edge.col[indi3] <-
                          rgb(
                            coul[1],
                            coul[2],
                            coul[3],
                            alpha = 255 - (p - 1) * 255 / (mms),
                            maxColorValue = 255
                          )
                      }
                      
                      if (i != T) {
                        #                        for(k in 1:ngrp){
                        for (k in 1:i) {
                          #                          for(k2 in (1:ngrp)[-k]){
                          for (k2 in min((i + 1), T):T) {
                            coul <- col2rgb(color.vertex[k])
                            indi3 <-
                              which(gr[Q[, 1]] == k & gr[Q[, 2]] == k2)
                            
                            
                            
                            
                            al <- round(mms / abs(k2 - k))
                            ald <- min(1 + al * abs(i - k), mms)
                            if (k2 == (i + 1)) {
                              alf <- mms
                            }
                            else{
                              alf <- min(al * (i - k + 1), mms)
                            }
                            #alf<-max(1,min(al*abs(i-k+1),mms))}
                            alseq <- seq(25, 255, length.out = mms)
                            #                            cat(al," ");cat(ald," ");cat(alf," ");cat(mms," ");cat(alseq," ");cat("\n")
                            alphaseq2 <-
                              seq(alseq[ald], alseq[alf], length.out = mms)
                            
                            Edge.col[indi3] <-
                              rgb(coul[1],
                                  coul[2],
                                  coul[3],
                                  alpha = alphaseq2[p],
                                  maxColorValue = 255)
                          }
                        }
                      }
                      
                      plot(
                        G
                        ,
                        layout = ini
                        ,
                        vertex.size = size
                        #,edge.width=size.ed*5*edge.thickness
                        #,edge.arrow.size=edge.arrow.size*size.ed*edge.thickness
                        #,edge.arrow.width=edge.arrow.size*size.ed*edge.thickness
                        ,
                        edge.arrow.size = 0.6 * (1 + size.ed)
                        ,
                        edge.width = size.ed * 5
                        ,
                        edge.arrow.size = 0.5 * (1 + size.ed)
                        ,
                        vertex.color = Ver.col
                        ,
                        vertex.label.cex = 1
                        ,
                        edge.color = Edge.col
                        ,
                        asp = 0
                        ,
                        vertex.label = nom2
                        ,
                        vertex.frame.color = frame.color
                      )
                      if (length(unique(gr[nom])) > 1) {
                        legend(
                          legend.position
                          ,
                          horiz = horiz
                          ,
                          pch = 21
                          ,
                          pt.bg = color.vertex[unique(gr[nom])[order(unique(gr[nom]))]]
                          ,
                          col = frame.color
                          ,
                          legend = paste(ani.group.legend, unique(gr[nom])[order(unique(gr[nom]))])
                        )
                      }
                    }
                    
                    
                    
                  }
                  
                  
                }, img.name = ani.img.name, imgdir = ani.imgdir, htmlfile = ani.htmlfile)
                
                
              }
              else{
                par(
                  mar = c(0, 0, 0, 0),
                  oma = c(0, 0, 0, 0),
                  mai = c(0, 0, 0, 0)
                )
                if (length(unique(gr[nom])) > 1 &
                    !is.null(color.edge)) {
                  trr <- color.edge[gr[Q[, 1]]]
                } else {
                  if (length(unique(gr[nom])) > 1 &
                      is.null(color.edge)) {
                    trr <-
                      color.vertex[gr[Q[, 2]]]
                    #cat(length(unique(gr[nom])) > 1 & is.null(color.edge))
                  } else {
                    trr <- "grey"
                    #cat(length(unique(gr[nom])) > 1 & is.null(color.edge))
                  }
                }
                
                if (is.null(ini)) {
                  L <- position(x, nv)[, 2:3]
                  #maxL<-c(max(abs(L[,1])),max(abs(L[,2])))
                  #matMaxL<-matrix(rep(maxL,nrow(L)),ncol=2,byrow=T)
                  #L<-L/matMaxL
                  if (coloring == 0) {
                    plot(
                      G
                      ,
                      layout = L
                      ,
                      vertex.size = size
                      #,edge.arrow.size=edge.arrow.size*size.ed*edge.thickness
                      #,edge.width=size.ed*5*edge.thickness
                      #,edge.arrow.width=edge.arrow.size*size.ed*edge.thickness
                      ,
                      edge.arrow.size = 0.6 * (1 + size.ed)
                      ,
                      edge.width = size.ed * 5
                      ,
                      edge.arrow.width = 0.5 * (1 + size.ed)
                      ,
                      vertex.color = color.vertex[nom]
                      ,
                      vertex.label.cex = 1
                      ,
                      edge.color = trr
                      ,
                      asp = 0
                      ,
                      vertex.label = nom2
                      ,
                      vertex.frame.color = frame.color
                    )
                  }
                  else{
                    plot(
                      G
                      ,
                      layout = L
                      ,
                      vertex.size = size
                      #,edge.width=size.ed*5*edge.thickness
                      #,edge.arrow.size=edge.arrow.size*size.ed*edge.thickness
                      #,edge.arrow.width=edge.arrow.size*size.ed*edge.thickness
                      ,
                      edge.arrow.size = 0.6 * (1 + size.ed)
                      ,
                      edge.width = size.ed * 5
                      ,
                      edge.arrow.width = 0.5 * (1 + size.ed)
                      ,
                      vertex.color = color.vertex[gr[nom]]
                      ,
                      vertex.label.cex = 1
                      ,
                      edge.color = trr
                      ,
                      asp = 0
                      ,
                      vertex.label = nom2
                      ,
                      vertex.frame.color = frame.color
                    )
                  }
                  
                  if (length(unique(gr[nom])) > 1) {
                    legend(
                      legend.position
                      ,
                      horiz = horiz
                      ,
                      pch = 21
                      ,
                      pt.bg = color.vertex[unique(gr[nom])[order(unique(gr[nom]))]]
                      ,
                      col = frame.color
                      ,
                      legend = paste(ani.group.legend, unique(gr[nom])[order(unique(gr[nom]))])
                    )
                  }
                }
                else{
                  L <- ini[, 2:3]
                  if (coloring == 0) {
                    plot(
                      G
                      ,
                      layout = L
                      ,
                      vertex.size = size
                      #,edge.arrow.size=edge.arrow.size*size.ed*edge.thickness
                      #,edge.width=size.ed*5*edge.thickness
                      #,edge.arrow.width=edge.arrow.size*size.ed*edge.thickness
                      ,
                      edge.arrow.size = 0.6 * (1 + size.ed)
                      ,
                      edge.width = size.ed * 5
                      ,
                      edge.arrow.width = 0.5 * (1 + size.ed)
                      ,
                      vertex.color = color.vertex[nom]
                      ,
                      vertex.label.cex = 1
                      ,
                      edge.color = trr
                      ,
                      asp = 0
                      ,
                      vertex.label = nom2
                      ,
                      vertex.frame.color = frame.color
                    )
                  }
                  else{
                    plot(
                      G
                      ,
                      layout = L
                      ,
                      vertex.size = size
                      #,edge.width=size.ed*5*edge.thickness
                      #,edge.arrow.size=edge.arrow.size*size.ed*edge.thickness
                      #,edge.arrow.width=edge.arrow.size*size.ed*edge.thickness
                      ,
                      edge.arrow.size = 0.6 * (1 + size.ed)
                      ,
                      edge.width = size.ed * 5
                      ,
                      edge.arrow.width = 0.5 * (1 + size.ed)
                      ,
                      vertex.color = color.vertex[gr[nom]]
                      ,
                      vertex.label.cex = 1
                      ,
                      edge.color = trr
                      ,
                      asp = 0
                      ,
                      vertex.label = nom2
                      ,
                      vertex.frame.color = frame.color
                    )
                  }
                  if (length(unique(gr[nom])) > 1) {
                    legend(
                      legend.position
                      ,
                      horiz = horiz
                      ,
                      pch = 21
                      ,
                      pt.bg = color.vertex[unique(gr[nom])[order(unique(gr[nom]))]]
                      ,
                      col = frame.color
                      ,
                      legend = paste(ani.group.legend, unique(gr[nom])[order(unique(gr[nom]))])
                    )
                  }
                  L <- ini
                }
              }
              invisible(G)
            }
          })