########################################################################
# MedlineR: an open source library in R 
#   for Medline literature data mining
#
########################################################################

## History of this library 
##  XharvestOAI # 7-11-03
##  R- caBIG    # 8-7-03    12-8-03
##  MedlineR  
##   written  by Simon Lin  
##                 # v00, v10     alpha release      # 1-15-04   2-21-04
##                 # v20          public release     # 3-25-04
##   reviewed by Patrick McConnell, Kim Johnson, and Jennifer Shoemaker 
##                 # 3-18-04
##   tested by   Simon Lin
##                 # 3-18-04  using R1.6.2

## Table of Contents
##   Initialize the package
##   Count the number of abstracts with both term1 and term2
##   pause between queries, according to the NCBI rule
##   return a co-occurrence matrix for a list of terms
##   output the co-occurrence matrix in Pajek format for visualization
##   Retrieve a MEDLINE abstract
##   query NCBI Pubmed
##   encode URI
##   An example Application

########################################################################
# Initialize the package
#
########################################################################

# DEBUG # setwd ('C:\\data\\Text-DM\\MedlineR\\R')

# Check if the required XML package is installed
  require (XML) || stop ("Need the XML package from CRAN!")

# Set up the base URL for NCBI Pubmed http service
# This service is provided by Entrez.
#
# The documentation of Entrez is at
# http://eutils.ncbi.nlm.nih.gov/entrez/query/static/eutils_help.html
#

# Set the URL for the NCBI Pubmed http service
#  This URL is stored as a
#  global variable in R called "serviceUrl.entrez".
#
#  Note: the global variable in R is 
#  called options. See help (options) for
#  details
#
#  Example:
# 
#    options   ("serviceUrl.entrez" = 
#               "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/")
#    getOption ("serviceUrl.entrez")
#
#  Note: The syntax of these operations is similar to setwd(), getwd()
  options("serviceUrl.entrez" = 
          "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/")

########################################################################
# Count the number of abstracts with both term1 and term2
#
########################################################################

#
# Examples:
#   countApair(term1="STE18", term2="STE12")
#   countApair(term1="SST2", term2="STE12",
#              termAdditional= "+AND+(Saccharomyces+OR+yeast")
#
# How does it work?
#  The actural http request is translated into
#  http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=pubmed&
#   term=STE12+AND+SST2+AND+Saccharomyces
#

countApair<- function (
    term1, term2,                           # a pair of terms
    termAdditional="",                      # additional modifying terms
    baseUrl=getOption("serviceUrl.entrez")  # URL of the Pubmed service
    ) {
    
    # QC: make sure the baseUrl is all right.
    if (is.null (baseUrl)) {
        stop ("Need to define the URL of the Pubmed service!")
    } 
     
  # Get the query string ready. This string should be in the 
  # Pubmed syntax. The Pubmed syntax is documented at
  # http://eutils.ncbi.nlm.nih.gov/entrez/query/static/esearch_help.html
    query<- paste (baseUrl, 
                   "esearch.fcgi?",
                   "db=pubmed&",             # database is Pubmed
                   "rettype=count&",         # simply return the counts
                   "term=",                  #  using these terms
                     term1, "+AND+", term2, termAdditional,
                   sep="")  
    result.xml<- try (xmlTreeParse(file=query, isURL=T))   # in XML
    # count is parsed from XML
    count<- as.numeric(xmlValue (xmlRoot (result.xml) [["Count"]]))
    return (count)
}

########################################################################
# pause between queries, according to the NCBI rule
#  
########################################################################

# Note:
# According to the NCBI rule, query larger than 100 requests
#  should be 'nicely' run off hours and weekends
#
# The NCBI rule is documented at
#  http://eutils.ncbi.nlm.nih.gov/entrez/query/static/eutils_help.html
# It prohibits WebBot abuse of the http service
#
# Note: How to sleep for 3 seconds in R?
#  Sys.sleep (time=3)
#

pauseBetweenQueries<- function (
   sleep.peak=15,                # pause (in seconds) during peak hours
   sleep.offpeak=3               # pause (in seconds) during off-peak 
  ) {

 # sleep.peak<-15; sleep.offpeak<-3

 # Date example:
 # "Thu"      "Jan"      "15"       "16:46:11" "2004" 
 result.date<- unlist (strsplit(
   date(), split=" "))
 hour<- as.numeric(unlist (strsplit (result.date[4], split=':'))[1])
 
 # off peak hours are Sat, Sun or anytime between 9 pm and 5 am
 if (
   (result.date[1]=="Sat") | (result.date[1]=="Sun") |
   (hour > 21) | (hour<5)
  ) {off.peak<-T} else {off.peak<-F}

 # perform the sleep
 if (off.peak) {
  Sys.sleep (sleep.offpeak)
 } else {
  Sys.sleep (sleep.peak)
 }
}

########################################################################
# return a co-occurrence matrix for a list of terms
#  
########################################################################

# Example: 
#  find the co-occurance matrix of these terms in yeast
#
#  termList<- c("STE18", "DIG1", "STE12", "SST2", "HOG1", 
#               "pheromone", "cell+cycle")
#  result.matrix<- getAmatrix(termList,
#    termAdditional= "+AND+(Saccharomyces+OR+yeast")
#  # after building the matrix you can get an image for visualization
#  image (result.matrix)
#  image (result.matrix^0.2, col=terrain.colors(16)) 
#
# Note:
#  must transform the raw data matrix for display, 
#  otherwise, the dynamic range of the color is not 
#  suitable for human recognition!

getAmatrix<- function (
  termList,                                # a list of terms
  termAdditional=""                        # additional modifier
 ) {
 # termList<- c("STE18", "DIG1", "STE12");
 # termAdditional<- "+AND+(Saccharomyces+OR+yeast"

 # initialize the co-occurance matrix
 n.terms<- length (termList)
 result.matrix<-                   
  matrix (0, ncol=n.terms, nrow=n.terms)

 # fill in the diagonal first for an optimization   
 # if a term is not found anywhere in medline
 # it certainly won't match any other terms
 for (i in 1:n.terms){               # i.e., the size of the nodes
  result.matrix[i,i]<- countApair (  #   in the network 
    term1=termList[i], 
    term2=termList[i],  
    termAdditional=termAdditional)
  pauseBetweenQueries()
 }

 # query the pairs                  # i.e., the connection of the nodes
 for (i in 1:(n.terms-1)){          #   in the network
  # check to see if the diagonal is 0
  # if so, don't process it
  if (result.matrix [i,i]==0) {next}

  for (j in (i+1):n.terms) {
     n.counts <-countApair (
         term1=termList[i], 
         term2=termList[j],  
         termAdditional=termAdditional)
     pauseBetweenQueries()
     result.matrix[i,j]<- n.counts  # the matrix is symetric
     result.matrix[j,i]<- n.counts  # i.e., uppper right v.s. lower left
                                    # copy for visualization purposes
   }
 }
 return (result.matrix)
}

########################################################################
# output the co-occurrence matrix in Pajek format for visualization
#  
########################################################################

# The output is a graph, where each node is a term connected to other
# terms by the co-occurance matrix.  The Pajek software should be used
# to visualize the output.
#
# documentation for the Pajek file format can be found at
# http://vlado.fmf.uni-lj.si/pub/networks/pajek/
#
# Note:
#   Pajek file name must end with .net
#   It is a requirement for Pajeck to recognize the file type.
#
# Example:
#  termList<- c("STE18", "DIG1", "STE12", "SST2", "HOG1", 
#               "pheromone", "cell+cycle")
#  result.matrix<- getAmatrix(termList,
#    termAdditional= "+AND+(Saccharomyces+OR+yeast")
#  write.pajek (result.matrix, termList, fileName="test.net")
#

write.pajek<- function (result.matrix,       # co-occurance matrix
                        termList=NULL,       # node name
                        nodeColorList=NULL,  # node color, default
                                             #   is yellow
                        fileName){           # must end with .net

  # Pajek files must end with .net
  # fileName<- "genes.net"

  # setup color list
  n.vertices<- dim (result.matrix)[1]
  if (is.null(nodeColorList)){               # default: yellow
    nodeColorList<- I(rep ("Yellow", n.vertices)) 
  }                                          # see help (I)

  # setup term list, default is numerical value                                     
  if (is.null(termList)){
    termList<- 1:n.vertices
  }

  # builds data frame
  df.vertices<- data.frame (
    id= 1:n.vertices,
    name= termList,
    color= I(nodeColorList)
  )

  zz <- file(fileName, "w")             # open an output file connection

  # write vertices to Pajek-formatted file
  cat ("*Vertices ", n.vertices,        # nodes ("Vertices")
       "\n", sep="", file=zz)
  for (i in 1: n.vertices) {
   cat (df.vertices$id[i], " \"", 
        as.character(df.vertices$name[i]), "\"",         # name
        " x_fact ", round (result.matrix[i,i]^0.2)+1,    # size
        " y_fact ", round (result.matrix[i,i]^0.2)+1,    #
      # " s_size ", round (result.matrix[i,i]^0.2)+1,    #
        " ic ", as.character(df.vertices$color[i]),      # color
        "\n",
        sep="", file=zz)
  }

  # write edges to Pajek-formatted file
  cat ("*Edges", "\n", sep="", file=zz) # conncetions ("Edges")

  for (i in 1:(n.vertices-1)){
   if (result.matrix [i,i]==0) {next}
   for (j in (i+1):n.vertices) {
       n.counts <-result.matrix[i,j]
       if (n.counts==0) {next}
       thickness.line<- n.counts^0.5   # transformation for visual effect
        cat (
         i, " ",                       # from 
         j, " ",                       # to
         round (thickness.line)+1,     # width
         "\n", sep="", file=zz)
    }
  }
  close(zz)
}

########################################################################
# Retrieve a MEDLINE abstract
#  
########################################################################

#
# Example:
#  fetchAnAbstract (pmID=134567)
# 

fetchAnAbstract<- function (
    pmID,                                    # PubMed ID
    baseUrl=getOption("serviceUrl.entrez")   # URL of the Pubmed service
   ){

    # pmID<- 12405121

    # QC: make sure the baseUrl is all right.
    if (is.null (baseUrl)) {
        stop ("Need to define the URL of the Pubmed service!")
    } 

  # Get the query string ready. This string should be in the Pubmed 
  # syntax. The Pubmed syntax is documented at
  # http://eutils.ncbi.nlm.nih.gov/entrez/query/static/esearch_help.html
    query<- paste (baseUrl, 
                   "efetch.fcgi?",
                   "db=pubmed&id=",
                   pmID,
                   "&retmode=xml", 
                   sep="")  
    # parse resulting XML into a tree to be returned to user
    result.xml<- try (xmlTreeParse(file=query, isURL=T))
    return (result.xml)
}

########################################################################
# query NCBI Pubmed to get the pmid (PubMed identifier)
#  
########################################################################

#
# Example:
#  queryPubmed (term="cancer+OR+diabetes")
# 

queryPubmed <- function (
    term,                                    # query term
    baseUrl=getOption("serviceUrl.entrez")   # URL of the Pubmed service
   ){
    
    # QC: make sure the baseUrl is all right.
    if (is.null (baseUrl)) {
        stop ("Need to define the URL of the Pubmed service!")
    } 
     
  # Get the query string ready. This string should be in the Pubmed 
  # syntax. The Pubmed syntax is documented at
  # http://eutils.ncbi.nlm.nih.gov/entrez/query/static/esearch_help.html
    query<- paste (baseUrl, 
                   "esearch.fcgi?",
                   "db=pubmed&term=",
                   term,
                   sep="")  
    
    # parse resulting XML into a tree to be returned to user
    result.xml<- try (xmlTreeParse(file=query, isURL=T))
    return (result.xml)
}

########################################################################
# encode URI
#  
########################################################################

# Example
#   encodeUri("sf&df:lk;lo #op")
# 
# encode a phase with strange characters in it
#
encodeUri<- function (x) {

    # x<- ""
    # input x is a PHRASE string
    # WARNING: x should not be the whole http://someCGI header!
    x<- as.character(x)

    # replcae with gsub (derived from grep!)
    x<- gsub (pattern="%",   replacement="%25", x) # must be done first!
    x<- gsub (pattern="/",   replacement="%2F", x)
    x<- gsub (pattern="\\?", replacement="%3F", x)
    x<- gsub (pattern="#",   replacement="%23", x)
    x<- gsub (pattern="=",   replacement="%3D", x)
    x<- gsub (pattern="&",   replacement="%26", x)
    x<- gsub (pattern=":",   replacement="%3A", x)
    x<- gsub (pattern=";",   replacement="%3B", x)
    x<- gsub (pattern=" ",   replacement="%20", x)
    x<- gsub (pattern="\\+", replacement="%2B", x)

    return (x)
}

########################################################################
# An Example Application:
#  (1) mines the co-occurrence association among a list of terms
#  (2) visualizes as an image
#  (3) writes a Pajek-formatted graph for visualization
# Note: Most the waiting time is spent on the mandatory waiting 
#       interval to query the Medline!
########################################################################

# Usuage:
#   findAssociation (termList=c("STE18", "DIG1", "STE12"), 
#                    fileName="test.net")
#   might take 5 minutes due to the mandatory waiting!
#

findAssociation<- function (
  termList,                         # a list of terms
  termAdditional="",                # additional modifier
  nodeColorList=NULL,                  # optional argument, node color
  fileName                          # output file name (in Pajek format)
 ) {

  # termList<- c("STE18", "DIG1", "STE12", "SST2", "HOG1")
  # termAdditional<- "+AND+(Saccharomyces+OR+yeast)"
  # fileName<- "test.net"

    association.matrix<- getAmatrix(termList)
    image (-log(association.matrix+1), col=grey(0:16/16))   # Note: log transf for visual
    write.pajek (association.matrix,           # co-occurance matrix
                 termList=termList,       # node name
                 nodeColorList=nodeColorList,   # node color
                 fileName)                # file name
}

#
# Figure 1 in the paper
#   findAssociation (termList=c("STE18", "DIG1", "HOG1", "cell+cycle", 
#                    "pheromone"), fileName="test1.net")