inst/scripts/import/stamlab/import.R
bbc543a9
 # stamlab/import.R
 #------------------------------------------------------------------------------------------------------------------------
 library (org.Hs.eg.db)
 library (org.Mm.eg.db)
 #------------------------------------------------------------------------------------------------------------------------
 printf <- function(...) print(noquote(sprintf(...)))
 #------------------------------------------------------------------------------------------------------------------------
72d70625
 run = function (dataDir)
bbc543a9
 {
0074c971
   dataDir <- file.path(dataDir, "stamlab")
bbc543a9
   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',
7f0d65d0
                     pubmedID="22955618")
bbc543a9
     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
4d5f4f4b
   tbl.md$geneIdType = rep (NA_character_, nrow (tbl.md))
bbc543a9
 
   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
 #------------------------------------------------------------------------------------------------------------------------