deeenes

some little help for you, lovely bioinformatics students :)

Jan 2nd, 2013
110
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # loading the package
  2. require(igraph)
  3.  
  4. # some little modification on the graph merging functions
  5. # can be found here: http://pastebin.com/hiu9dbEH
  6. source("merge_functions.r")
  7.  
  8. # reading input files
  9. notch <- read.graph("notch1.xml",format="graphml")
  10. smad <- read.graph("smad2.xml",format="graphml")
  11. axin <- read.graph("axin.xml",format="graphml")
  12.  
  13. # those graph's vertices do have `name' attribute:
  14. # V(axin)$name
  15. # that's required for graph.intersect.by.name()
  16. # and graph.union.by.name() functions to work
  17.  
  18. # merge just two:
  19. ns <- graph.union.by.name(notch,smad)
  20. # and then the union of two with the third:
  21. nsa <- graph.union.by.name(axin,ns)
  22.  
  23. # let's make membership matrix
  24. # iterating on the common graph's vertices:
  25. i <- 1
  26. n <- a <- s <- m <- q <- NULL
  27. for(v in V(nsa)$name){
  28.     ifelse(length(which(V(axin)$name==v))==0,aa <- FALSE,aa <- TRUE)
  29.     ifelse(length(which(V(smad)$name==v))==0,ss <- FALSE,ss <- TRUE)   
  30.     ifelse(length(which(V(notch)$name==v))==0,nn <- FALSE,nn <- TRUE)
  31.     a <- c(a,aa)
  32.     s <- c(s,ss)
  33.     n <- c(n,nn)
  34. }
  35.  
  36. m <- cbind(a,s,n)
  37.  
  38. # this is a little trick,
  39. # we can do it only because we have not more than three networks
  40. # this is how it works:
  41. # if the node is member of the axin network, we add 1 to it's value
  42. # if it is in the smad2 network, then we add 2
  43. # if in case of notch1 network membership, then the increment is 4
  44. # so all the 7 possibilities gonna be represented by an individual number
  45. # e.g. 7 means the node is member in all the three networks
  46. # or 6 means it's member of smad2 and notch1, but missing from the axin network
  47. # yeah, this is the same how unix file privilege notation works, btw :)
  48. for(v in V(nsa)$name){
  49.     ifelse(length(which(V(axin)$name==v))==0,aa <- 0,aa <- 1)
  50.     ifelse(length(which(V(smad)$name==v))==0,ss <- 0,ss <- 2)  
  51.     ifelse(length(which(V(notch)$name==v))==0,nn <- 0,nn <- 4)
  52.     q <- c(q,sum(aa,ss,nn))
  53. }
  54.  
  55. colnames(m) <- c("axin1","smad2","notch1")
  56.  
  57. # as palettes are vectors with integer indexes,
  58. # this integer vector `q' gives an easy way to  
  59. # select different colors for the 7 groups:
  60. rcol <- rainbow(7)
  61. V(nsa)$color <-rcol[q]
  62. V(nsa)$label <- V(nsa)$name
  63. nsa$layout <- layout.kamada.kawai
  64. E(nsa)$curved <- TRUE
  65. cairo_pdf(filename="merged.pdf",
  66.     onefile=TRUE,width=12,height=12,pointsize=12)
  67.     plot(nsa)
  68. dev.off()
  69.  
  70. # you should add some legend with the legend() function
  71. # of course before devoff()
  72. # please read the manual: ?legend or help(legend)
  73. # and don't copy, please understand well every bit of R functions u use!
  74. # denes »» [email protected] »»
Advertisement
Add Comment
Please, Sign In to add comment