inst/scripts/import/stamlab/import.R

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

# stamlab/import.R
#------------------------------------------------------------------------------------------------------------------------
library (org.Hs.eg.db)
library (org.Mm.eg.db)
#------------------------------------------------------------------------------------------------------------------------
printf <- function(...) print(noquote(sprintf(...)))
#------------------------------------------------------------------------------------------------------------------------
run = function (dataDir)
{
  dataDir <- file.path(dataDir, "stamlab")
  rawMatrixList <- readRawMatrices (dataDir)
  novels <- readNovelStatus (dataDir)
  matrices <- extractAndNormalizeMatrices (rawMatrixList)
  tbl.md <- createMetadataTable (matrices, novels)
  matrices <- renameMatrices (matrices, tbl.md)

  serializedFile <- "stamlab.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
#------------------------------------------------------------------------------------------------------------------------
readRawMatrices = function (dataDir)
{
  filename <- file.path(dataDir, "de.novo.pwm")
  all.lines = scan (filename, what=character(0), sep='\n', quiet=TRUE)
  title.lines = grep ('UW.Motif', all.lines)
  title.line.count <<- length (title.lines)
  max = title.line.count - 1

  pwms = list ()
  
  for (i in 1:max) {
    start.line = title.lines [i]
    end.line = title.lines [i+1] - 1
    new.pwm = parsePwm (all.lines [start.line:end.line])
    pwms = c (pwms, list (new.pwm))
    } # for i

  
  invisible (pwms)

} # readRawMatrices
#------------------------------------------------------------------------------------------------------------------------
readNovelStatus = function (dataDir)
{
  filename <- file.path(dataDir, "novels.txt")
  novels = scan (filename, what=character(0), sep='\n', quiet=TRUE)
  all.names = sprintf ('UW.Motif.%04d', 1:683)
  status = rep (FALSE, length (all.names))
  true.novels = match (novels, all.names)
  status [true.novels] = TRUE
  names (status) = all.names
  return (status)
  
} # readNovelStatus
#------------------------------------------------------------------------------------------------------------------------
extractAndNormalizeMatrices = function (pwm.list)
{
  ms = sapply (pwm.list, function (element) element$matrix)
  nms = normalizeMatrices (ms)
  names (nms) = sapply (pwm.list, function (element) element$title)
  return (nms)

} # extractAndNormalizeMatrices
#------------------------------------------------------------------------------------------------------------------------
#  matrices = sapply (list.pwms, function (pwm) pwm$matrix)
#  matrix.names = sapply (list.pwms, function (pwm) pwm$title)
#  names (matrices) = matrix.names
convertRawMatricesToStandard = function (tbl.rmat)
{
  matrix.ids = unique (tbl.rmat$id)
  result =  vector ('list', length (matrix.ids))

  i = 1
  for (matrix.id in matrix.ids) {
    tbl.sub = subset (tbl.rmat, id == matrix.id)
      # sanity check this rather unusual representation of a position count matrix
    base.count = as.data.frame (table (tbl.sub$base))
    stopifnot (base.count$Var1 == c ('A', 'C', 'G', 'T'))
      # conservative length check.  actually expect sequence lengths of 6 - 20 bases
    if  (base.count$Freq [1] < 4 && base.count$Freq [1] > 30) {
      printf ('matrix.id %s has sequence of length %d', matrix.id, base.count$Freq [1])
      }
    stopifnot (all (base.count$Freq == base.count$Freq [1]))
    nucleotide.counts = tbl.sub$count
    row.names = c('A', 'C', 'G', 'T'); 
    col.names = 1:(nrow (tbl.sub) / 4)
    m = matrix (nucleotide.counts, byrow=TRUE, nrow=4, dimnames=list(row.names, col.names))
    result [[i]] = m
    i = i + 1
    } # for matrix.id

  names (result) = matrix.ids
  return (result)

} # convertRawMatricesToStandard 
#------------------------------------------------------------------------------------------------------------------------
createAnnotationTable = function ()
{
  tbl.matrix =  read.table ('MATRIX.txt', sep='\t', header=F, as.is=TRUE)
  colnames (tbl.matrix) = c ('id', 'category', 'mID', 'version', 'binder')

  tbl.protein = read.table ('MATRIX_PROTEIN.txt', sep='\t', header=F, as.is=TRUE)
  colnames (tbl.protein) =  c ('id', 'proteinID')

  tbl.species = read.table ('MATRIX_SPECIES.txt', sep='\t', header=F, as.is=TRUE)
  colnames (tbl.species) = c ('id', 'speciesID')

  tbl.anno = read.table ('MATRIX_ANNOTATION.txt', sep='\t', header=F, as.is=TRUE, quote="")
  colnames (tbl.anno) = c ('id', 'attribute', 'value')

  tbl.family  = subset (tbl.anno, attribute=='family') [, -2];   
  colnames (tbl.family) = c ('id', 'family')

  tbl.tax     = subset (tbl.anno, attribute=='tax_group') [,-2]; 
  colnames (tbl.tax) = c ('id', 'tax')

  tbl.class   = subset (tbl.anno, attribute=='class') [,-2];     
  colnames (tbl.class) = c ('id', 'class')

  tbl.comment = subset (tbl.anno, attribute=='comment')[,-2];    
  colnames (tbl.comment) = c ('id', 'comment')

  tbl.pubmed  = subset (tbl.anno, attribute=='medline')[,-2];    
  colnames (tbl.pubmed) = c ('id', 'pubmed')

  tbl.type    = subset (tbl.anno, attribute=='type')[,-2];       
  colnames (tbl.type) = c ('id', 'type')


  tbl.md = merge (tbl.matrix, tbl.species, all.x=TRUE)
  tbl.md = merge (tbl.md, tbl.protein, all.x=TRUE)
  tbl.md = merge (tbl.md, tbl.family, all.x=TRUE)
  tbl.md = merge (tbl.md, tbl.tax, all.x=TRUE)
  tbl.md = merge (tbl.md, tbl.class, all.x=TRUE)
  tbl.md = merge (tbl.md, tbl.pubmed, all.x=TRUE)
  tbl.md = merge (tbl.md, tbl.type, all.x=TRUE)

  fullID = paste (tbl.md$mID, tbl.md$version, sep='.')
  tbl.md = cbind (fullID, tbl.md, stringsAsFactors=FALSE)

  invisible (tbl.md)

} # createAnnotationTable
#------------------------------------------------------------------------------------------------------------------------
# assemble these columns:
#                      names=character(),                    # species-source-gene: stamlab-Hsapiens-UW.Motif.0001
#                      nativeNames=character(),              # UW.Motif.0001
#                      geneSymbols=character(),              # NA
#                      sequenceCounts=integer(),             # NA
#                      organisms=character(),                # Hsapiens
#                      bindingMolecules=character(),         # NA
#                      bindingMoleculeIdTypes=character(),   # NA
#                      bindingDomainTypes=character(),       # NA
#                      dataSources=character(),              # stamlab
#                      experimentTypes=character(),          # digital genomic footprinting
#                      pubmedIDs=character(),                # 22959076
#                      tfFamilies=character())               # NA
#
# from these
#
createMetadataTable = function (matrices, novels)
{
  options (stringsAsFactors=FALSE)
  tbl.md = data.frame ()
  matrix.ids = names (matrices) 
  
  for (matrix.id in matrix.ids) {
    matrix = matrices [[matrix.id]]
    taxon.code = 'Hsapiens'
    geneId.info = NA
    new.row = list (providerName=matrix.id,
                    providerId=matrix.id,
                    dataSource='stamlab',
                    geneSymbol=NA,
                    geneId=NA,
                    geneIdType=NA,
                    proteinId=NA,
                    proteinIdType=NA,
                    organism='Hsapiens',
                    sequenceCount=NA,
                    bindingSequence=NA_character_,
                    bindingDomain=NA,
                    tfFamily=NA,
                    experimentType='digital genomic footprinting',
                    pubmedID="22955618")
    tbl.md = rbind (tbl.md, data.frame (new.row, stringsAsFactors=FALSE))
    full.name = sprintf ('%s-%s-%s', 'Hsapiens', 'stamlab', matrix.id)
    rownames (tbl.md) [nrow (tbl.md)] = full.name
    } # for i

  novelPFM = rep ('knownMotif', nrow (tbl.md))
  novels.ordered = novels [tbl.md$providerName]  # make sure we follow the order in the tbl
  novelPFM [which (novels.ordered)] = 'novelMotif'
  tbl.md$geneId = novelPFM
  tbl.md$geneIdType = rep (NA_character_, nrow (tbl.md))

  invisible (tbl.md)

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

} # renameMatrices
#------------------------------------------------------------------------------------------------------------------------
convertTaxonCode = function (ncbi.code)
{
  #if (is.na (ncbi.code))  
  #  return (NA_character_)
  ncbi.code = as.character (ncbi.code)
  if (ncbi.code %in% c ('-', NA_character_, 'NA'))
    return ('Vertebrata')

  tbl = data.frame (code= c('10090', '10116', '10117', '3702', '3888', '4094', '4102',
                            '4151', '4513', '4565', '4577', '4932', '6239', '7227', '7729',
                            '7742', '8022', '8355', '8364', '9031', '9606', '9913', '9986'),
                    name=c ('Mmusculus', 'Rnorvegicus', 'Rrattus', 'Athaliana', 'Psativum', 
                            'Nsylvestris', 'Phybrida', 'Amajus', 'Hvulgare', 'Taestivam',
                            'Zmays', 'Scerevisiae', 'Celegans', 'Dmelanogaster',
                            'Hroretzi', 'Vertebrata', 'Omykiss', 'Xlaevis', 'Xtropicalis', 
                            'Gallus', 'Hsapiens', 'Btaurus', 'Ocuniculus'),
                    stringsAsFactors=FALSE)

  ncbi.code = as.character (ncbi.code)
  index = which (tbl$code == ncbi.code)
  if (length (index) == 1)
    return (tbl$name [index])
  else {
    write (sprintf (" unable to map organism code |%s|", ncbi.code), stderr ())
    return (NA_character_)
    }

} # convertTaxonCode
#------------------------------------------------------------------------------------------------------------------------
# an empirical and not altogether trustworthy solution to identifying identifier types.
guessProteinIdentifierType = function (moleculeName)
{
  if (nchar (moleculeName) == 0)
    return (NA_character_)
  if (is.na (moleculeName))
    return (NA_character_) 

  first.char = substr (moleculeName, 1, 1)

  if (first.char == 'Y')
    return ('SGD')

  if (first.char %in% c ('P', 'Q', 'O', 'A', 'C'))
    return ('UNIPROT')

  if (length (grep ('^EAW', moleculeName)) == 1)
    return ('NCBI')

  if (length (grep ('^EAX', moleculeName)) == 1)
    return ('NCBI')

  if (length (grep ('^NP_', moleculeName)) == 1)
    return ('REFSEQ')

  if (length (grep ('^BAD', moleculeName)) == 1)
    return ('EMBL')

   return (NA_character_)

} # guessProteinIdentifierType
#------------------------------------------------------------------------------------------------------------------------
normalizeMatrices = function (matrices)
{
  mtx.normalized = sapply (matrices, function (mtx) apply (mtx, 2, function (colvector) colvector / sum (colvector)))
  invisible (mtx.normalized)

} # normalizeMatrices
#------------------------------------------------------------------------------------------------------------------------
assignGeneId = function (proteinId)
{
  if (!exists ('id.uniprot.human')) {

    tbl = toTable (org.Hs.egUNIPROT)
    id.uniprot.human <<- as.character (tbl$gene_id)
    names (id.uniprot.human) <<- tbl$uniprot_id

    tbl = toTable (org.Hs.egREFSEQ)
    tbl = tbl [grep ('^NP_', tbl$accession),]
    id.refseq.human <<- as.character (tbl$gene_id)
    names (id.refseq.human) <<- tbl$accession

    tbl = toTable (org.Mm.egUNIPROT)
    id.uniprot.mouse <<- as.character (tbl$gene_id)
    names (id.uniprot.mouse) <<- tbl$uniprot_id

    tbl = toTable (org.Mm.egREFSEQ)
    tbl = tbl [grep ('^NP_', tbl$accession),]
    id.refseq.mouse <<- as.character (tbl$gene_id)
    names (id.refseq.mouse) <<- tbl$accession
    }

  proteinId = strsplit (proteinId, '\\.')[[1]][1]   # remove any trailing '.*'

  if (proteinId %in% names (id.uniprot.human))
    return (list (geneId=as.character (id.uniprot.human [proteinId]), type='ENTREZ'))

  if (proteinId %in% names (id.uniprot.mouse))
    return (list (geneId=as.character (id.uniprot.mouse [proteinId]), type='ENTREZ'))

  if (proteinId %in% names (id.refseq.human))
    return (list (geneId=as.character (id.refseq.human [proteinId]), type='ENTREZ'))

  if (proteinId %in% names (id.refseq.mouse))
    return (list (geneId=as.character (id.refseq.mouse [proteinId]), type='ENTREZ'))

  found.leading.Y = length (grep ("^Y", proteinId, perl=TRUE))

  if (found.leading.Y == 1)
    return (list (geneId=proteinId, type='SGD'))

  return (list (geneId=NA_character_, type=NA_character_))

} # assignGeneId
#------------------------------------------------------------------------------------------------------------------------
parsePwm = function (text)
{
  #printf ('parsing pwm %s', text [1])
  lines = strsplit (text, '\t')
  title = lines [[1]][1]
  consensus.sequence = lines [[1]][2]
  line.count = length (lines)
  #printf ('%s: %s', title, consensus.sequence)

  result = matrix (nrow=line.count-1, ncol=4, dimnames=list(1:(line.count-1), c ('A','C','G','T')))  
  row = 1
  for (line in lines [2:line.count]) {
    result [row,] = as.numeric (line)
    row = row + 1
    } # for line

  result = t (result)
    
  return (list (title=title, consensus.sequence=consensus.sequence, matrix=result))

} # parsePwm
#------------------------------------------------------------------------------------------------------------------------