Entrez Gene IDからPubMedを検索して共起関係をGraphで表示する

id:dakiyamaさんの日記から引用した関数。既にはてなにいないようでリンク切れしていた。

# Entrez Gene IDからPubMedを検索して共起関係をGraphで表示する
gene2pubnet <- function(ID=c(3651,389692),fname="tmp",pdf=TRUE,out=TRUE,width=10,height=10){
 library(RCurl)
 library(XML)
 library(Rgraphviz)
 library(biocGraph)
 library(org.Hs.eg.db)
 ID2symbol <- as.list(org.Hs.egSYMBOL)
 # EntrezIDとSymbolの対応vector
 gene_name <- sort(unlist(ID2symbol[names(ID2symbol)%in%ID]))
 url  <- paste("http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=gene&mode=xml&id",ID,sep="=")
 # EntrezIDを含むPubMedの検索式
 gene <- list()
 for(i in 1:length(ID)){
  x <- getURL(url[i])
  x <- xmlRoot(xmlTreeParse(x))
  x <- xmlElementsByTagName(x,"PubMedId",recursive=TRUE)
  x <- unique(sapply(x,xmlValue))
  # PubMedIDの抜きだし
  y <- paste("http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&term=",x,"[PMID]+AND+human[ORGN]",sep="")
  # EntrezIDを含むPubMedに記載される全Symbolの検索式
  y <- getURL(y)
  y <- lapply(y,function(x) xmlRoot(xmlTreeParse(x)))
  y <- lapply(y,function(x) xmlSApply(x[["IdList"]], xmlValue))
  y <- lapply(y,function(x) unlist(ID2symbol[x]))
  # Symbolの抜きだし
  names(y) <- x
  # Symbolの名前にPubMedIDをつける -> 文献に含まれるGeneを明示
  a <- sapply(y,length)
  ok <- a > 1 & a < 5
  # 共起数が2以上を抽出
  y <- y[ok]
  gene[[i]] <- y
 }
 x <- sapply(unique(unlist(lapply(gene,names))),list)
 y <- lapply(x,function(x,y) unique(y[sub("\\..*","",names(y))%in%x]),y=unlist(gene))
 node <- sort(c(unique(unlist(y)),names(y)))
 # Symbolと文献の一覧 -> nodeになる
 node <- c(node[!node%in%gene_name],node[node%in%gene_name])
 edge <- lapply(node,function(z)list(edges=NULL))
 # Symbolと文献の数分のedgeを準備
 names(edge) <- node
 for(j in names(y)){
  edge[[j]] <- list(edges=match(y[[j]],node))
 }
 graph <- new("graphNEL",nodes=node,edgeL=edge,edgemode="directed")
  
 colx <- c(rep(2,length(y)),rep(1,length(unique(unlist(y)))))
 # Symbolを1、文献を2
 colx <- c("pink","skyblue")[colx]
 # Symbolをpink、文献をskyblue
 names(colx) <- node
 colx[names(colx)%in%gene_name] <- "red"
 # 検索したEntrezIDをred
 nodeRenderInfo(graph) <- list(fill=colx)
 if(out){
  out <- edgeNames(graph)
  out <- matrix(unlist(strsplit(out,split="~")),ncol=2,byrow=TRUE)
  colnames(out) <- c("node","target")
  write.table(out, paste(fname,"out",sep="_"),sep="\t",row.names=FALSE)
 }
 if(pdf){
  pdf(paste(fname,".pdf",sep=""),width=width,height=height)
  graph.par(list(nodes=list(fontsize=14)))
  renderGraph(layoutGraph(graph,layoutType="neato"))
  dev.off()
 }
 png(paste(fname,".png",sep=""),width=width*480/7,height=height*480/7)
 graph.par(list(nodes=list(fontsize=14)))
 graph <- layoutGraph(graph,layoutType="neato")
 graph <- renderGraph(graph)
 dev.off()
 con <- openHtmlPage(paste(fname,".html",sep=""),title=fname)
 gnodes <- nodes(graph)[colx!="skyblue"]
 gnodes_id <- unlist(as.list(org.Hs.egSYMBOL2EG[gnodes]))
 pnodes <- nodes(graph)[colx=="skyblue"]
 gpnodes <- c(gnodes,pnodes)
 links <- character(length(gpnodes))
 tooltips <- gpnodes
 gene_links <- paste("http://www.ncbi.nlm.nih.gov/sites/entrez?db=gene&cmd=Retrieve&dopt=full_report&list_uids=",gnodes_id,sep="")
 pubmed_links <- paste("http://www.ncbi.nlm.nih.gov/pubmed/",pnodes,sep="")
 links <- c(gene_links,pubmed_links)
 names(links) <- gpnodes
 names(tooltips) <- gpnodes
 tags <- list(HREF=links,TITLE=tooltips)
 imageMap(graph,con=con,tags=tags,imgname=paste(fname,".png",sep=""))
 closeHtmlPage(con)
 return(y)
}

どこにいますか?お礼が言いたい。活用してますよ。