inst/scripts/import/uniprobe/import.R

In MotifDb: An Annotated Collection of Protein-DNA Binding Sequence Motifs

# uniprobe/import.R
#------------------------------------------------------------------------------------------------------------------------
#library (RMySQL)
library (RSQLite)
library (org.Hs.eg.db)
library (org.Mm.eg.db)
library (org.Sc.sgd.db)
library (org.Ce.eg.db)
library(tools)   # for md5sum
#------------------------------------------------------------------------------------------------------------------------
printf <- function(...) print(noquote(sprintf(...)))
#------------------------------------------------------------------------------------------------------------------------
run = function (dataDir)
{
  dataDir <- file.path(dataDir, "uniprobe")
  stopifnot(file.exists(dataDir))

  all.files = identifyFiles (file.path(dataDir, 'All_PWMs'))
  matrices = readAndParse (all.files)
  tbl.pubRef = createPublicationRefTable ()
  tbl.geneRef = createGeneRefTable (dataDir)
  tbl.md = createMetadata (matrices, tbl.pubRef, tbl.geneRef)
  stopifnot (length (matrices) == nrow (tbl.md))
  matrices = renameMatrices (matrices, tbl.md)

  serializedFile <- "uniprobe.RData"
  save (matrices, tbl.md, file=serializedFile)
  printf("saved %d matrices to %s", length(matrices), serializedFile)
  printf("next step: copy %s to <packageRoot>/MotifDb/inst/extdata, rebuild package", serializedFile)

} # run
#------------------------------------------------------------------------------------------------------------------------
createMatrixNameUniqifier = function (matrix)
{
  temporary.file <- tempfile ()
  write (as.character (matrix), file=temporary.file)
  md5sum.string <- as.character (md5sum (temporary.file))
  stopifnot (nchar (md5sum.string) == 32)
  md5.6chars = substr (md5sum.string, 29, 32)
  #unlink (temporary.file)

} # createMatrixNameUniqifier
#------------------------------------------------------------------------------------------------------------------------
parsePWMfromText = function (lines.of.text)
{
  z = lines.of.text
  stopifnot (sum (sapply (c ('A', 'C', 'T', 'G'), function (token) length (grep (token, lines.of.text)))) == 4)

    # determine the number of columns
  token.count.first.line = length (strsplit (lines.of.text [1], '\t')[[1]])
  column.count = token.count.first.line - 1  # subtract off the 'A:\t' token
  result = matrix (nrow=4, ncol=column.count, byrow=TRUE, dimnames=list (c ('A', 'C', 'G', 'T'), 1:column.count))

  for (line in lines.of.text) {
    tokens = strsplit (line, '\\s*[:\\|]') [[1]]
    nucleotide = tokens [1]
    numbers.raw = tokens [2]
    number.tokens = strsplit (numbers.raw, '\\s+', perl=T)[[1]]
    while (nchar (number.tokens [1]) == 0)
      number.tokens = number.tokens [-1]
    numbers = as.numeric (number.tokens)
    #printf ('adding %s: %s', nucleotide, list.to.string (numbers))
    result [nucleotide,] = numbers
    }

  return (result)

} # parsePWMfromText
#------------------------------------------------------------------------------------------------------------------------
extractPWMfromFile = function (filename)
{
  text = scan (filename, sep='\n', what=character (0), quiet=TRUE)
  matrix.start.lines = grep ('A\\s*[:\\|]', text)
  stopifnot (length (matrix.start.lines) == 1)
  #printf ('%50s: %s', filename, list.to.string (matrix.start.lines))
  start.line = matrix.start.lines [1]
  end.line = start.line + 3
  lines = text [start.line:end.line]
  pwm.matrix = parsePWMfromText (lines)
  return (pwm.matrix)

} # extractPWMfromFile
#------------------------------------------------------------------------------------------------------------------------
createPublicationRefTable = function ()
{
  options (stringsAsFactors=FALSE)
  tbl.ref = data.frame (folder=c('CR09', 'Cell08', 'Cell09', 'EMBO10', 'GD09', 'GR09', 'MMB08', 'NBT06', 'PNAS08', 'SCI09'))
  tbl.ref = cbind (tbl.ref, author=c('Scharer', 'Berger', 'Grove', 'Wei', 'Lesch', 'Zhu', 'Pompeani', 'Berger', 'De Silva', 'Badis'))
  tbl.ref = cbind (tbl.ref, pmid=c('19147588','18585359','19632181','20517297','19204119','19158363','18681939','16998473','18541913','19443739'))
  tbl.ref = cbind (tbl.ref, organism=c('Hsapiens', 'Mmusculus', 'Celegans', 'Mmusculus', 'Celegans', 'Scerevisiae', 'Vharveyi',
                                       'Scerevisiae;Hsapiens;Mmusculus', 'Apicomplexa', 'Mmusculus'))
  tbl.ref = cbind (tbl.ref, count=c(1, 168, 34, 25, 1, 89, 1, 5, 3, 104))
  titles = c ('Genome-wide promoter analysis of the SOX4 transcriptional network in prostate cancer cells',
              'Variation in homeodomain DNA binding revealed by high-resolution analysis of sequence preferences',
              'A Multiparameter Network Reveals Extensive Divergence between C. elegans bHLH Transcription Factors',
              'Genome-wide analysis of ETS-family DNA-binding in vitro and in vivo',
              'Transcriptional regulation and stabilization of left right neuronal identity in C. elegans',
              'High-resolution DNA-binding specificity analysis of yeast transcription factors',
              'The Vibrio harveyi master quorum-sensing regulator, LuxR...',
              'Compact, universal DNA microarrays...',
              'Specific DNA-binding by Apicomplexan AP2 transcription factors',
              'Diversity and complexity in DNA recognition by transcription factors')
  tbl.ref = cbind (tbl.ref, title=titles)

  tbl.ref

} # createPublicationRefTable
#------------------------------------------------------------------------------------------------------------------------
identifyFiles = function (filePath)
{
  stopifnot(file.exists(filePath))
  
  cmd = sprintf ('find %s -name "*pwm*"', filePath)

  files.raw = system (cmd, intern=TRUE)

      # all legit files end in ".pwm" or ".txt"
  files = c (grep (".pwm$", files.raw, value=T, ignore.case=T), 
             grep (".txt$", files.raw, value=T, ignore.case=T))

  reverseComplementMatrices = grep ('RC', files)
  if (length (reverseComplementMatrices) > 0)
    files = files [-reverseComplementMatrices]

      # don't know why these were excluded.  leave them in for now
  embo10.excluders = grep ('EMBO', files)
  if (length (embo10.excluders) > 0)
    files = files [-embo10.excluders]

  cell09.excluders = grep ('Cell09', files)   # include these once the sql tables are updated
  if (length (cell09.excluders) > 0)
    files = files [-cell09.excluders]
  secondary.excluders = grep ('secondary', files)
  if (length (secondary.excluders) > 0)
    files = files [-secondary.excluders]

  invisible (files)

} # identifyFiles
#------------------------------------------------------------------------------------------------------------------------
readAndParse = function (file.list)
{
  matrices = list ()

  for (file in file.list) {
    #printf ('read and parse %s', file)
    text = scan (file, sep='\n', what=character (0), quiet=TRUE)
    matrix.start.lines = grep ('A:', text)
    stopifnot (length (matrix.start.lines) == 1)
    start.line = matrix.start.lines [1]
    end.line = start.line + 3
    lines.of.text = text [start.line:end.line]
    pwm.matrix = parsePWMfromText (lines.of.text)
    name.tokens <- strsplit(file,"/")[[1]]
    token.count <- length(name.tokens)

    matrix.name <- paste(name.tokens[(token.count-1):token.count], collapse="/")
    matrices [[matrix.name]] = pwm.matrix
    }

  invisible (matrices)

} # readAndParse
#------------------------------------------------------------------------------------------------------------------------
# eg, ./All_PWMs/GD09/Nsy-7.pwm to simply 'Nsy-7'
#
translateFileNameToGeneName = function (long.name)
{
  if (length (grep ('/', long.name)) > 0) {
    tokens = strsplit (long.name, '/')[[1]]
    count = length (tokens)
    gene.name.raw = tokens [count]  # get the last one
    }
  else {
    gene.name.raw = long.name
    }

  gene.name = strsplit (gene.name.raw, '\\.')[[1]][1]   # remove any file suffix

  gene.name

} # test.translateFileNameToGeneName
#------------------------------------------------------------------------------------------------------------------------
# and update matrix names
createMetadata = function (matrices, tbl.pubRef, tbl.geneRef)
{
  options (stringsAsFactors=FALSE)
  dataSource = 'UniPROBE'
  trim = function (s) {sub (' +$', '', sub ('^ +', '', s))}
  tbl.md = data.frame ()
  
  removers = list (Cart1=110, Cutl1=115, Hoxa7=156, Irx3=185)
 
  for (m in 1: length (matrices)) {
    matrix.name = names (matrices) [m]
    #printf ('%d: %s', m, matrix.name)
    native.name.raw = gsub ('All_PWMs/', '', matrix.name)
    native.name = extractNativeNames (native.name.raw)

       # extractNativeNames needs to be rethought.  but for now, just hack in some fixes
    if (native.name == 'Cgd2') native.name='Cgd2_3490'   # inconsistency at uniprobe, accomodated here
    if (native.name == 'PF14') native.name='PF14_0633'
    if (native.name == 'Uncx4') native.name='Uncx4.1'

    if (!native.name %in% tbl.geneRef$name) {
      browser (text='native.name not in tbl.geneRef')
      }
    
    experiment.folder = strsplit (native.name.raw, '/')[[1]][1]
    up.id.number = subset (tbl.geneRef, name==native.name & folder_name==experiment.folder)$id
    # todo BUG here!
    full.uniprobe.id = sprintf ('UP%05d', up.id.number)
    if (length (full.uniprobe.id) != 1) {
      browser (text='full.uniprobe.id has wrong length')
      }

    geneId = subset (tbl.geneRef, name==native.name & folder_name==experiment.folder)$stdID
    if (is.na (geneId) | geneId == '')
       geneId = NA_character_

    bindingDomain = trim (subset (tbl.geneRef, name==native.name & folder_name==experiment.folder)$domain)
    organism = subset (tbl.geneRef, name==native.name & folder_name==experiment.folder)$species

      # a native name (a gene symbol) may have a dash, but for the long name, we want dashes to separate
      # the organism-dataSource-geneIdentifier matrix name.
      # so covert this here, but do not eclipse any dash-including native.names

    native.name.no.dashes = gsub ('-', '_', native.name)
    native.name.uniq = sprintf ('%s-%s-%s.%s', organism, dataSource, native.name.no.dashes, full.uniprobe.id)
    if (native.name.uniq %in% rownames (tbl.md)) {
      matrix = matrices [[m]]
      uniqifier = createMatrixNameUniqifier (matrix)
      #printf ('before, not unique: %s', native.name.uniq)
      native.name.uniq = paste (native.name.uniq, uniqifier, sep='.')
      #printf ('after, not unique: %s', native.name.uniq)
      }
      
    sequenceCount = NA_integer_
    bindingSequence = NA_character_
    tfFamily = NA_character_

    experimentType = 'protein binding microarray'
    pubmedID = subset (tbl.pubRef, folder==experiment.folder)$pmid
    #printf ('%12s: %20s', native.name, organism)

    if (is.na (geneId))
      geneIdType = NA
    else if (organism == 'Scerevisiae')
      geneIdType = 'SGD'
    else if (organism %in% c ('Mmusculus', 'Hsapiens', 'Celegans'))
      geneIdType = 'ENTREZ'
    else
      geneIdType = 'todo'

    proteinId = NA_character_
    proteinIdType = NA_character_

    proteinId.tmp = subset (tbl.geneRef, name==native.name & folder_name==experiment.folder)$uniprot
    if (!is.na (proteinId.tmp) & nchar (proteinId.tmp) > 1) {
      #printf ('getting uniprot id for %s: %s', native.name, proteinId.tmp)
      proteinId = proteinId.tmp
      proteinIdType = 'UNIPROT'
      }  # found good id in the uniprot column

    if (is.na (proteinId.tmp) | nchar (proteinId.tmp) == 0) {
      proteinId.tmp = subset (tbl.geneRef, name==native.name & folder_name==experiment.folder)$refseq_id
      if (!is.na (proteinId.tmp) & nchar (proteinId.tmp) > 1) {
        #printf ('getting refseq id for %s: %s', native.name, proteinId.tmp)
        proteinId = proteinId.tmp
        proteinIdType = 'REFSEQ'
        }
      } # need to look in the refseq column


    bindingSequence =  subset (tbl.geneRef, name==native.name & folder_name==experiment.folder)$bindingSequence
    #printf ('%12s: %40s', native.name, bindingSequence)

    new.row = list (providerName=native.name.raw,
                    providerId=full.uniprobe.id,
                    dataSource=dataSource,
                    geneSymbol=native.name,
                    geneId=geneId,
                    geneIdType=geneIdType,
                    proteinId=proteinId,
                    proteinIdType=proteinIdType,
                    organism=organism,
                    sequenceCount=NA,
                    bindingSequence=bindingSequence,
                    bindingDomain=bindingDomain,
                    tfFamily=tfFamily,
                    experimentType=experimentType,
                    pubmedID=pubmedID)
    if (native.name.uniq %in% rownames (tbl.md)) browser (text='dup row name')
    tbl.md = rbind (tbl.md, data.frame (new.row, stringsAsFactors=FALSE))
    if (length (native.name.uniq) != 1) browser (text='native.name.unique length != 1')
    rownames (tbl.md) [m] = native.name.uniq
    }

  invisible (tbl.md)

} # createMetadata
#------------------------------------------------------------------------------------------------------------------------
extractNativeNames = function (native.names.raw)
{
  name.count = length (native.names.raw)
  result = vector ('character', name.count)

  for (i in 1:name.count) {
    native.name.raw = native.names.raw [i]
    tokens = strsplit (native.name.raw, '/') [[1]]
    count = length (tokens)
    cooked.1 = native.name.raw
    if (count > 1)
      cooked.1 = tokens [length (tokens)]
    tokens = strsplit (cooked.1, '[_\\.]')[[1]]
    cooked.2 = tokens [1]
    result [i] = cooked.2
    }

  invisible (result)

} # extractNativeNames
#------------------------------------------------------------------------------------------------------------------------
# 'standard' is usually entrez gene ID.  for yeast it is orf. for fly ...
uniprotToStandardID = function (organism, uniprot.ids)
{
#  uniprot.ids = unique (uniprot.ids)

  if (!exists ('lst.yeast.uniprot')) {
    tbl.tmp = toTable (org.Sc.sgdUNIPROT)
    yeast.uniprot <<- tbl.tmp$systematic_name
    names (yeast.uniprot) <<- tbl.tmp$uniprot_id
    }

  if (!exists ('mouse.uniprot')) {
    tbl.tmp = toTable (org.Mm.egUNIPROT)
    mouse.uniprot <<- tbl.tmp$gene_id
    names (mouse.uniprot) <<- tbl.tmp$uniprot_id
    }

  if (!exists ('human.uniprot')) {
    tbl.tmp = toTable (org.Hs.egUNIPROT)
    human.uniprot <<- tbl.tmp$gene_id
    names (human.uniprot) <<- tbl.tmp$uniprot_id
    }

  if (!exists ('worm.uniprot')) {
    tbl.tmp = toTable (org.Ce.egUNIPROT)
    worm.uniprot <<- tbl.tmp$gene_id
    names (worm.uniprot) <<- tbl.tmp$uniprot_id
    }

  organism = unique (organism)
  stopifnot (length (unique (organism)) == 1)
  stopifnot (organism %in% (c ('human', 'mouse', 'yeast', 'worm')))
 
  if (organism == 'human') {
    ids = human.uniprot [uniprot.ids]
    }
  else if (organism == 'mouse') {
    ids = mouse.uniprot [uniprot.ids]
    }
  else if (organism == 'yeast') {
    ids = yeast.uniprot [uniprot.ids]
    }
  else if (organism == 'worm') {
    ids = worm.uniprot [uniprot.ids]
    }


    # embarrassing but true:  could not get this to work by creating a list directly
    # round-about solution:  make a 1-column data.frame, then construct a named list from that
  df = data.frame (uniprot=uniprot.ids, std=as.character (ids), stringsAsFactors=FALSE)
  #rownames (df) = uniprot.ids
  #result = df$stdID
  #names (result) = uniprot.ids
  return (df)

} # uniprotToStandardID
#------------------------------------------------------------------------------------------------------------------------
# see ~/v/snotes/log "* load uniprobe sql dump file (11 may 2012)" for info on the uniprobe mysql database 
# used here.  
createGeneRefTable = function (dataDir)
{
  if (!exists ('db')){
      dbFile <- file.path(dataDir, "uniprobe.sqlite")
      stopifnot(file.exists(dbFile))
      db <<- dbConnect (dbDriver("SQLite"), dbFile)
      }

  tbl.pubmed = data.frame (folder=c('CR09', 'Cell08', 'Cell09', 'EMBO10', 'GD09', 'GR09', 'MMB08', 'NBT06', 'PNAS08', 'SCI09'),
                 pmid=c('19147588', '18585359', '19632181', '20517297', '19204119', '19158363', '18681939', '16998473', '18541913', '19443739'),
                 stringsAsFactors=FALSE)
       #CR09	1	Scharer 	19147588	human	Genome-wide promoter analysis of the SOX4 transcriptional 
       #                                                        network in prostate cancer cells
       #Cell08	168	Berger  	18585359	mouse	Variation in homeodomain DNA binding revealed by high-resolution 
       #                                                        analysis of sequence preferences
       #Cell09	34	Grove   	19632181	worm	A Multiparameter Network Reveals Extensive Divergence between 
       #                                                        C. elegans bHLH Transcription Factors
       #EMBO10	25	Wei     	20517297	mouse	Genome-wide analysis of ETS-family DNA-binding in vitro and in vivo
       #GD09	1	Lesch   	19204119	worm	Transcriptional regulation and stabilization of left right neuronal 
       #                                                         identity in C. elegans
       #GR09	89	Zhu     	19158363	yeast	High-resolution DNA-binding specificity analysis of yeast transcription factors
       #MMB08	1	Pompeani	18681939	Vibrio	The Vibrio harveyi master quorum-sensing regulator, LuxR...
       #NBT06	5	Berger  	16998473	yeast;human;mouse	Compact, universal DNA microarrays...
       #PNAS08	3	De Silva	18541913	apicomplexa	Specific DNA-binding by Apicomplexan AP2 transcription factors
       #SCI09	104	Badis   	19443739	mouse	Diversity and complexity in DNA recognition by transcription factors


  tbl.gene = dbGetQuery (db, 'select * from gene_ids_public')
  colnames (tbl.gene) = c ('id', 'name', 'species', 'pub', 'type')
  tbl.genomic = dbGetQuery (db, 'select * from genomic_info') [, -c(2,3,8:11)]  # remove the longish 'description' field
  tbl.pub = dbGetQuery (db, 'select * from publication_ids') [,-3]  # remove 'full_ref' for legibility
  tbl = merge (tbl.gene, tbl.pub [, c (1,4)], by.x='pub', by.y='publication_id', all.x=TRUE)
  tbl = merge (tbl, tbl.pubmed, by.x='folder_name', by.y='folder', all.x=TRUE)
  tbl = merge (tbl, tbl.genomic, all.x=TRUE, by.x='name', by.y='gene_name')
  redundant.column = grep ('species.y', colnames (tbl))
  stopifnot (length (redundant.column) == 1)
  tbl = tbl [, -redundant.column]
  column.to.rename = grep ('species.x', colnames (tbl))
  stopifnot (length (column.to.rename) == 1)
  colnames (tbl) [column.to.rename] = 'species'
  

   # now add 'standard IDs' -- orf names for yeast, entrez geneIDs for everything else

  stdID = getAllStandardIDs (tbl)
  tbl = cbind (tbl, stdID)
  duplicates = which (duplicated (tbl [, c (1,2)]))
  if (length (duplicates) > 0)
    tbl = tbl [-duplicates,]
  
    # fix the organism (species) name, from (e.g.) "Homo sapiens" to "Hsapiens"
  tbl$species = standardizeSpeciesNames (tbl$species)
  invisible (tbl)

    # now add bindingSequence, from DBDs:  gene_name + seq, apparently the binding sequence when known, of 171 of 372 genes
  tbl.dbds = dbGetQuery (db, 'select * from DBDs')  
  dbds.list = tbl.dbds$seq
  names (dbds.list) = tbl.dbds$gene_name
  bindingSequence = rep (NA_character_, nrow (tbl))
  names (bindingSequence) = tbl$name
  shared.names = unique (intersect (tbl.dbds$gene_name, tbl$name))
  bindingSequence [shared.names] = dbds.list [shared.names]
  tbl = cbind (tbl, bindingSequence)
  invisible (tbl)
  
} # createGeneRefTable
#------------------------------------------------------------------------------------------------------------------------
# make successive species-specific calls to 'uniprotToStandardID', assembling a new column to be added to tbl.geneRef
getAllStandardIDs = function (tbl.geneRef)
{
  mouse.rows  = grep ('Mus musculus', tbl.geneRef$species)
  yeast.rows = grep ('Saccharomyces cerevisiae', tbl.geneRef$species)
  human.rows = grep ('Homo sapiens', tbl.geneRef$species)
  worm.rows = grep ('Caenorhabditis elegans', tbl.geneRef$species)

  stdID = rep (NA_character_, nrow (tbl.geneRef))

  mouse.uids = tbl.geneRef$uniprot [mouse.rows]
  yeast.uids = tbl.geneRef$uniprot [yeast.rows]
  human.uids = tbl.geneRef$uniprot [human.rows]
  worm.uids = tbl.geneRef$uniprot [worm.rows]

    # add mouse geneIDs
  tbl.mouse =  uniprotToStandardID ('mouse', mouse.uids)
  stopifnot (length (mouse.rows) == nrow (tbl.mouse))
  stdID [mouse.rows] = tbl.mouse$std
  
    # add yeast orfs
  tbl.yeast =  uniprotToStandardID ('yeast', yeast.uids)
  stopifnot (length (yeast.rows) == nrow (tbl.yeast))
  stdID [yeast.rows] = tbl.yeast$std
  
    # add human geneIDs
  tbl.human =  uniprotToStandardID ('human', human.uids)
  stopifnot (length (human.rows) == nrow (tbl.human))
  stdID [human.rows] = tbl.human$std

    # add worm geneIDs
  tbl.worm =  uniprotToStandardID ('worm', worm.uids)
  stopifnot (length (worm.rows) == nrow (tbl.worm))
  stdID [worm.rows] = tbl.worm$std

  invisible (stdID)

} # getAllStandardIDs
#------------------------------------------------------------------------------------------------------------------------
# change, e.g., "Homo sapiens" to "Hsapiens"
standardizeSpeciesNames = function (names)
{
   fix = function (name) {
     tokens = strsplit (name, ' ')[[1]]
     stopifnot (length (tokens) == 2)
     genus = tokens [1]
     species = tokens [2]
     return (paste (substr (genus, 1, 1), species, sep=''))
     } 

   fixed.names = as.character (sapply (names, fix))
   invisible (fixed.names)

} # standardizeSpeciesNames
#------------------------------------------------------------------------------------------------------------------------
renameMatrices = function (matrices, tbl.md)
{
  stopifnot (length (matrices) == nrow (tbl.md))
  names (matrices) = rownames (tbl.md)
  invisible (matrices)

} # renameMatrices
#------------------------------------------------------------------------------------------------------------------------