# R/network.diffusion.R In SIMLR: Single-cell Interpretation via Multi-kernel LeaRning (SIMLR)

# perform network diffusion of K steps over the network A
"network.diffusion" <- function( A, K ) {

# set the values of the diagonal of A to 0
diag(A) = 0

# compute the sign matrix of A
sign_A = A
sign_A[which(A>0,arr.ind=TRUE)] = 1
sign_A[which(A<0,arr.ind=TRUE)] = -1

# compute the dominate set for A and K
P = dominate.set(abs(A),min(K,nrow(A)-1)) * sign_A

# sum the absolute value of each row of P
DD = apply(abs(P),MARGIN=1,FUN=sum)

# set DD+1 to the diagonal of P
diag(P) = DD + 1

# compute the transition field of P
P = transition.fields(P)

# compute the eigenvalues and eigenvectors of P
eigen_P = eigen(P)
U = eigen_P\$vectors
D = eigen_P\$values

# set to d the real part of the diagonal of D
d = Re(D + .Machine\$double.eps)

# perform the diffusion
alpha = 0.8
beta = 2
d = ((1-alpha)*d)/(1-alpha*d^beta)

# set to D the real part of the diagonal of d
D = array(0,c(length(Re(d)),length(Re(d))))
diag(D) = Re(d)

# finally compute W
W = U %*% D %*% t(U)
diagonal_matrix = array(0,c(nrow(W),ncol(W)))
diag(diagonal_matrix) = 1
W = (W * (1-diagonal_matrix)) / apply(array(0,c(nrow(W),ncol(W))),MARGIN=2,FUN=function(x) {x=(1-diag(W))})
diag(D) = diag(D)[length(diag(D)):1]
W = diag(DD) %*% W
W = (W + t(W)) / 2

W[which(W<0,arr.ind=TRUE)] = 0

return(W)

}

# compute the dominate set for the matrix aff.matrix and NR.OF.KNN
"dominate.set" <- function( aff.matrix, NR.OF.KNN ) {

# create the structure to save the results
PNN.matrix = array(0,c(nrow(aff.matrix),ncol(aff.matrix)))

# sort each row of aff.matrix in descending order and saves the sorted
# array and a collection of vectors with the original indices
res.sort = apply(t(aff.matrix),MARGIN=2,FUN=function(x) {return(sort(x, decreasing = TRUE, index.return = TRUE))})
sorted.aff.matrix = t(apply(as.matrix(1:length(res.sort)),MARGIN=1,function(x) { return(res.sort[[x]]\$x) }))
sorted.indices = t(apply(as.matrix(1:length(res.sort)),MARGIN=1,function(x) { return(res.sort[[x]]\$ix) }))

# get the first NR.OF.KNN columns of the sorted array
res = sorted.aff.matrix[,1:NR.OF.KNN]

# create a matrix of NR.OF.KNN columns by binding vectors of
# integers from 1 to the number of rows/columns of aff.matrix
inds = array(0,c(nrow(aff.matrix),NR.OF.KNN))
inds = apply(inds,MARGIN=2,FUN=function(x) {x=1:nrow(aff.matrix)})

# get the first NR.OF.KNN columns of the indices of aff.matrix
loc = sorted.indices[,1:NR.OF.KNN]

# assign to PNN.matrix the sorted indices
PNN.matrix[(as.vector(loc)-1)*nrow(aff.matrix)+as.vector(inds)] = as.vector(res)

# compute the final results and return them
PNN.matrix = (PNN.matrix + t(PNN.matrix))/2

return(PNN.matrix)

}

# compute the transition field of the given matrix
"transition.fields" <- function( W ) {

# get any index of columns with all 0s
zero.index = which(apply(W,MARGIN=1,FUN=sum)==0)

# compute the transition fields
W = dn(W,'ave')

w = sqrt(apply(abs(W),MARGIN=2,FUN=sum)+.Machine\$double.eps)
W = W / t(apply(array(0,c(nrow(W),ncol(W))),MARGIN=2,FUN=function(x) {x=w}))
W = W %*% t(W)

# set to 0 the elements of zero.index
W[zero.index,] = 0
W[,zero.index] = 0

return(W)

}

# normalizes a symmetric kernel
"dn" = function( w, type ) {

# compute the sum of any column
D = apply(w,MARGIN=2,FUN=sum)

# type "ave" returns D^-1*W
if(type=="ave") {
D = 1 / D
D_temp = matrix(0,nrow=length(D),ncol=length(D))
D_temp[cbind(1:length(D),1:length(D))] = D
D = D_temp
wn = D %*% w
}
# type "gph" returns D^-1/2*W*D^-1/2
else if(type=="gph") {
D = 1 / sqrt(D)
D_temp = matrix(0,nrow=length(D),ncol=length(D))
D_temp[cbind(1:length(D),1:length(D))] = D
D = D_temp
wn = D %*% (w %*% D)
}
else {
stop("Invalid type!")
}

return(wn)

}

## Try the SIMLR package in your browser

Any scripts or data that you put into this service are public.

SIMLR documentation built on Nov. 8, 2020, 5:40 p.m.