setGeneric('query', signature='object', function(object, queryString, ignore.case=TRUE)
standardGeneric ('query'))
setGeneric('motifToGene', signature='object', function(object, motifs, source) standardGeneric('motifToGene'))
setGeneric('geneToMotif', signature='object', function(object, geneSymbols, source) standardGeneric('geneToMotif'))
setGeneric('associateTranscriptionFactors', signature='object',
function(object, tbl.withMotifs, source, expand.rows) standardGeneric('associateTranscriptionFactors'))
#------------------------------------------------------------------------------------------------------------------------
setClass ('MotifList',
contains='SimpleList',
representation (elementMetadata='DataFrame',
manuallyCuratedGeneMotifAssociationTable="data.frame"),
prototype(elementType='matrix')
)
setValidity('MotifList', function(object) {
msg = NULL
## what makes for a valid MotifList?
if (length (object) == 0)
return (TRUE)
if (is.null(names(object)))
msg = c(msg, '"names()" must not be NULL')
else if (any(duplicated(names(object))))
msg = c(msg, 'all "names()" must be unique')
if (!all(sapply(object, is.matrix)))
msg = c(msg, 'all matrices must be of class "matrix"')
if (!all(sapply(object, nrow) == 4))
msg = c(msg, 'all matrices must have 4 rows')
# all columns of all matrices should be normalized, summing to one.
# in fact, 2/2086 matrices
# Cparvum-UniPROBE-Cgd2_3490.UP00395
# Pfalciparum-UniPROBE-PF14_0633.UP00394
# fail to pass this test. rounding to one decimal place allows these matrices,
# also, to pass. round (0.98, digits=1) --> 1.0
# see the UniPROBE manpage for the full story on these two matrices
ok = sapply(object, function(elt) all (round (colSums(elt), digits=1) == 1))
if (!all(ok))
msg = c(msg, 'all elements must have colSums equal to 1')
if (is.null(msg)) TRUE else msg
})
#-------------------------------------------------------------------------------
MotifList = function (matrices=list(), tbl.metadata=data.frame ())
{
if (nrow (tbl.metadata) == 0)
df = DataFrame ()
else
df = DataFrame (tbl.metadata, row.names = rownames (tbl.metadata))
tbl.tfClass.filename <- system.file(package="MotifDb", "extdata", "tfClass.tsv")
stopifnot(file.exists(tbl.tfClass.filename))
tbl.tfClass <- read.table(tbl.tfClass.filename, header=TRUE, as.is=TRUE, sep="\t")
object = new ('MotifList', listData=matrices, elementMetadata = df,
manuallyCuratedGeneMotifAssociationTable=tbl.tfClass)
if (length (matrices) == 0)
names (object) = character ()
else
names (object) = rownames (tbl.metadata)
object
} # ctor
#-------------------------------------------------------------------------------
setMethod ('subset', signature = 'MotifList',
function (x, subset, select, drop = FALSE, ...) {
if (missing (subset))
i = TRUE
else {
i = eval(substitute (subset), elementMetadata (x), parent.frame ())
i = try(as.logical(i), silent=TRUE)
if (inherits(i, 'try-error'))
stop('"subset" must be coercible to logical')
i = i & !is.na(i)
} # else
return (x [i])
})
#-------------------------------------------------------------------------------
# transpose 4-row matrices to 4-column, so that there are as many rows as
# there are nucleotides in the motif sequence. meme requires normalized matrices
# exactly as we store them
transformMatrixToMemeRepresentation = function (m)
{
return (t (m))
} # transformMatrixToMemeRepresentation
#-------------------------------------------------------------------------------
# http://stuff.mit.edu/afs/athena/software/meme_v3.5.4/etc/meme-explanation.html
# The motif itself is a position-specific probability matrix giving, for each
# position in the pattern, the observed frequency ('probability') of each
# possible letter. The probability matrix is printed 'sideways'--columns
# correspond to the letters in the alphabet (in the same order as shown in
# the simplified motif) and rows corresponding to the positions of the motif,
# position one first. The motif is preceded by a line starting with
# 'letter-probability matrix:' and containing the length of the alphabet,
# width of the motif, number of occurrences of the motif, and the E-value
# of the motif.
matrixToMemeText = function (matrices)
{
matrix.count = length (matrices)
# incoming matrices have nucleotide rows, position columns. meme
# format, however, requires position-rows, and nucleotide-columns
# calculate the number of lines of text by counting columns
total.transposed.matrix.rows = sum (as.integer (sapply (matrices, ncol)))
predicted.line.count = 12 + (3 * length (matrices)) +
total.transposed.matrix.rows
#s = vector ('character', predicted.line.count)
s = character (predicted.line.count)
s [1] = 'MEME version 4'
s [2] = ''
s [3] = 'ALPHABET= ACGT'
s [4] = ''
s [5] = 'strands: + -'
s [6] = ''
s [7] = 'Background letter frequencies'
s [8] = 'A 0.250 C 0.250 G 0.250 T 0.250 '
s [9] = ''
index = 10
for (name in names (matrices)) {
# transpose the frequency matrix version of the incoming matrix,
# hence 'tfMat'
tfMat = transformMatrixToMemeRepresentation (matrices [[name]])
# meme output may be used by tomtom, which uses matrix names as
# part of image filenames. removed all file-system-unfriendly
# characters here
fixed.name = gsub ('\\/', '_', name)
s [index] = sprintf ('MOTIF %s', fixed.name)
index = index + 1
new.line =
sprintf ('letter-probability matrix: alength= 4 w= %d nsites= %d E=8.1e-020',
nrow (tfMat), 45, 8.1e-020)
s [index] = new.line
index = index + 1
for (r in 1:nrow (tfMat)) {
new.row = sprintf (' %12.10f %12.10f %12.10f %12.10f', tfMat [r,1],
tfMat [r,2], tfMat [r,3], tfMat [r,4])
s [index] = new.row
index = index + 1
}
s [index] = ''
index = index + 1
} # for name
invisible (s)
} # matrixToMemeText
#-------------------------------------------------------------------------------
# connection is a character string, create a file by that name, open the file.
# dispatch to export which dispatches on con='connection'
setMethod ('export', signature=c(object='MotifList', con='character',
format='character'),
function (object, con, format, ...) {
## do minimum work unique to this method, then dispatch to avoid
## code duplication
con = file (con, 'w')
on.exit(close(con))
export(object, con, format, ...)
})
#-------------------------------------------------------------------------------
# write to connection with specified format
# format includes TRANSFAC, meme (also good for tomtom), and tsv
setMethod ('export', signature=c(object='MotifList', con='connection',
format='character'),
function (object, con, format, ...) {
fmt = match.arg (tolower (format), c ('meme', 'transfac','jaspar'))
## match.arg fails if !fmt %in% c('meme', 'transfac'), so no need
## for test
## let the user manage opened cons
if (!isOpen(con)) {
open(con)
on.exit(close(con))
}
if (fmt == 'meme') {
text = matrixToMemeText (object)
} else if (fmt == 'jaspar') {
text = matrixToJasparText (object)
}
cat (text, sep='\n', file=con)
})
#-------------------------------------------------------------------------------
# write to connection, using default format, ??? for matrix list, tsv for
# metadata
setMethod ('export', signature=c(object='MotifList', con='missing',
format='character'),
function (object, con, format, ...) {
fmt = match.arg (tolower (format), c ('meme','jaspar')) # , 'transfac'
if (fmt == 'meme') {
text = paste (matrixToMemeText (object), collapse='\n')
cat (text)
invisible (text)
} else if (fmt == 'jaspar') {
text = paste (matrixToJasparText (object), collapse='\n')
cat (text)
invisible (text)
}
})
#-------------------------------------------------------------------------------
setMethod('show', 'MotifList',
function(object) {
msg = sprintf ('MotifDb object of length %d', length (object))
cat (msg, '\n', sep='')
if (length (object) == 0)
return ()
cat ('| Created from downloaded public sources: 2013-Aug-30', '\n', sep='')
tbl.dataSource = as.data.frame (table (mcols (object)$dataSource))
tbl.org = as.data.frame (table (mcols (object)$organism))
tbl.org = head (tbl.org [order (tbl.org$Freq, decreasing=TRUE),])
totalMatrixCount = length (object)
totalOrganismCount = length (unique (mcols (object)$organism))
dataSourceCount = nrow (tbl.dataSource)
source.singular.or.plural = 'sources'
if (dataSourceCount == 1)
source.singular.or.plural = 'source'
msg = sprintf ('| %d position frequency matrices from %d %s:',
totalMatrixCount, dataSourceCount, source.singular.or.plural)
cat (msg, '\n', sep='')
for (r in 1:nrow (tbl.dataSource)) {
dataSource = tbl.dataSource$Var1 [r]
matrixCount = tbl.dataSource$Freq [r]
msg = sprintf ('| %18s: %4d', dataSource, matrixCount)
cat (msg, '\n', sep='')
} # for r
msg = sprintf ('| %d organism/s', totalOrganismCount)
cat (msg, '\n', sep='')
for (r in 1:nrow (tbl.org)) {
organism = tbl.org$Var1 [r]
orgCount = tbl.org$Freq [r]
msg = sprintf ('| %18s: %4d', organism, orgCount)
cat (msg, '\n', sep='')
} # for r
otherOrgCount = totalMatrixCount - sum (tbl.org$Freq)
if (otherOrgCount > 0) {
msg = sprintf ('| %18s: %4d', 'other', otherOrgCount)
cat (msg, '\n', sep='')
}
if (!is.null (names (object))) {
all.names = names (object)
count = length (all.names)
if (count <= 10)
cat (paste (all.names, '\n'), sep='')
else {
cat (paste (all.names [1:5], '\n'), sep='')
cat ('...', '\n', sep='')
end = length (all.names)
start = end - 4
cat (paste (all.names [start:end], '\n'), sep='')
}
}
})
#-------------------------------------------------------------------------------
setMethod ('query', 'MotifList',
function (object, queryString, ignore.case=TRUE) {
indices = unique (as.integer (unlist (sapply (colnames (mcols (object)),
function (colname)
grep (queryString, mcols (object)[, colname],
ignore.case=ignore.case)))))
object [indices]
})
#-------------------------------------------------------------------------------
# Addition on 2017/06/15 from Matt Richards
# This will not exactly match JASPAR because units are PFM and JASPAR uses PCM
# General JASPAR Format:
# > "Motif Name"\t"Transcription Factor"
# A [ PCMS ]
# C [ PCMS ]
# G [ PCMS ]
# T [ PCMS ]
#
# ...
# Note: the PCMs are space-delimited
matrixToJasparText <- function (matrices)
{
matrix.count <- length (matrices)
# Incoming matrices have nucleotide rows, position columns.
# This is the correct orientation for JASPAR; however, we need to also
# add brackets and letters to them
# Calculate the number of lines of text by counting matrices and assuming
# 6 lines per matrix
predicted.line.count <- 6*matrix.count
#s = vector ('character', predicted.line.count)
s <- character (predicted.line.count)
index <- 1
for (name in names (matrices)) {
# Print the name with an arrow, follwed by the motif
s[index] <- sprintf('>%s',name)
index <- index + 1
# For each line of the matrix, print the correct letter and the
# matrix row surrounded by brackets
motif.matrix <- matrices[name][[1]]
for (r in 1:nrow(motif.matrix)) {
s[index] <- sprintf("%s [ %s ]",
rownames(motif.matrix)[r],
paste(motif.matrix[r,],collapse=" "))
index <- index + 1
}
s[index] <- ""
index <- index + 1
} # for name
# Remove the last blank line
s <- s[-length(s)]
invisible (s)
} # matrixToJasparText
#-------------------------------------------------------------------------------
# returns a data.frame with motif, geneSymbol, source, pubmedID columns
setMethod ('motifToGene', 'MotifList',
function (object, motifs, source) {
source <- tolower(source)
stopifnot(source %in% c("motifdb", "tfclass"))
tbl <- data.frame()
if(source %in% c("motifdb")){
providerId <- NULL # avoid R CMD check note
tbl <- as.data.frame(subset(mcols(object), providerId %in% motifs))
if(nrow(tbl) == 0)
return(data.frame())
tbl <- unique(tbl [, c("geneSymbol", "providerId", "dataSource", "organism", "pubmedID")])
colnames(tbl) <- c("geneSymbol", "motif", "dataSource", "organism", "pubmedID")
tbl <- tbl[, c("motif", "geneSymbol", "dataSource", "organism", "pubmedID")]
tbl$source <- "MotifDb"
}
if(source %in% c("tfclass")){
motif <- NULL
tbl <- subset(object@manuallyCuratedGeneMotifAssociationTable, motif %in% motifs)
if(nrow(tbl) == 0)
return(data.frame())
tbl <- unique(tbl[, c("motif", "tf.gene", "pubmedID")])
tbl <- tbl[order(tbl$motif),]
rownames(tbl) <- NULL
colnames(tbl) <- c("motif", "geneSymbol", "pubmedID")
tbl$source <- "TFClass"
}
tbl
})
#-------------------------------------------------------------------------------
# returns a data.frame with motif, geneSymbol, source, pubmedID columns
setMethod ('geneToMotif', 'MotifList',
function (object, geneSymbols, source) {
source <- tolower(source)
stopifnot(source %in% c("motifdb", "tfclass"))
extract.mdb <- function(gene){
geneSymbol <- NULL # workaround the R CMD check "no visible binding for global variable"
tbl <- as.data.frame(subset(mcols(object), geneSymbol == gene))
tbl <- unique(tbl [, c("geneSymbol", "providerId", "dataSource", "organism", "pubmedID")])
colnames(tbl) <- c("geneSymbol", "motif", "dataSource", "organism", "pubmedID")
tbl
}
if(source %in% c("motifdb")){
tbls <- lapply(geneSymbols, extract.mdb)
result <- do.call(rbind, tbls)
result$source <- "MotifDb"
}
if(source %in% c("tfclass")){
tbl <- subset(object@manuallyCuratedGeneMotifAssociationTable, tf.gene %in% geneSymbols)
tf.gene <- NULL; motif <- NULL # workaround R CMD CHECK "no visible binding ..." bogus error
tbl <- unique(tbl[, c("motif", "tf.gene", "pubmedID")])
tbl <- tbl[order(tbl$tf.gene),]
rownames(tbl) <- NULL
colnames(tbl) <- c("motif", "geneSymbol", "pubmedID")
result <- tbl[, c("geneSymbol", "motif", "pubmedID")]
result$source <- "TFClass"
}
result
})
#-------------------------------------------------------------------------------
setMethod('associateTranscriptionFactors', 'MotifList',
function(object, tbl.withMotifs, source, expand.rows){
source <- tolower(source)
stopifnot(source %in% c("motifdb", "tfclass"))
tbl.out <- data.frame()
if(source %in% c("motifdb")){
# lookup up in the object metadata, expect one TF geneSymbol per matrix name
pfm.ids <- tbl.withMotifs[, "motifName"]
matched.rows <- match(pfm.ids, names(as.list(object)))
#if(length(matched.rows) == nrow(tbl.withMotifs)) {
tbl.new <- mcols(object)[matched.rows, c("geneSymbol", "pubmedID")]
tbl.new$geneSymbol[nchar(tbl.new$geneSymbol)==0] <- NA
tbl.new$pubmedID[nchar(tbl.new$pubmedID)==0] <- NA
tbl.out <- as.data.frame(cbind(tbl.withMotifs, tbl.new))
} # direct
if(source %in% c("tfclass")){
if(! "shortMotif" %in% colnames(tbl.withMotifs)){
stop("MotifDb::assoicateTranscriptionFactors needs a 'shortMotif' column with the TFClass source")
}
tbl.tfClass <- read.table(system.file(package="MotifDb", "extdata", "tfClass.tsv"), sep="\t", as.is=TRUE, header=TRUE)
motif.ids <- tbl.withMotifs[, "shortMotif"]
geneSymbols <- lapply(motif.ids, function(id)
paste(tbl.tfClass$tf.gene[grep(id, tbl.tfClass$motif, fixed=TRUE)], collapse=";"))
geneSymbols <- unlist(geneSymbols)
pubmedIds <- lapply(motif.ids, function(id)
unique(tbl.tfClass$pubmedID[grep(id, tbl.tfClass$motif, fixed=TRUE)]))
pubmedIds <- as.character(pubmedIds)
pubmedIds <- gsub("integer(0)", "", pubmedIds, fixed=TRUE)
tbl.new <- data.frame(geneSymbol=geneSymbols, pubmedID=pubmedIds, stringsAsFactors=FALSE)
tbl.new$geneSymbol[nchar(tbl.new$geneSymbol)==0] <- NA
tbl.new$pubmedID[nchar(tbl.new$pubmedID)==0] <- NA
tbl.out <- as.data.frame(cbind(tbl.withMotifs, tbl.new))
if(expand.rows){
rows.with.na <- which(is.na(tbl.out$geneSymbol))
rows.with.geneSymbol <- setdiff(1:nrow(tbl.out), rows.with.na)
tbl.asIs <- tbl.out[rows.with.na,]
tbl.toExpand <- tbl.out[rows.with.geneSymbol,]
geneSymbols.split <- strsplit(tbl.toExpand$geneSymbol, ";")
counts <- unlist(lapply(geneSymbols.split, length))
geneSymbols.split.vec <- unlist(geneSymbols.split)
tbl.expanded <- splitstackshape::expandRows(tbl.toExpand, counts, count.is.col=FALSE, drop=FALSE)
stopifnot(length(geneSymbols.split.vec) == nrow(tbl.expanded))
tbl.expanded$geneSymbol <- geneSymbols.split.vec
tbl.out <- rbind(tbl.expanded, tbl.asIs)
}
} # indirect
tbl.out
})
#-------------------------------------------------------------------------------
