2

I am inexperienced with iteration in R and am hoping to speed up a process as I am implementing some analysis in a website. I realise that this may be a Q for SO but I couldn't get an answer from them...

I found a very useful tutorial that allows me to iterate through a matrix, pick out the criteria data above a certain threshold (>0.01) and then permeate three vectors: "source, target and corr" with these values to eventually make a nicely organised data frame.

adjacency<-adjacency(dat)
m<-adjacency
source=c()
target=c()
corr<-c()
g1<-rownames(adjacency)[1:dim(m)[1]]
g2<-g1
for(gene in g1){
  for(gen in g2){
    if(m[gene,gen]>0.001){
      source<-c(source,gene)
      target<-c(target,gen)
      corr<-c(corr,m[gene,gen])
    }
  }
}
network<-data.frame(source,target,corr)

However for big data frames this code is quite slow and I switched to vectorised format. And yet when I try it in vectorised format its output is different...

network2 = which(m > 0.001, arr.ind = TRUE)
source = network2[,1]
target  = network2[,2]
corr = m[m > 0.001]
network2<-cbind(network2, corr)
colnames(network2)<-c("source", "target", "corr")

My problem is that when I run these two pieces of code on the same vector

>network
   source target        corr
1      n1     n1 1.000000000
2      n1     n3 0.001466542
3      n1     n7 0.001164927
4      n1     n9 0.004027818
5      n2     n2 1.000000000
6      n2     n9 0.002580978

The above is the expected output from the for loop code:

However for big data frames this code is quite slow and I switched to vectorised format. And yet when I try it in vectorised format its output is different...

network2 = which(m > 0.001)
source = network2[,1]
target  = network2[,2]
corr = m[m > 0.001]
network2<-cbind(network2, corr)
colnames(network2)<-c("source", "target", "corr")

The problem is that I get the following error:

matching_indices = which(m > 0.01)
source = matching_indices[,1]

Error in matching_indices[, 1] : incorrect number of dimensions

This is because I have not used the argument 'are.ind=TRUE' in the 'match' command. But specifying this gives undesired (as it takes only the self-interacting nodes into account) output... so I'm wondering what I should do about this?

The undesired output is as follows (notice that 'corr' is unanimously 1, only taking self replicating nodes into account):

> network2
    source target corr
n1       1      1    1
n2       2      2    1
n3       3      3    1
n4       4      4    1
n5       5      5    1
n6       6      6    1

How do I get the output of the vectorised format code to be the same as the the for loop kind?

example data can be stored in the variable 'dat' with the following code:

library('dplyr')
library('igraph')
library('RColorBrewer')

set.seed(1)

# generate a couple clusters
nodes_per_cluster <- 30
n <- 10

nvals <- nodes_per_cluster * n

# cluster 1 (increasing) 
cluster1 <- matrix(rep((1:n)/4, nodes_per_cluster) + 
                   rnorm(nvals, sd=1),
                   nrow=nodes_per_cluster, byrow=TRUE)

# cluster 2 (decreasing)
cluster2 <- matrix(rep((n:1)/4, nodes_per_cluster) + 
                   rnorm(nvals, sd=1),
                   nrow=nodes_per_cluster, byrow=TRUE)

# noise cluster
noise <- matrix(sample(1:2, nvals, replace=TRUE) +
                rnorm(nvals, sd=1.5),
                nrow=nodes_per_cluster, byrow=TRUE)

dat <- rbind(cluster1, cluster2, noise)
colnames(dat) <- paste0('n', 1:n)
rownames(dat) <- c(paste0('cluster1_', 1:nodes_per_cluster), 
                   paste0('cluster2_', 1:nodes_per_cluster),
                   paste0('noise_',    1:nodes_per_cluster))
0

I couldn't find adjacency, so I'll use a sample matrix:

set.seed(1)
m <- matrix(runif(100), nr=10)
network2 <- which(m < 0.1, arr.ind = TRUE)

I think this may be the crux of your difference:

corr <- m[network2]
cbind(network2, corr)
#      row col       corr
# [1,]  10   1 0.06178627
# [2,]   7   3 0.01339033
# [3,]   7   5 0.02333120
# [4,]   5   6 0.07067905
# [5,]   6   6 0.09946616
# [6,]   9   7 0.08424691
# [7,]   2  10 0.05893438

You can rename and restructure this however you need.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.