R/miRsponge.R

Defines functions moduleSurvival moduleFEA moduleDEA netModule spongeValidate spongeMethod integrateMethod cernia muTaME hermes ppc sppc pc miRHomology mcode mcode.post.process mcode.fluff.complex mcode.find.complexex mcode.find.complex mcode.vertex.weighting cluster.save cluster dtHybrid recommendation graphWeights parMM predCor combpvalue calCMI querymiRTargetbinding

Documented in integrateMethod moduleDEA moduleFEA moduleSurvival netModule querymiRTargetbinding spongeMethod spongeValidate

## Query miRNA-target interactions by combining expression data and putative miRNA-target
## interactions
querymiRTargetbinding <- function(ExpData, miRTarget) {

    ExpDataNames <- c(as.matrix(ExpData[1, ]))
    miRTarget <- as.matrix(miRTarget)
    miRTargetCandidate <- miRTarget[intersect(which(miRTarget[, 1] %in% ExpDataNames),
        which(miRTarget[, 2] %in% ExpDataNames)), ]

    return(miRTargetCandidate)
}

## Internal functions (calCMI, combpvalue, predCor) of hermes method are from the website:
## http://califano.c2b2.columbia.edu/hermes Copyright Columbia University in the City of New
## York. You may not use this file except in compliance with the License. You may obtain a copy
## of the License at http://califano.c2b2.columbia.edu/hermes-license

## calCMI: Calculate conditional mutual information
calCMI <- function(dmat) {

    # rank data
    N <- nrow(dmat)
    cdim <- ncol(dmat)
    idat <- apply(dmat, 2, order)

    ydat <- idat
    sapply(seq_len(cdim), function(i) ydat[idat[, i], i] <- seq_len(N))

    ddim <- 2^cdim
    dim2 <- 2 * cdim

    # init
    rm(dmat, idat)
    poc <- vector()
    kon <- vector()
    cmi <- 0
    npar <- 1
    poc[1] <- 1
    kon[1] <- N
    poradi <- seq_len(N)
    NN <- matrix(0, nrow = 1, ncol = ddim)
    marg <- matrix(0, nrow = 8 * ddim, ncol = dim2)
    marg[1, ] <- c(rep(1, cdim), rep(N, cdim))
    Imm <- matrix(c(0, 1), ncol = 1)  # Binary matrix for all combinations
    for (d in seq(2, cdim)) {
        Imm <- rbind(cbind(matrix(0, nrow = nrow(Imm), ncol = 1), Imm),
            cbind(matrix(1, nrow = nrow(Imm), ncol = 1), Imm))
    }

    chi2 <- c(0, 7.81, 13.9, 25, 42)
    run <- 0

    # partition
    while (npar > 0) {
        run <- run + 1
        apoc <- poc[npar]
        akon <- kon[npar]
        apor <- poradi[seq(apoc, akon)]
        Nex <- length(apor)
        margave <- floor((marg[npar, seq_len(cdim)] + marg[npar, seq(cdim + 1, dim2)])/2)
        J <- (ydat[apor, ] <= (matrix(1, nrow = Nex, ncol = 1) %*% margave)) * 1

        cI <- matrix(0, nrow = Nex, ncol = ddim)
        amarg <- matrix(1, nrow = ddim, ncol = 1) %*%
            subset(marg, subset = seq_len(nrow(marg)) %in% npar)

        for (d in seq_len(ddim)) {
            cI[, d] <- matrix(1, nrow = Nex, ncol = 1)
            for (k in seq_len(cdim)) {
                if (Imm[d, k] == 1) {
                  cI[, d] <- 1 * (cI[, d] & (!J[, k]))
                  amarg[d, k] <- margave[k] + 1
                } else {
                  cI[, d] <- 1 * (cI[, d] & J[, k])
                  amarg[d, k + cdim] <- margave[k]
                }
            }
        }

        NN <- colSums(cI)
        tst <- ddim * sum((NN - Nex/ddim * matrix(1, nrow = 1, ncol = ddim))^2)/Nex

        if ((tst > chi2[cdim]) | (run == 1)) {
            # decide partition or not
            npar <- npar - 1
            for (ind in seq_len(ddim)) {
                if (NN[ind] > ddim) {
                  npar <- npar + 1
                  akon <- apoc + NN[ind] - 1
                  poc[npar] <- apoc
                  kon[npar] <- akon
                  marg[npar, ] <- amarg[ind, ]
                  poradi[seq(apoc, akon)] <- apor[which(cI[, ind] != 0, arr.ind = TRUE)]
                  apoc <- akon + 1
                } else {
                  if (NN[ind] > 0) {
                    Nxx <- apply(amarg[ind, seq(cdim + 1, dim2)] - amarg[ind, seq_len(cdim)] + matrix(1,
                      nrow = 1, ncol = cdim), 2, prod)

                    Nz <- amarg[ind, 6] - amarg[ind, 3] + 1
                    Jx <- 1 * ((ydat[, 1] >= amarg[ind, 1]) & (ydat[, 1] <= amarg[ind, 4]))
                    Jy <- 1 * ((ydat[, 2] >= amarg[ind, 2]) & (ydat[, 2] <= amarg[ind, 5]))
                    Jz <- 1 * ((ydat[, 3] >= amarg[ind, 3]) & (ydat[, 3] <= amarg[ind, 6]))

                    Nxz <- sum(1 * (Jx & Jz))
                    Nyz <- sum(1 * (Jy & Jz))
                    cond <- (NN[ind] * Nz)/(Nxz * Nyz)
                    if (is.infinite(cond))
                      cond <- 1
                    if (cond == 0)
                      cond <- 1
                    cmi <- cmi + NN[ind] * log(cond)
                  }
                }
            }
        } else {

            Nxx <- apply(marg[npar, seq(cdim + 1, dim2)] - marg[npar, seq_len(cdim)] + matrix(1,
                nrow = 1, ncol = cdim), 2, prod)

            Nz <- marg[npar, 6] - marg[npar, 3] + 1
            Jx <- 1 * ((ydat[, 1] >= marg[npar, 1]) & (ydat[, 1] <= marg[npar, 4]))
            Jy <- 1 * ((ydat[, 2] >= marg[npar, 2]) & (ydat[, 2] <= marg[npar, 5]))
            Jz <- 1 * ((ydat[, 3] >= marg[npar, 3]) & (ydat[, 3] <= marg[npar, 6]))

            Nxz <- sum(1 * (Jx & Jz))
            Nyz <- sum(1 * (Jy & Jz))
            cond <- (Nex * Nz)/(Nxz * Nyz)
            if (is.infinite(cond))
                cond <- 1
            if (cond == 0)
                cond <- 1
            cmi <- cmi + Nex * log(cond)
            npar <- npar - 1
        }
    }
    # normalize
    cmi <- cmi/N
    return(cmi)
}

## combpvalue: Combine p-values using Fisher's method
combpvalue <- function(p_values) {

    # calculate chi-square statistic and combined p-value
    Q <- -2 * sum(log(p_values))
    Degree <- length(p_values)
    Combp <- pchisq(Q, 2 * Degree, lower.tail = FALSE)
    return(Combp)
}

## predCor: Predict competing endogenous RNAs via evidence for competition for miRNA regulation
## based on conditional mutual information (CMI) or partial pearson correlation (PPC)
predCor <- function(expr, num_perm, method = c("CMI", "PPC")) {

    exprRownames <- rownames(expr)
    expr <- as.matrix(expr)

    # generate random sequence
    num_sample <- ncol(expr)
    perm_seq <- t(sapply(seq_len(num_perm), function(i) sample(num_sample)))

    # identify mediator candidate
    num_cand <- nrow(expr) - 2

    # evaluate significance of triplet
    tri_id <- matrix(0, nrow = 2 * num_cand, ncol = 3)
    tri_cor <- matrix(0, nrow = 2 * num_cand, ncol = 1)
    tri_pval <- rep(1, 2 * num_cand)

    for (i in seq_len(2)) {
        for (j in seq_len(num_cand)) {

            idx_tar <- i  # target index
            idx_miR <- j + 2  # miRNA index
            idx_mod <- 2/i  # modulator index
            idx_tri <- (i - 1) * num_cand + j  # triplet index

            tri_id[idx_tri, ] <- c(exprRownames[idx_tar], exprRownames[idx_miR], exprRownames[idx_mod])

            # calculate Cor( target ; miRNA | modulator )
            data <- t(expr[c(idx_tar, idx_miR, idx_mod), ])

            # construct null distribution
            rand_exp <- lapply(seq_len(num_perm), function(i) expr[idx_mod, ][perm_seq[i, ]])
            nulldata <- lapply(seq_len(num_perm), function(i) cbind(t(expr[c(idx_tar, idx_miR),
                ]), rand_exp[[i]]))

            if (method == "CMI") {
                tri_cor[idx_tri, ] <- calCMI(data)
                null <- unlist(lapply(seq_len(num_perm), function(i) calCMI(nulldata[[i]])))
            } else if (method == "PPC") {
                tri_cor[idx_tri, ] <- corpcor::pcor.shrink(data, verbose = FALSE)[1, 2]
                null <- unlist(lapply(seq_len(num_perm), function(i) corpcor::pcor.shrink(nulldata[[i]],
                  verbose = FALSE)[1, 2]))
            }

            # calculate p-value
            tri_pval[idx_tri] <- max(1, sum(tri_cor[idx_tri, ] <= null))/num_perm
        }
    }

    # evaluate significance of interaction
    tri_idx <- order(tri_pval)
    tri_pval <- tri_pval[tri_idx]

    tri_id <- tri_id[tri_idx, ]
    tri_cor <- tri_cor[tri_idx]
    pcomb <- as.vector(sapply(seq_along(tri_pval), function(i) combpvalue(tri_pval[seq_len(i)])[1]))

    # identify final mediators
    min_pval <- min(pcomb)
    return(min_pval)
}

## Internal functions (parMM, graphWeights, recommendation, dtHybrid) of cernia method are from
## the website: https://github.com/dsardina/cernia Copyright 2016 Rosalba Giugno Licensed under
## the Apache License, Version 2.0 (the 'License'); you may not use this file except in
## compliance with the License. You may obtain a copy of the License at
## http://www.apache.org/licenses/LICENSE-2.0

## Unless required by applicable law or agreed to in writing, software distributed under the
## License is distributed on an 'AS IS' BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
## either express or implied. See the License for the specific language governing permissions and
## limitations under the License.

## Matrix multiplication in parallel
parMM <- function(cl, A, B) {

    if (!all(is.na(cl)) && is.object(cl)) {

        nA <- nrow(A)
        ncl <- length(cl)
        # Split an indice equally
        i <- seq_len(nA)
        if (ncl == 1) {
            splitIndices <- i
        } else {
            fuzz <- min((nA - 1)/1000, 0.4 * nA/ncl)
            breaks <- seq(1 - fuzz, nA + fuzz, length = ncl + 1)
            splitIndices <- structure(split(i, cut(i, breaks)), names = NULL)
        }

        # Split a matrix equally according to the rows of the matrix
        splitRows <- lapply(splitIndices, function(i) A[i, , drop = FALSE])

        fun <- get(as.character("rbind"), envir = .GlobalEnv, mode = "function")
        R <- do.call("fun", lapply(clusterApply(cl = cl, x = splitRows, get("%*%"), B), enquote))
    } else {
        R <- A %*% B
    }
    return(R)
}

## The first step of DT-Hybrid recommendation algorithm: generating the weights for each pair of
## target nodes
graphWeights <- function(n, m, A, lambda = 0.5, alpha = 0.5, S = NA, S1 = NA, cl = NA) {

    if (nrow(A) != n || ncol(A) != m) {
        stop("The matrix A should be an n by m matrix.")
    }

    has.similarity <- (!all(is.na(S)) && is.matrix(S) && !all(is.na(S1)) && is.matrix(S1))

    if (has.similarity) {
        if (nrow(S1) != m || ncol(S1) != m) {
            stop("The matrix S1 should be an m by m matrix.")
        }
        if (nrow(S) != n || ncol(S) != n) {
            stop("The matrix S should be an n by n matrix.")
        }
    }

    Ky <- diag(1/colSums(A))
    Ky[is.infinite(Ky) | is.na(Ky)] <- 0  #BugFix: 1/0=Infinite replaced with 0

    kx <- rowSums(A)
    Nx <- 1/(matrix(kx, nrow = n, ncol = n, byrow = TRUE)^(lambda) * matrix(kx, nrow = n, ncol = n,
        byrow = FALSE)^(1 - lambda))

    Nx[is.infinite(Nx) | is.na(Nx)] <- 0  #BugFix: 1/0=Infinite replaced with 0
    kx[is.infinite(kx) | is.na(kx)] <- 0  #BugFix: 1/0=Infinite replaced with 0

    W <- t(parMM(cl, A, Ky))
    W <- parMM(cl, A, W)
    W <- Nx * W
    rownames(W) <- rownames(A)
    colnames(W) <- rownames(A)

    if (has.similarity) {
        X5 <- parMM(cl, A, S1)
        X6 <- parMM(cl, X5, t(A))
        X7 <- parMM(cl, A, matrix(1, nrow = m, ncol = m))
        X8 <- parMM(cl, X7, t(A))
        S2 <- X6/X8
        W <- W * (1 + (alpha * S) + ((1 - alpha) * S2))
    }

    W[is.nan(W)] <- 0  #This should never happen
    return(W)
}

## The second step of DT-Hybrid recommendation algorithm: generating ecommendation scores of each
## RNA-RNA pair
recommendation <- function(A, lambda = 0.5, alpha = 0.5, S = NA, S1 = NA, cl = NA) {

    n <- nrow(A)
    m <- ncol(A)
    W <- graphWeights(n = n, m = m, A = A, lambda = lambda, alpha = alpha, S = S, S1 = S1, cl = cl)

    R <- parMM(cl, W, A)
    return(R)
}

## Make projection from bipartite network using DT-hybrid sources
dtHybrid <- function(miRTarget) {

    # Extract miRs and their targets
    mir <- unique(miRTarget[, 1])
    tar <- unique(miRTarget[, 2])

    # Create the matrix of the miRTarget
    A <- matrix(nrow = length(tar), ncol = length(mir), data = 0)
    colnames(A) <- mir
    rownames(A) <- tar

    for (i in seq_len(nrow(miRTarget))) {
        A[which(tar %in% as.character(miRTarget[i, 2])),
            which(mir %in% as.character(miRTarget[i, 1]))] <- 1
    }

    # Make projection from bipartite network using DT-hybrid sources
    cl <- makeCluster(detectCores() - 2)
    M <- recommendation(A, cl = cl)
    W <- graphWeights(nrow(M), ncol(M), M, cl = cl)

    stopCluster(cl)

    return(W)
}

## Internal function cluster from ProNet package 
## (https://github.com/cran/ProNet) with GPL-2 license.
cluster <- function(graph, method="MCL", expansion = 2, inflation = 2, 
                  hcmethod = "average", directed = FALSE, outfile = NULL, ...)
{
  method<-match.arg(method)
  if(method=="FN"){
	  graph <- simplify(graph)
	  fc <- fastgreedy.community(graph, merges = TRUE, modularity = TRUE)
	  membership <- membership(fc)
	  if(!is.null(V(graph)$name)){
              names(membership) <- V(graph)$name
	  }	  
	  if(!is.null(outfile)){
	      cluster.save(cbind(names(membership),membership),outfile=outfile)		
	  }else{
	      return(membership)
	  }
  }else if(method=="LINKCOMM"){
	  edgelist <- get.edgelist(graph)
	  if(!is.null(E(graph)$weight)){
              edgelist <- cbind(edgelist,E(graph)$weight)
	  }
	  lc <- getLinkCommunities(edgelist,plot=FALSE,directed=directed,hcmethod=hcmethod)
	  if(!is.null(outfile)){
		  cluster.save(lc$nodeclusters,outfile=outfile)		  
	  }else{		  
		  return(lc$nodeclusters)
	  }
  }else if(method=="MCL"){
        adj <- matrix(rep(0,length(V(graph))^2),nrow=length(V(graph)),ncol=length(V(graph)))
        for(i in seq_along(V(graph))){
            neighbors <- neighbors(graph,v=V(graph)$name[i],mode="all")
            j <- match(neighbors$name,V(graph)$name,nomatch=0)
            adj[i,j] = 1
        }
        lc <- mcl(adj,addLoops=TRUE,expansion=expansion,inflation=inflation,allow1=TRUE,max.iter=100,ESM=FALSE)
        lc$name <- V(graph)$name
        lc$Cluster <- lc$Cluster
        
        if(!is.null(outfile)){
            cluster.save(cbind(lc$name,lc$Cluster),outfile=outfile)		
        }else{
        result <- lc$Cluster
        names(result) <- V(graph)$name
        return(result)
        }
  }else if(method=="MCODE"){
	  compx <- mcode(graph,vwp=0.9,haircut=T,fluff=T,fdt=0.1)
	  index <- which(!is.na(compx$score))
	  membership <- rep(0,vcount(graph))
	  for(i in seq_along(index)){
	      membership[compx$COMPLEX[[index[i]]]]<-i
	  }
	      if(!is.null(V(graph)$name)) names(membership)<-V(graph)$name	  
	  if(!is.null(outfile)){
		  cluster.save(cbind(names(membership),membership),outfile=outfile)
		  invisible(NULL)
	  }else{
		  return(membership)
	  }
  }
}

## Internal function cluster.save from ProNet package
## (https://github.com/cran/ProNet) with GPL-2 license.
cluster.save <- function(membership, outfile){
	wd <- dirname(outfile)
	wd <- ifelse(wd==".",paste(wd,"/",sep=""),wd)
	filename <- basename(outfile)
	if((filename=="")||(grepl(":",filename))){
		filename <- "membership.txt"
	}else if(grepl("\\.",filename)){
		filename <- sub("\\.(?:.*)",".txt", filename)
	}
	write.table(membership,file=paste(wd,filename,sep="/"),
              row.names=FALSE,col.names=c("node","cluster"),quote =FALSE)
}

## Internal function mcode.vertex.weighting from ProNet package
## (https://github.com/cran/ProNet) with GPL-2 license.
mcode.vertex.weighting<-function(graph, neighbors){	
	stopifnot(is.igraph(graph))  
	weight <- lapply(seq_len(vcount(graph)),
                 function(i){
		              subg<-induced.subgraph(graph,neighbors[[i]])
		              core<-graph.coreness(subg)
		              k<-max(core)
				          ### k-coreness
				          kcore<-induced.subgraph(subg,which(core==k))
				          if(vcount(kcore)>1){
					          if(any(is.loop(kcore))){
						          k*ecount(kcore)/choose(vcount(kcore)+1,2)						
					          }else{
						          k*ecount(kcore)/choose(vcount(kcore),2)
					          }
				          }else{
                                             0
				          }
				         }
			 )
  
	return(unlist(weight))
}

## Internal function mcode.find.complex from ProNet package
## (https://github.com/cran/ProNet) with GPL-2 license.
mcode.find.complex <- function(neighbors, neighbors.indx, vertex.weight,
                             vwp, seed.vertex, seen)
{
	  res<-.C("complex",as.integer(neighbors),as.integer(neighbors.indx),
          as.single(vertex.weight),as.single(vwp),as.integer(seed.vertex),
          seen=as.integer(seen),COMPLEX=as.integer(rep(0,length(seen))), PACKAGE = "miRsponge"
          )
	
	  return(list(seen=res$seen,COMPLEX=which(res$COMPLEX!=0)))
}

## Internal function mcode.find.complexex from ProNet package
## (https://github.com/cran/ProNet) with GPL-2 license.
mcode.find.complexex <- function(graph, neighbors, vertex.weight, vwp)
{
	seen<-rep(0,vcount(graph))

	neighbors<-lapply(neighbors,function(item){item[-1]})
	neighbors.indx<-cumsum(unlist(lapply(neighbors,length)))
	
	neighbors.indx<-c(0,neighbors.indx)
	neighbors<-unlist(neighbors)-1
	
	COMPLEX<-list()
	n<-1
        w.order<-order(vertex.weight,decreasing=TRUE)
	for(i in w.order){
		if(!(seen[i])){
			res<-mcode.find.complex(neighbors,neighbors.indx,vertex.weight,vwp,i-1,seen)
			if(length(res$COMPLEX)>1){
				COMPLEX[[n]]<-res$COMPLEX
				seen<-res$seen
				n<-n+1
			}
		}
	}	
	rm(neighbors)	
	return(list(COMPLEX=COMPLEX,seen=seen))
}

## Internal function mcode.fluff.complex from ProNet package
## (https://github.com/cran/ProNet) with GPL-2 license.
mcode.fluff.complex <- function(graph, vertex.weight, fdt=0.8, complex.g, seen)
{
	seq_complex.g<-seq_along(complex.g)
	for(i in seq_complex.g){
	    node.neighbor<-unlist(neighborhood(graph,1,complex.g[i]))
	    if(length(node.neighbor)>1){
                subg<-induced.subgraph(graph,node.neighbor)
            if(graph.density(subg, loops=FALSE)>fdt){
                complex.g<-c(complex.g,node.neighbor)
       }
		 }
	}
  
	return(unique(complex.g))
}

## Internal function mcode.post.process from ProNet package
## (https://github.com/cran/ProNet) with GPL-2 license.
mcode.post.process<-function(graph, vertex.weight, haircut, fluff, fdt=0.8,
                             set.complex.g, seen)
{
	indx<-unlist(lapply(set.complex.g,
                      function(complex.g){
		          if(length(complex.g)<=2)
			      0
			  else
		              1
			  }
		    ))
	set.complex.g<-set.complex.g[indx!=0]
	set.complex.g<-lapply(set.complex.g,
                        function(complex.g){
			    coreness<-graph.coreness(induced.subgraph(graph,complex.g))						
			    if(fluff){
				complex.g<-mcode.fluff.complex(graph,vertex.weight,fdt,complex.g,seen)
			    if(haircut){
				## coreness needs to be recalculated
				coreness<-graph.coreness(induced.subgraph(graph,complex.g))
				complex.g<-complex.g[coreness>1]
				}
				}else if(haircut){
				complex.g<-complex.g[coreness>1]
				}
				return(complex.g)
				})
	set.complex.g<-set.complex.g[lapply(set.complex.g,length)>2]
	return(set.complex.g)
}

## Internal function mcode from ProNet package
## (https://github.com/cran/ProNet) with GPL-2 license.
mcode <- function(graph, vwp=0.5, haircut=FALSE, fluff=FALSE, fdt=0.8, loops=TRUE)
{
	stopifnot(is.igraph(graph))
	if(vwp>1 | vwp <0){
            stop("vwp must be between 0 and 1")
	}
	if(!loops){
            graph<-simplify(graph,remove.multiple=FALSE,remove.loops=TRUE)
	}
	neighbors<-neighborhood(graph,1)
	W<-mcode.vertex.weighting(graph,neighbors)
	res<-mcode.find.complexex(graph,neighbors=neighbors,vertex.weight=W,vwp=vwp)
	COMPLEX<-mcode.post.process(graph,vertex.weight=W,haircut=haircut,fluff=fluff,
                              fdt=fdt,res$COMPLEX,res$seen)		
	score<-unlist(lapply(COMPLEX,
                       function(complex.g){
		           complex.g<-induced.subgraph(graph,complex.g)
			   if(any(is.loop(complex.g)))
			   score<-ecount(complex.g)/choose(vcount(complex.g)+1,2)*vcount(complex.g)
			   else
			   score<-ecount(complex.g)/choose(vcount(complex.g),2)*vcount(complex.g)
			   return(score)
			}
		    ))
	order_score<-order(score,decreasing=TRUE)
	return(list(COMPLEX=COMPLEX[order_score],score=score[order_score]))
}

## Utility methods for identifying miRNA sponge interactions For input expression data, the
## columns are genes and the rows are samples.  For input miRTarget, the miRNA-target
## interactions could be miRNA-mRNA, miRNA-lncRNA, miRNA-circRNA, miRNA-pseudogene, etc.  For
## input mres, each row contains five elements: Mirna, Target, energy, gap_l, gap_r.

## 1. miRHomology
miRHomology <- function(miRTarget, minSharedmiR = 3, padjustvaluecutoff = 0.01, padjustmethod = "BH") {

    miRTarget <- as.matrix(miRTarget)

    m1 <- nrow(miRTarget)
    n1 <- ncol(miRTarget)

    miR <- miRTarget[, 1]
    tar <- miRTarget[, 2]

    miRSym <- unique(miR)
    targetSym <- unique(tar)

    m2 <- length(targetSym)

    # Initialize variables
    ceRInt <- matrix(NA, m2 * (m2 - 1)/2, 2)
    C <- matrix(NA, m2 * (m2 - 1)/2, 2)

    for (i in seq_len(m2 - 1)) {
        for (j in seq(i + 1, m2)) {

            Interin1 <- miRTarget[which(miRTarget[, 2] %in% targetSym[i]), 1]
            Interin2 <- miRTarget[which(miRTarget[, 2] %in% targetSym[j]), 1]

            M1 <- length(Interin1)
            M2 <- length(Interin2)
            M3 <- length(intersect(Interin1, Interin2))
            M4 <- length(miRSym)
            M5 <- 1 - phyper(M3 - 1, M2, M4 - M2, M1)

            if (M3 >= minSharedmiR & M5 < padjustvaluecutoff) {

                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- targetSym[i]
                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- targetSym[j]

                C[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- M3
                C[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- M5

            }

        }
    }

    # Extract RNA-RNA pair based on the homology of miRNA binding sites present in the 3'UTR
    # sequence
    ceRInt <- ceRInt[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff) == "TRUE"), ]

    C <- C[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff) == "TRUE"), ]

    if (is.vector(C)) {
        miRHomologyceRInt <- c(ceRInt, C)
        names(miRHomologyceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs")
    } else {
        miRHomologyceRInt <- cbind(ceRInt, C)
        colnames(miRHomologyceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs")
    }

    return(miRHomologyceRInt)
}

## 2. Positive Correlation (PC) method
pc <- function(miRTarget, ExpData, minSharedmiR = 3, poscorcutoff = 0, padjustvaluecutoff = 0.01,
    padjustmethod = "BH") {

    miRTargetCandidate <- querymiRTargetbinding(ExpData, miRTarget)
    miRTargetCandidate <- as.matrix(miRTargetCandidate)
    m1 <- nrow(miRTargetCandidate)
    n1 <- ncol(miRTargetCandidate)

    ExpDataNames <- c(as.matrix(ExpData[1, ]))
    ExpData <- unfactor(ExpData[-1, ])    
    colnames(ExpData) <- ExpDataNames

    miR <- miRTargetCandidate[, 1]
    tar <- miRTargetCandidate[, 2]

    miRSym <- unique(miR)
    targetSym <- unique(tar)

    m2 <- length(targetSym)

    # Initialize variables
    ceRInt <- matrix(NA, m2 * (m2 - 1)/2, 2)
    C <- matrix(NA, m2 * (m2 - 1)/2, 4)

    for (i in seq_len(m2 - 1)) {
        for (j in seq(i + 1, m2)) {

            Interin1 <- miRTargetCandidate[which(miRTargetCandidate[, 2] %in% targetSym[i]), 1]
            Interin2 <- miRTargetCandidate[which(miRTargetCandidate[, 2] %in% targetSym[j]), 1]

            M1 <- length(Interin1)
            M2 <- length(Interin2)
            M3 <- length(intersect(Interin1, Interin2))
            M4 <- length(miRSym)
            M5 <- 1 - phyper(M3 - 1, M2, M4 - M2, M1)

            if (M3 >= minSharedmiR & M5 < padjustvaluecutoff) {

                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- targetSym[i]
                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- targetSym[j]

                tarExpIdx1 <- which(ExpDataNames %in% targetSym[i])
                tarExpIdx2 <- which(ExpDataNames %in% targetSym[j])

                # Calculate Pearson correlation of each RNA-RNA pair
                M6 <- cor.test(ExpData[, tarExpIdx1], ExpData[, tarExpIdx2])$estimate
                M7 <- cor.test(ExpData[, tarExpIdx1], ExpData[, tarExpIdx2])$p.value

                C[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- M3
                C[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- M5
                C[(i - 1) * m2 + j - sum(seq_len(i)), 3] <- M6
                C[(i - 1) * m2 + j - sum(seq_len(i)), 4] <- M7
            }
        }
    }

    # Extract positive correlated RNA-RNA pairs.
    ceRInt <- ceRInt[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & C[,
        3] > poscorcutoff & p.adjust(C[, 4], method = padjustmethod) < padjustvaluecutoff) == "TRUE"),
        ]

    C <- C[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & C[, 3] > poscorcutoff &
        p.adjust(C[, 4], method = padjustmethod) < padjustvaluecutoff) == "TRUE"), ]

    if (is.vector(C)) {
        PCceRInt <- c(ceRInt, C)
        names(PCceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "correlation", "p.adjusted_value of positive correlation")
    } else {
        PCceRInt <- cbind(ceRInt, C)
        colnames(PCceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "correlation", "p.adjusted_value of positive correlation")
    }

    return(PCceRInt)
}

## 3. Sensitivity Partial Pearson Correlation (SPPC) method
sppc <- function(miRTarget, ExpData, minSharedmiR = 3, poscorcutoff = 0, padjustvaluecutoff = 0.01,
    padjustmethod = "BH", senscorcutoff = 0.3) {

    miRTargetCandidate <- querymiRTargetbinding(ExpData, miRTarget)
    miRTargetCandidate <- as.matrix(miRTargetCandidate)
    m1 <- nrow(miRTargetCandidate)
    n1 <- ncol(miRTargetCandidate)

    ExpDataNames <- c(as.matrix(ExpData[1, ]))
    ExpData <- unfactor(ExpData[-1, ])
    colnames(ExpData) <- ExpDataNames

    miR <- miRTargetCandidate[, 1]
    tar <- miRTargetCandidate[, 2]

    miRSym <- unique(miR)
    targetSym <- unique(tar)

    m2 <- length(targetSym)

    # Initialize variables
    ceRInt <- matrix(NA, m2 * (m2 - 1)/2, 2)
    C <- matrix(NA, m2 * (m2 - 1)/2, 5)

    for (i in seq_len(m2 - 1)) {
        for (j in seq(i + 1, m2)) {

            Interin1 <- miRTargetCandidate[which(miRTargetCandidate[, 2] %in% targetSym[i]), 1]
            Interin2 <- miRTargetCandidate[which(miRTargetCandidate[, 2] %in% targetSym[j]), 1]

            M1 <- length(Interin1)
            M2 <- length(Interin2)
            M3 <- length(intersect(Interin1, Interin2))
            M4 <- length(miRSym)
            M5 <- 1 - phyper(M3 - 1, M2, M4 - M2, M1)

            if (M3 >= minSharedmiR & M5 < padjustvaluecutoff) {

                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- targetSym[i]
                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- targetSym[j]

                tarExpIdx1 <- which(ExpDataNames %in% targetSym[i])
                tarExpIdx2 <- which(ExpDataNames %in% targetSym[j])
                miRExpIdx <- which(ExpDataNames %in% intersect(Interin1, Interin2))

                # Calculate sensitivity partial pearson correlation of each RNA-RNA pair
                M6 <- cor.test(ExpData[, tarExpIdx1], ExpData[, tarExpIdx2])$estimate
                M7 <- cor.test(ExpData[, tarExpIdx1], ExpData[, tarExpIdx2])$p.value
                M8 <- M6 - corpcor::pcor.shrink(cbind(ExpData[, tarExpIdx1], ExpData[, tarExpIdx2],
                  ExpData[, miRExpIdx]), verbose = FALSE)[1, 2]

                C[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- M3
                C[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- M5
                C[(i - 1) * m2 + j - sum(seq_len(i)), 3] <- M6
                C[(i - 1) * m2 + j - sum(seq_len(i)), 4] <- M7
                C[(i - 1) * m2 + j - sum(seq_len(i)), 5] <- M8
            }
        }
    }

    # Extract RNA-RNA pairs with sensitivity correlation more than senscorcutoff.
    ceRInt <- ceRInt[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & C[,
        3] > poscorcutoff & p.adjust(C[, 4], method = padjustmethod) < padjustvaluecutoff & C[,
        5] > senscorcutoff) == "TRUE"), ]

    C <- C[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & C[, 3] > poscorcutoff &
        p.adjust(C[, 4], method = padjustmethod) < padjustvaluecutoff & C[, 5] > senscorcutoff) ==
        "TRUE"), ]

    if (is.vector(C)) {
        SPPCceRInt <- c(ceRInt, C)
        names(SPPCceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "correlation", "p.adjusted_value of positive correlation", "sensitivity partial pearson correlation")
    } else {
        SPPCceRInt <- cbind(ceRInt, C)
        colnames(SPPCceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "correlation", "p.adjusted_value of positive correlation", "sensitivity partial pearson correlation")
    }

    return(SPPCceRInt)
}

## 4. Partial Pearson Correlation (PPC) method
ppc <- function(miRTarget, ExpData, minSharedmiR = 3, num_perm = 100, padjustvaluecutoff = 0.01,
    padjustmethod = "BH") {

    miRTargetCandidate <- querymiRTargetbinding(ExpData, miRTarget)
    miRTargetCandidate <- as.matrix(miRTargetCandidate)
    m1 <- nrow(miRTargetCandidate)
    n1 <- ncol(miRTargetCandidate)

    ExpDataNames <- c(as.matrix(ExpData[1, ]))
    ExpData <- unfactor(ExpData[-1, ])
    colnames(ExpData) <- ExpDataNames

    miR <- miRTargetCandidate[, 1]
    tar <- miRTargetCandidate[, 2]

    miRSym <- unique(miR)
    targetSym <- unique(tar)

    m2 <- length(targetSym)

    # Initialize variables
    ceRInt <- matrix(NA, m2 * (m2 - 1)/2, 2)
    C <- matrix(NA, m2 * (m2 - 1)/2, 3)

    for (i in seq_len(m2 - 1)) {
        for (j in seq(i + 1, m2)) {

            Interin1 <- miRTargetCandidate[which(miRTargetCandidate[, 2] %in% targetSym[i]), 1]
            Interin2 <- miRTargetCandidate[which(miRTargetCandidate[, 2] %in% targetSym[j]), 1]
            miRExpIdx <- which(ExpDataNames %in% intersect(Interin1, Interin2))

            M1 <- length(Interin1)
            M2 <- length(Interin2)
            M3 <- length(intersect(Interin1, Interin2))
            M4 <- length(miRSym)
            M5 <- 1 - phyper(M3 - 1, M2, M4 - M2, M1)

            if (M3 >= minSharedmiR & M5 < padjustvaluecutoff) {

                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- targetSym[i]
                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- targetSym[j]

                tarExpIdx1 <- which(ExpDataNames %in% targetSym[i])
                tarExpIdx2 <- which(ExpDataNames %in% targetSym[j])

                inputdata <- t(cbind(ExpData[, tarExpIdx1], ExpData[, tarExpIdx2], ExpData[, miRExpIdx]))

                C[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- M3
                C[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- M5
                C[(i - 1) * m2 + j - sum(seq_len(i)), 3] <- predCor(inputdata, num_perm, method = "PPC")

            }
        }
    }

    # Extract significant RNA-RNA pairs.
    ceRInt <- ceRInt[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & p.adjust(C[,
        3], method = padjustmethod) < padjustvaluecutoff) == "TRUE"), ]

    C <- C[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & p.adjust(C[, 3],
        method = padjustmethod) < padjustvaluecutoff) == "TRUE"), ]

    if (is.vector(C)) {
        PPCceRInt <- c(ceRInt, C)
        names(PPCceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "p.adjusted_value of RNA competition")

    } else {
        PPCceRInt <- cbind(ceRInt, C)
        colnames(PPCceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "p.adjusted_value of RNA competition")
    }

    return(PPCceRInt)
}

## 5. Hermes method
hermes <- function(miRTarget, ExpData, minSharedmiR = 3, num_perm = 100, padjustvaluecutoff = 0.01,
    padjustmethod = "BH") {

    miRTargetCandidate <- querymiRTargetbinding(ExpData, miRTarget)
    miRTargetCandidate <- as.matrix(miRTargetCandidate)
    m1 <- nrow(miRTargetCandidate)
    n1 <- ncol(miRTargetCandidate)

    ExpDataNames <- c(as.matrix(ExpData[1, ]))
    ExpData <- unfactor(ExpData[-1, ])    
    colnames(ExpData) <- ExpDataNames

    miR <- miRTargetCandidate[, 1]
    tar <- miRTargetCandidate[, 2]

    miRSym <- unique(miR)
    targetSym <- unique(tar)

    m2 <- length(targetSym)

    # Initialize variables
    ceRInt <- matrix(NA, m2 * (m2 - 1)/2, 2)
    C <- matrix(NA, m2 * (m2 - 1)/2, 3)

    for (i in seq_len(m2 - 1)) {
        for (j in seq(i + 1, m2)) {

            Interin1 <- miRTargetCandidate[which(miRTargetCandidate[, 2] %in% targetSym[i]), 1]
            Interin2 <- miRTargetCandidate[which(miRTargetCandidate[, 2] %in% targetSym[j]), 1]
            miRExpIdx <- which(ExpDataNames %in% intersect(Interin1, Interin2))

            M1 <- length(Interin1)
            M2 <- length(Interin2)
            M3 <- length(intersect(Interin1, Interin2))
            M4 <- length(miRSym)
            M5 <- 1 - phyper(M3 - 1, M2, M4 - M2, M1)

            if (M3 >= minSharedmiR & M5 < padjustvaluecutoff) {

                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- targetSym[i]
                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- targetSym[j]

                tarExpIdx1 <- which(ExpDataNames %in% targetSym[i])
                tarExpIdx2 <- which(ExpDataNames %in% targetSym[j])

                inputdata <- t(cbind(ExpData[, tarExpIdx1], ExpData[, tarExpIdx2], ExpData[, miRExpIdx]))

                C[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- M3
                C[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- M5
                C[(i - 1) * m2 + j - sum(seq_len(i)), 3] <- predCor(inputdata, num_perm, method = "CMI")
            }
        }
    }

    # Extract significant RNA-RNA pairs.
    ceRInt <- ceRInt[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & p.adjust(C[,
        3], method = padjustmethod) < padjustvaluecutoff) == "TRUE"), ]

    C <- C[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & p.adjust(C[, 3],
        method = padjustmethod) < padjustvaluecutoff) == "TRUE"), ]

    if (is.vector(C)) {
        HermesceRInt <- c(ceRInt, C)
        names(HermesceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "p.adjusted_value of RNA competition")
    } else {
        HermesceRInt <- cbind(ceRInt, C)
        colnames(HermesceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "p.adjusted_value of RNA competition")
    }

    return(HermesceRInt)
}

## 6. MuTaME method
muTaME <- function(miRTarget, mres, minSharedmiR = 3, padjustvaluecutoff = 0.01, padjustmethod = "BH",
    scorecutoff = 0.5) {

    miRTarget <- as.matrix(miRTarget)
    m1 <- nrow(miRTarget)
    n1 <- ncol(miRTarget)

    miR <- miRTarget[, 1]
    tar <- miRTarget[, 2]

    miRSym <- unique(miR)
    targetSym <- unique(tar)

    m2 <- length(targetSym)

    # Initialize variables
    ceRInt <- matrix(NA, m2 * (m2 - 1)/2, 2)
    C <- matrix(NA, m2 * (m2 - 1)/2, 8)

    for (i in seq_len(m2 - 1)) {
        for (j in seq(i + 1, m2)) {

            Interin1 <- miRTarget[which(miRTarget[, 2] %in% targetSym[i]), 1]
            Interin2 <- miRTarget[which(miRTarget[, 2] %in% targetSym[j]), 1]
            cm <- intersect(Interin1, Interin2)
            SharedMREs <- mres[mres[, 2] %in% c(targetSym[i], targetSym[j]) & mres[, 1] %in% cm, ]

            M1 <- length(Interin1)
            M2 <- length(Interin2)
            M3 <- length(intersect(Interin1, Interin2))
            M4 <- length(miRSym)
            M5 <- 1 - phyper(M3 - 1, M2, M4 - M2, M1)

            if (M3 >= minSharedmiR & M5 < padjustvaluecutoff & nrow(SharedMREs) > 0) {

                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- targetSym[i]
                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- targetSym[j]

                C[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- M3
                C[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- M5

                # Score 1 for the fraction of coomon miRNAs
                C[(i - 1) * m2 + j - sum(seq_len(i)), 3] <- log(M3/min(M1, M2))

                # Score 2 for the density of the MREs for all shared miRNAs
                C[(i - 1) * m2 + j - sum(seq_len(i)), 4] <- sum(sapply(cm, function(miR) {
                  MREs <- SharedMREs[SharedMREs[, 1] == miR, ]
                  if (nrow(MREs) <= 0) return(1)
                  MREslr <- MREs[, c("gap_l", "gap_r")]
                  D <- abs(max(MREslr[, 1]) - min(MREslr[, 2]))
                  return(log(nrow(MREs)/D))
                }))

                # Score 3 for the distribution of MREs of the putative RNA-RNA pairs
                C[(i - 1) * m2 + j - sum(seq_len(i)), 5] <- sum(sapply(cm, function(miR) {
                  positions <- SharedMREs[SharedMREs[, 1] == miR, c("gap_l", "gap_r")]
                  if (nrow(positions) <= 0) return(1)
                  return(log(abs(max(positions[, 1]) - min(positions[, 2]))^2/sum((positions[, 2] -
                    positions[, 1])^2)))
                }))

                # Score 4 for the relation between the overall number of MREs for a putative miRNA sponge,
                # compared with the number of miRNAs that yield these MREs
                B <- nrow(SharedMREs)
                if (B == length(unique(SharedMREs[, 1]))) {
                  C[(i - 1) * m2 + j - sum(seq_len(i)), 6] <- log(1/B)
                } else {
                  C[(i - 1) * m2 + j - sum(seq_len(i)), 6] <- log((B - length(unique(SharedMREs[,
                    1])) - 1)/B)
                }

                C[(i - 1) * m2 + j - sum(seq_len(i)), 7] <- C[(i - 1) * m2 + j - sum(seq_len(i)),
                  3] + C[(i - 1) * m2 + j - sum(seq_len(i)), 4] + C[(i - 1) * m2 + j - sum(seq_len(i)),
                  5] + C[(i - 1) * m2 + j - sum(seq_len(i)), 6]
            }
        }
    }

    # Extract RNA-RNA pair based on four scores.
    ceRInt <- ceRInt[apply(ceRInt, 1, function(x) !all(is.na(x))), ]
    C <- C[apply(C, 1, function(x) !all(is.na(x))), ]
    C[, 8] <- (C[, 7] - min(C[, 7]))/(max(C[, 7]) - min(C[, 7]))

    ceRInt <- ceRInt[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & C[,
        8] > scorecutoff) == "TRUE"), ]

    C <- C[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & C[, 8] > scorecutoff) ==
        "TRUE"), ]

    if (is.vector(C)) {
        MuTaMEceRInt <- c(ceRInt, C)
        names(MuTaMEceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "Score 1", "Score 2", "Score 3", "Score 4", "Combined score", "Normalized score")
    } else {
        MuTaMEceRInt <- cbind(ceRInt, C)
        colnames(MuTaMEceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "Score 1", "Score 2", "Score 3", "Score 4", "Combined score", "Normalized score")
    }
    return(MuTaMEceRInt)
}

## 7. CERNIA method
cernia <- function(miRTarget, ExpData, mres, minSharedmiR = 3, poscorcutoff = 0, padjustvaluecutoff = 0.01,
    padjustmethod = "BH", scorecutoff = 0.5) {

    miRTargetCandidate <- querymiRTargetbinding(ExpData, miRTarget)
    miRTargetCandidate <- as.matrix(miRTargetCandidate)
    m1 <- nrow(miRTargetCandidate)
    n1 <- ncol(miRTargetCandidate)

    ExpDataNames <- c(as.matrix(ExpData[1, ]))
    ExpData <- unfactor(ExpData[-1, ])    
    colnames(ExpData) <- ExpDataNames

    miR <- miRTargetCandidate[, 1]
    tar <- miRTargetCandidate[, 2]

    miRSym <- unique(miR)
    targetSym <- unique(tar)

    m2 <- length(targetSym)

    # Initialize variables
    ceRInt <- matrix(NA, m2 * (m2 - 1)/2, 2)
    C <- matrix(NA, m2 * (m2 - 1)/2, 11)

    for (i in seq_len(m2 - 1)) {
        for (j in seq(i + 1, m2)) {

            Interin1 <- miRTargetCandidate[which(miRTargetCandidate[, 2] %in% targetSym[i]), 1]
            Interin2 <- miRTargetCandidate[which(miRTargetCandidate[, 2] %in% targetSym[j]), 1]
            cm <- intersect(Interin1, Interin2)
            SharedMREs <- mres[mres[, 2] %in% c(targetSym[i], targetSym[j]) & mres[, 1] %in% cm, ]

            M1 <- length(Interin1)
            M2 <- length(Interin2)
            M3 <- length(intersect(Interin1, Interin2))
            M4 <- length(miRSym)
            M5 <- 1 - phyper(M3 - 1, M2, M4 - M2, M1)

            tarExpIdx1 <- which(ExpDataNames %in% targetSym[i])
            tarExpIdx2 <- which(ExpDataNames %in% targetSym[j])

            M6 <- cor.test(ExpData[, tarExpIdx1], ExpData[, tarExpIdx2])$estimate

            if (M3 >= minSharedmiR & M5 < padjustvaluecutoff & M6 > poscorcutoff &
                nrow(SharedMREs) > 0) {

                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- targetSym[i]
                ceRInt[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- targetSym[j]

                C[(i - 1) * m2 + j - sum(seq_len(i)), 1] <- M3
                C[(i - 1) * m2 + j - sum(seq_len(i)), 2] <- M5

                # Score 1 for the fraction of coomon miRNAs
                C[(i - 1) * m2 + j - sum(seq_len(i)), 3] <- log(M3/min(M1, M2))

                # Score 2 for the density of the MREs for all shared miRNAs
                C[(i - 1) * m2 + j - sum(seq_len(i)), 4] <- sum(sapply(cm, function(miR) {
                  MREs <- SharedMREs[SharedMREs[, 1] == miR, ]
                  if (nrow(MREs) <= 0) return(1)
                  MREslr <- MREs[, c("gap_l", "gap_r")]
                  D <- abs(max(MREslr[, 1]) - min(MREslr[, 2]))
                  return(log(nrow(MREs)/D))
                }))

                # Score 3 for the distribution of MREs of the putative RNA-RNA pairs
                C[(i - 1) * m2 + j - sum(seq_len(i)), 5] = sum(sapply(cm, function(miR) {
                  positions <- SharedMREs[SharedMREs[, 1] == miR, c("gap_l", "gap_r")]
                  if (nrow(positions) <= 0) return(1)
                  return(log(abs(max(positions[, 1]) - min(positions[, 2]))^2/sum((positions[, 2] -
                    positions[, 1])^2)))
                }))

                # Score 4 for the relation between the overall number of MREs for a putative miRNA sponge,
                # compared with the number of miRNAs that yield these MREs
                B <- nrow(SharedMREs)
                if (B == length(unique(SharedMREs[, 1]))) {
                  C[(i - 1) * m2 + j - sum(seq_len(i)), 6] <- log(1/B)
                } else {
                  C[(i - 1) * m2 + j - sum(seq_len(i)), 6] <- log((B - length(unique(SharedMREs[,
                    1])) - 1)/B)
                }

                # Score 5 for the density of the hybridization energies related to MREs for all shared miRNAs
                SharedMREs <- mres[mres[, 2] %in% c(targetSym[i], targetSym[j]) & mres[, 1] %in%
                  cm, ]
                C[(i - 1) * m2 + j - sum(seq_len(i)), 7] <- sum(sapply(cm, function(miR) {
                  MREs <- SharedMREs[SharedMREs[, 1] == miR, ]
                  if (nrow(MREs) <= 0) return(1)
                  MREslr <- MREs[, c("gap_l", "gap_r")]
                  D <- abs(max(MREslr[, 1]) - min(MREslr[, 2]))
                  return(log(sum(abs(MREs[, 3]))/D))
                }))

                # Score 6 for the DT-Hybrid recommendation scores
                cerna_recommendations <- dtHybrid(miRTargetCandidate)
                C[(i - 1) * m2 + j - sum(seq_len(i)), 8] <- cerna_recommendations[targetSym[i],
                  targetSym[j]]

                # Score 7 for the pairwise Peason correlation between putative RNA-RNA pair expression data
                C[(i - 1) * m2 + j - sum(seq_len(i)), 9] <- log(M6)

                C[(i - 1) * m2 + j - sum(seq_len(i)), 10] <- C[(i - 1) * m2 + j - sum(seq_len(i)),
                  3] + C[(i - 1) * m2 + j - sum(seq_len(i)), 4] + C[(i - 1) * m2 + j - sum(seq_len(i)),
                  5] + C[(i - 1) * m2 + j - sum(seq_len(i)), 6] + C[(i - 1) * m2 + j - sum(seq_len(i)),
                  7] + C[(i - 1) * m2 + j - sum(seq_len(i)), 8] + C[(i - 1) * m2 + j - sum(seq_len(i)),
                  9]
            }
        }
    }

    # Extract RNA-RNA pair based on seven scores
    ceRInt <- ceRInt[apply(ceRInt, 1, function(x) !all(is.na(x))), ]
    C <- C[apply(C, 1, function(x) !all(is.na(x))), ]
    C[, 11] <- (C[, 10] - min(C[, 10]))/(max(C[, 10]) - min(C[, 10]))

    ceRInt <- ceRInt[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & C[,
        11] > scorecutoff) == "TRUE"), ]

    C <- C[which((p.adjust(C[, 2], method = padjustmethod) < padjustvaluecutoff & C[, 11] > scorecutoff) ==
        "TRUE"), ]

    if (is.vector(C)) {
        CERNIAceRInt <- c(ceRInt, C)
        names(CERNIAceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "Score 1", "Score 2", "Score 3", "Score 4", "Score 5", "Score 6", "Score 7", "Combined score",
            "Normalized score")
    } else {
        CERNIAceRInt <- cbind(ceRInt, C)
        colnames(CERNIAceRInt) <- c("sponge_1", "sponge_2", "#shared miRNAs", "p.adjusted_value of shared miRNAs",
            "Score 1", "Score 2", "Score 3", "Score 4", "Score 5", "Score 6", "Score 7", "Combined score",
            "Normalized score")
    }

    return(CERNIAceRInt)

}

## 8. Integrate method for miRNA sponge interactions by integrating different methods.
integrateMethod <- function(Interlist, Intersect_num) {

    if (length(Interlist) >= 2 & length(Interlist) >= Intersect_num) {
        Combcase <- t(combn(length(Interlist), Intersect_num))
        Combnum <- dim(Combcase)[1]
        Integrate_Inter <- list()

        for (i in seq_len(Combnum)) {
            Interin <- do.call(rbind, lapply(Combcase[i, ], function(i) Interlist[[i]]))
            Interin_paste <- paste(Interin[, 1], Interin[, 2], sep = "-")
            Interin_table <- table(Interin_paste)
            Interin_names <- names(Interin_table)[which(Interin_table == Intersect_num)]
            Integrate_Inter[[i]] <- Interin[which(Interin_paste %in% Interin_names), ]
        }

        Integrate_res <- unique(do.call(rbind, Integrate_Inter))
        return(Integrate_res)
    } else {
        stop("Please check your input!\n")
    }

}

## Consolidating seven functions: miRHomology, pc, sppc, ppc, hermes, muTaME and cernia.
spongeMethod <- function(miRTarget, ExpData = NULL, mres = NULL, minSharedmiR = 3, poscorcutoff = 0,
    num_perm = 100, padjustvaluecutoff = 0.01, padjustmethod = "BH", senscorcutoff = 0.3, scorecutoff = 0.5,
    method = c("miRHomology", "pc", "sppc", "ppc", "hermes", "muTaME", "cernia")) {

    if (method == "miRHomology") {
        ceRInt <- miRHomology(miRTarget, minSharedmiR = minSharedmiR, padjustvaluecutoff = padjustvaluecutoff,
            padjustmethod = padjustmethod)
    } else if (method == "pc") {
        ceRInt <- pc(miRTarget, ExpData, minSharedmiR = minSharedmiR, poscorcutoff = poscorcutoff,
            padjustvaluecutoff = padjustvaluecutoff, padjustmethod = padjustmethod)
    } else if (method == "sppc") {
        ceRInt <- sppc(miRTarget, ExpData, minSharedmiR = minSharedmiR, poscorcutoff = poscorcutoff,
            padjustvaluecutoff = padjustvaluecutoff, padjustmethod = padjustmethod, senscorcutoff = senscorcutoff)
    } else if (method == "ppc") {
        ceRInt <- ppc(miRTarget, ExpData, minSharedmiR = minSharedmiR, num_perm = num_perm, padjustvaluecutoff = padjustvaluecutoff,
            padjustmethod = padjustmethod)
    } else if (method == "hermes") {
        ceRInt <- hermes(miRTarget, ExpData, minSharedmiR = minSharedmiR, num_perm = num_perm, padjustvaluecutoff = padjustvaluecutoff,
            padjustmethod = padjustmethod)
    } else if (method == "muTaME") {
        ceRInt <- muTaME(miRTarget, mres, minSharedmiR = minSharedmiR, padjustvaluecutoff = padjustvaluecutoff,
            padjustmethod = padjustmethod, scorecutoff = scorecutoff)
    } else if (method == "cernia") {
        ceRInt <- cernia(miRTarget, ExpData, mres, minSharedmiR = minSharedmiR, poscorcutoff = poscorcutoff,
            padjustvaluecutoff = padjustvaluecutoff, padjustmethod = padjustmethod, scorecutoff = scorecutoff)
    }
    return(ceRInt)
}

## Validation of computationally predicted miRNA sponge interactions.
## Validation of computationally predicted miRNA sponge interactions.
spongeValidate <- function(spongenetwork, directed = FALSE, Groundtruth) {
  
  spongenetwork_graph <- graph_from_data_frame(spongenetwork, directed = directed)
  Groundtruth_graph <- graph_from_data_frame(Groundtruth, directed = directed)
  Validated_interactions <- as_data_frame(spongenetwork_graph %s% Groundtruth_graph)
  colnames(Validated_interactions) <- c("sponge_1", "sponge_2")
  
  return(Validated_interactions)
}

## netModule function for identifying miRNA sponge modules from network.  Possible methods
## include FN, MCL, LINKCOMM and MCODE.
netModule <- function(spongenetwork, method = "MCL", directed = FALSE, modulesize = 3, save = FALSE) {
  
  if (method == "FN" | method == "MCL" | method == "MCODE") {
    spongenetwork_Cluster <- cluster(graph_from_data_frame(spongenetwork, directed = directed),
                                             method = method, directed = directed)
  } else if (method == "LINKCOMM") {
    edgelist <- get.edgelist(graph_from_data_frame(spongenetwork, directed = directed))
    spongenetwork_Cluster <- getLinkCommunities(edgelist, directed = directed)$nodeclusters
  }
  
  if (method == "FN" | method == "MCL") {
    spongenetwork_Cluster_result <- lapply(seq_len(max(spongenetwork_Cluster)),
                                           function(i) rownames(as.matrix(spongenetwork_Cluster))[which(spongenetwork_Cluster == i)])
    size <- unlist(lapply(seq_len(max(spongenetwork_Cluster)), function(i) length(spongenetwork_Cluster_result[[i]])))
    spongenetwork_Cluster_result <- lapply(which(size >= modulesize), function(i) spongenetwork_Cluster_result[[i]])
  } else if (method == "LINKCOMM") {
    spongenetwork_Cluster_result <- lapply(seq_len(max(c(spongenetwork_Cluster$cluster))), 
                                           function(i) as.character(spongenetwork_Cluster$node[which(c(spongenetwork_Cluster$cluster) == i)]))
    size <- unlist(lapply(seq_len(max(c(spongenetwork_Cluster$cluster))), function(i) length
                          (spongenetwork_Cluster_result[[i]])))
    spongenetwork_Cluster_result <- lapply(which(size >= modulesize), function(i) spongenetwork_Cluster_result[[i]])
  } else if (method == "MCODE") {
    spongenetwork_Cluster <- spongenetwork_Cluster + 1
    spongenetwork_Cluster_result <- lapply(seq_len(max(spongenetwork_Cluster)),
                                           function(i) rownames(as.matrix(spongenetwork_Cluster))[which(spongenetwork_Cluster == i)])
    size <- unlist(lapply(seq_len(max(spongenetwork_Cluster)), function(i) length(spongenetwork_Cluster_result[[i]])))
    spongenetwork_Cluster_result <- lapply(which(size >= modulesize), function(i) spongenetwork_Cluster_result[[i]])
  }
  
  if (save) {
    if (method == "FN" | method == "MCL" | method == "LINKCOMM") {
      res <- spongenetwork_Cluster
    } else if (method == "MCODE") {
      res <- spongenetwork_Cluster + 1
    }
    fileName <- paste("spongenetwork_Cluster_", method, ".txt", sep = "")
    spongenetwork_Cluster_Name <- list()
    k <- 0
    
    for (i in seq_len(max(res))) {
      k <- k + 1
      spongenetwork_Cluster_Name[[k]] <- rownames(as.matrix(res))[which(res == i)]
      cat(c(k, "\t", length(which(res == i))), file = fileName, sep = "", append = TRUE)
      for (j in which(res == i)) {
        cat(c("\t", rownames(as.matrix(res))[j]), file = fileName, sep = "", append = TRUE)
      }
      if (i != max(res)) {
        cat("\n", file = fileName, sep = "", append = TRUE)
      }
    }
  }
  
  return(spongenetwork_Cluster_result)
  
}


## Disease enrichment analysis of modules
moduleDEA <- function(Modulelist, OrgDb = "org.Hs.eg.db", ont = "DO", padjustvaluecutoff = 0.05,
    padjustedmethod = "BH") {

    entrezIDs <- lapply(seq_along(Modulelist), function(i) bitr(Modulelist[[i]], fromType = "SYMBOL",
        toType = "ENTREZID", OrgDb = OrgDb)$ENTREZID)

    entrezIDs <- lapply(seq_along(Modulelist), function(i) as.character(entrezIDs[[i]]))

    enrichDOs <- lapply(seq_along(Modulelist), function(i) enrichDO(entrezIDs[[i]], ont = ont, pvalueCutoff = padjustvaluecutoff,
        pAdjustMethod = padjustedmethod))

    enrichDGNs <- lapply(seq_along(Modulelist), function(i) enrichDGN(entrezIDs[[i]], pvalueCutoff = padjustvaluecutoff,
        pAdjustMethod = padjustedmethod))

    enrichNCGs <- lapply(seq_along(Modulelist), function(i) enrichNCG(entrezIDs[[i]], pvalueCutoff = padjustvaluecutoff,
        pAdjustMethod = padjustedmethod))

    return(list(enrichDOs, enrichDGNs, enrichNCGs))
}

## Functional GO, KEGG and Reactome enrichment analysis of modules
moduleFEA <- function(Modulelist, ont = "BP", KEGGorganism = "hsa", Reactomeorganism = "human",
    OrgDb = "org.Hs.eg.db", padjustvaluecutoff = 0.05, padjustedmethod = "BH") {

    entrezIDs <- lapply(seq_along(Modulelist), function(i) bitr(Modulelist[[i]], fromType = "SYMBOL",
        toType = "ENTREZID", OrgDb = OrgDb)$ENTREZID)

    entrezIDs <- lapply(seq_along(Modulelist), function(i) as.character(entrezIDs[[i]]))

    enrichGOs <- lapply(seq_along(Modulelist), function(i) enrichGO(entrezIDs[[i]], OrgDb = OrgDb,
        ont = ont, pvalueCutoff = padjustvaluecutoff, pAdjustMethod = padjustedmethod))

    enrichKEGGs <- lapply(seq_along(Modulelist), function(i) enrichKEGG(entrezIDs[[i]], organism = KEGGorganism,
        pvalueCutoff = padjustvaluecutoff, pAdjustMethod = padjustedmethod))

    enrichReactomes <- lapply(seq_along(Modulelist), function(i) enrichPathway(entrezIDs[[i]], organism = Reactomeorganism,
        pvalueCutoff = padjustvaluecutoff, pAdjustMethod = padjustedmethod))

    return(list(enrichGOs, enrichKEGGs, enrichReactomes))

}

## Survival analysis of modules
moduleSurvival <- function(Modulelist, ExpData, SurvData, devidePercentage = 0.5, plot = FALSE) {

    ExpDataNames <- c(as.matrix(ExpData[1, ]))
    ExpData <- unfactor(ExpData[-1, ])    
    colnames(ExpData) <- ExpDataNames
    myfit <- list()
    LogRank <- list()

    for (i in seq_along(Modulelist)) {
        Interin_Data <- cbind(SurvData[, seq(2, 3)], ExpData[, which(ExpDataNames %in% Modulelist[[i]])])
        Interin_Data <- na.omit(Interin_Data)

        try_mm <- try(coxph(survival::Surv(time, status) ~ ., data = data.frame(Interin_Data)),
            silent = TRUE)
        if ("try-error" %in% class(try_mm))
            next

        mm <- coxph(survival::Surv(time, status) ~ ., data = data.frame(Interin_Data))

        Risk_score <- predict(mm, newdata = data.frame(Interin_Data), type = "risk")

        group <- rep("NA", dim(Interin_Data)[1])
        group[Risk_score > quantile(Risk_score, probs = devidePercentage)] <- "High"
        group[Risk_score <= quantile(Risk_score, probs = devidePercentage)] <- "Low"

        Data <- cbind(Interin_Data[, seq_len(2)], group)
        myfit[[i]] <- survfit(survival::Surv(time, status) ~ group, data = Data)

        sdf <- survdiff(survival::Surv(time, status) ~ group, data = Data)
        sdf.p.val <- 1 - pchisq(sdf$chisq, length(sdf$n) - 1)
        HR <- (sdf$obs[1]/sdf$exp[1])/(sdf$obs[2]/sdf$exp[2])
        HRlow95 <- exp(log(HR) - qnorm(0.975) * sqrt(1/sdf$exp[1] + 1/sdf$exp[2]))
        HRup95 <- exp(log(HR) + qnorm(0.975) * sqrt(1/sdf$exp[1] + 1/sdf$exp[2]))

        LogRank[[i]] <- c(sdf$chisq, sdf.p.val, HR, HRlow95, HRup95)
    }

    if (plot) {
        for (i in seq_along(myfit)) {
            if (!is.null(LogRank[[i]])) {
                dev.new()
                plot(myfit[[i]], lty = 1, col = c("red", "green"), main = paste("Module", i), xlab = "Time (Months)",
                    ylab = "Probability of survival")

                legend("topright", legend = c("High risk group", "Low risk group"), lty = seq_len(2),
                    col = c("red", "green"))
            }
        }
    }

    LogRank_res <- do.call(rbind, LogRank)

    if (length(myfit) >= 1) {
        colnames(LogRank_res) <- c("Chi-square", "p-value", "HR", "HRlow95", "HRup95")
        names(LogRank) <- seq_along(Modulelist)
        LogRank[sapply(LogRank, is.null)] <- NULL
        rownames(LogRank_res) <- paste("Module", names(LogRank))
    }

    return(LogRank_res)
}
zhangjunpeng411/miRsponge documentation built on Aug. 13, 2018, 4:46 a.m.