add procedure to calculate the average rank of genes over different threshold

5-create-community-network.R

@@ -22,10 +22,10 @@
   library(Hmisc)
   
 # determine the stage and transband ####
-  dev.stage <- 'im'
+  dev.stage <- 'rp'
   chr <- 5
-  start <- 6000000
-  end <- 8000000
+  start <- 24000000
+  end <- 26000000
   n.perm <- 1000
   
 # load required data ####
@@ -291,7 +291,30 @@
   plot(threshold.m$threshold, 
        threshold.m$node, 
        xlab = "rank threshold", 
-       ylab = 'number of node') # 1552 for RP5 and 1248 for PD3
+       ylab = 'number of node') # 1552 for RP4 and 1248 for PD2
+  
+  # determine the stability of rank (RP4)
+  
+  degree.table <- data.frame(candidates = candidates)
+  
+  for(i in max(round(edges$interaction, 0)):1) {
+    print(i)
+    edge.thresholded <- filter(edges, interaction <= i)
+    degree.table.tmp <- as.data.frame(table(edge.thresholded$source))
+    colnames(degree.table.tmp) <- c("candidates", i)
+    rank.degree.table <- merge(data.frame(candidates = candidates),
+                               degree.table.tmp, all.x = T)
+    rank.degree.table$rank <- rank(-rank.degree.table[, 2], 
+                                   na.last = T, ties.method = 'average')
+    degree.table[, paste0(i)] <- rank(-rank.degree.table[, 2], 
+                                      na.last = T, ties.method = 'average')
+  }
+  
+  all.rank <- data.frame(candidates,
+                         mean = round(apply(degree.table[2:ncol(degree.table)], 1, mean), 2),
+                         sd = round(apply(degree.table[2:ncol(degree.table)], 1, sd), 2))
+  all.rank <- all.rank[order(all.rank$mean), ]
+  write.csv(all.rank, paste0('networks/', dev.stage, chr, '-ranks.csv'))
   
   # visualize the network
   edges.th <- filter(edges, interaction <= 288)