## Code for reading in Data from Marla Feller's group.
## 2009-07-21
## The layout for EJ's array can be found in ejc_ElectrodeNumbering.jpg
## This file is installed as part of the R package.
## system.file("examples/ejc_ElectrodeNumbering.jpg",package="sjemea")
make.ejc.layout <- function(positions) {
## make the layout for SANGER MEA (cf. make.sanger1.layout)
## This is a hexagonal grid.
xlim <- ylim <- c(-300, 300) #TODO SPOT.
spacing <- 60
cols.an <- toupper(substring(positions,1,2))
columns <- match(cols.an, ejcmealayout$name)
## round columns to nearest integer - i.e. um.
pos <- cbind(x=as.integer(ejcmealayout$x[columns]),
y=as.integer(ejcmealayout$y[columns]),
electrode.num=ejcmealayout$number[columns])
rownames(pos) <- positions
array <- 'EJC1_hex_60um'
layout <- list(xlim=xlim, ylim=ylim, spacing=spacing,
pos=pos, array=array)
class(layout) <- "mealayout"
layout
}
make.ejcmealayout <- function(file=NULL,sep=60,
write.table=FALSE) {
## Create the data file that has the array layout for all electrodes.
## This will not normally need to be called, as the data file is stored
## in the package, and can be accessed using:
## data(ejcmealayout)
##
## FILE: location of the array layout file.
## SEP: electrode separation.
## write.table: set to TRUE if you want the file written to the /tmp
## directory.
## if file is not specified, get the file from the data directory.
if (is.null(file)) {
file <- system.file("data", "ejc_layout.txt", package = "sjemea")
}
table <- read.csv(file, sep='\t')
## convert the row, col, fields into x,y position, assuming a
## hexagonal grid layout; odd-numbered rows need shifting left.
dx <- 2 * cos(pi/6) * sep
x <- table$c * dx
x <- x - ifelse( (table$r %%2) == 1, dx/2, 0)
y <- table$r * sep * cos(pi/3)
table <- cbind(table, x=round(x,2), y=round(y,2))
table
if (write.table) {
write.table(table, file='/tmp/ejcmealayout.txt', sep='\t',
quote=FALSE, row.names=FALSE)
}
table
}
ejcmealayout.check <- function() {
## Check the layout looks okay.
data(ejcmealayout)
with(ejcmealayout, {
plot(x, y, pch=19, cex=.2, asp=1, main='ejcmealayout',
xlab='', ylab='', bty='n')
text(x, y, paste(number, name))
})
}
feller.spiketimes <- function(dir) {
## Read all the spike data contained in DIR, and return as a list.
filenames <- dir(path=dir)
## TODO:
## Assume last three characters contain the electrode position. How
## valid is this? This will be something like "c3a" or "c3b" with
## the first two chars representing the location.
electrodes <- substring(filenames, first=nchar(filenames)-2)
## Check if there any duplicate electrodes; this would indicate that
## there maybe several experiments patched together into same file.
if (any(duplicated(electrodes)))
stop("You have duplicate electrodes in ", dir)
## perhaps best to order electrodes by their number, rather than 2
## char id, as the number varies fairly smoothly with position on the
## array.
## Convert each electrode name in the filename into upper case
## (dropping any trailing a/b/c for multiple cells on a channel) and
## find equivalent channel number.
data(ejcmealayout)
## take first two letters of electrode name to get the AN value - A
## is an alphabetic char [A-G] and N is a number [1-8].
electrode.an <- toupper(substring(electrodes, 1,2))
channel.num <- match( electrode.an, ejcmealayout$name)
channel.num <- ejcmealayout$number[channel.num]
## If we load in the data, sorted by channel number, nearby neurons
## should look most similar.
## When ranking by channel.nu, there will be ties! Need to check that these
## are resolved okay. TODO
rank <- rank(channel.num)
channels <- data.frame(file=filenames, electrode=electrodes,
num=channel.num, rank=rank)
channels.ordered <- channels[order(rank), ]
## channels.ordered is a nice dataframe that shows you what is about
## to be read in.
print(channels.ordered)
## set to 1 if you want to check that the data is read in okay. Else
## divide by 20,000 to get time in seconds.
to.secs <- 1/20000
spikes <- lapply(as.character(channels.ordered$file), function(f) {
file <- paste(dir, f, sep='/')
vals <- scan(file, quiet=TRUE) * to.secs
vals
})
names(spikes) <- channels.ordered$electrode
## todo: could check that last 3 chars of filename match names(spikes)
spikes
}
remove.empty.channels <- function(spikes) {
## Remove any spike trains that are empty, i.e. zero spikes in them.
## This can happen if the beg, end range is too narrow, or if a
## datafile is empty, which happens sometime for the Feller data.
## TODO: currentlly only used by the Feller reader, perhaps it could
## also be used for other routines too?
nspikes <- sapply(spikes, length)
empty <- which(nspikes==0)
if ( any(empty) ) {
spikes <- spikes[-empty]
}
spikes
}
names.to.indexes <- function(names, elems, allow.na=FALSE, allow.regex=TRUE) {
## Return the indexes of where each element of ELEMS is within NAMES.
## If the first element of ELEMS is '-', then return all indexes except
## those matching ELEMS. If ELEMS is NULL, then 1:n is returned, where n is
## the length of NAMES.
## Example:
## names = c('a', 'b', 'c', 'd', 'e')
## names.to.indexes(names, c('d', 'b', 'a')) ## 4 2 1
## names.to.indexes(names, c( '-', 'c', 'a')) ## 2 4 5
## names.to.indexes(names, NULL)
## to check if first element is "-", we have to use this more
## complex expression, as elems[1] == "-" is an error if the first element
## by chance is NA.
if (is.null(elems)) {
return (1:length(names))
}
if ( isTRUE(all.equal("-", elems[1])) ) {
invert = TRUE
elems = elems[-1]
} else {
invert = FALSE
}
indexes = match(elems, names)
if (allow.regex) {
## see if any elements returned NA, in which case try them individually
## as regular expressions.
which.na <- which(is.na(indexes))
if (any(which.na)) {
regex.elems <- elems[which.na]
new.indexes <- lapply(regex.elems, function(r) {grep(r, names)})
new.indexes <- unique(unlist(new.indexes))
indexes <- indexes[-which.na]
indexes <- c(indexes, new.indexes) #TODO, preserve order?
}
allow.na <- TRUE #allow NAs through now.
}
if (!allow.na) {
if (any(is.na(indexes)))
stop('some indexes not found.')
}
if (invert)
indexes = setdiff(1:(length(names)), indexes)
indexes
}
filter.channel.names <- function(spikes, ids) {
## Filter out some channel names.
## Keep only the channels mentioned in IDS.
## If the elements of IDS are numeric, they are assumed to be the
## indexes of the spike trains; otherwise, they are assumed to be the
## names of cells.
## e.g.
## spikes2 <- filter.channel.names(spikes, c('-', 'g4a', 'a6a'))
## spikes2 <- filter.channel.names(spikes, c('g4a', 'a6a'))
## spikes2 <- filter.channel.names(spikes, c(5, 3, 1))
## first call throws away two channels; second call keeps just two channels.
## third just keeps the three trains mentioned.
if (any(is.character(ids)))
ids = names.to.indexes(names(spikes), ids)
spikes[ids]
}
##' Read in a directory of spike times from Marla Feller and create a "spikes"
##' data structure.
##'
##' Marla Feller has provided MEA data from an hexagonal array. The data from
##' each electrode is stored in a separate file. This function reads those text
##' files and creates a "spikes" data structure.
##'
##' The data from each electrode is not in seconds; divide by 20000 to get the
##' time in seconds.
##'
##' @param filename Name of the text file to be read in.
##' @param ids Which electrodes should be kept.
##' @param time.interval How often to estimate the firing rate.
##' @param beg Start time of the recording.
##' @param end End time of the recording.
##' @return Return the data structure 's'.
##' @author Stephen Eglen
##' @seealso \code{\link{jay.read.spikes}} \code{\link{sanger.read.spikes}}
##' @examples
##'
##' data.file <- system.file('examples', '2002-11-26-1', package='sjemea')
##' s <- feller.read.spikes(data.file)
##' ejcmealayout.check()
##' fourplot(s)
##'
##' ## This file contains the picture of the array for comparison
##' system.file("examples/ejc_ElectrodeNumbering.jpg",package="sjemea")
##'
##' ## next two examples show electrodes being rejected/accepted.
##' s <- feller.read.spikes(data.file, ids=c('-', 'a7a', 'b1a', 'h2a'))
##' s <- feller.read.spikes(data.file, ids=c( 'a7a', 'b1a', 'h2a'))
##'
##'
##' @export feller.read.spikes
feller.read.spikes <- function(filename, ids=NULL,
time.interval=1, beg=NULL, end=NULL) {
## Read in data from Marla Feller.
## FILENAME: directory that contains the spike times, one
## channel per file.
## IDS: an optional vector of cell numbers that should be analysed
## -- the other channels are read in but then ignored.
spikes <- feller.spiketimes(filename)
## Now remove spikes outside of range required.
spikes.range <- range(unlist(spikes))
if (!is.null(end)) {
spikes <- lapply(spikes, jay.filter.for.max, max=end)
} else {
end <- spikes.range[2]
}
if (!is.null(beg)) {
spikes <- lapply(spikes, jay.filter.for.min, min=beg)
} else {
beg <- spikes.range[1]
}
## Remove any channels that have zero spikes (which can happen if
## the beg, end range is too narrow, or if a datafile is empty,
## which happens sometime for the Feller data.
spikes <- remove.empty.channels(spikes)
if (!is.null(ids)) {
## Filter out some channel names, either inclusively or exclusively.
spikes <- filter.channel.names(spikes, ids)
}
rec.time <- c(beg, end)
channels <- names(spikes)
## Count the number of spikes per channel, and label them.
nspikes <- sapply(spikes, length)
## meanfiring rate is the number of spikes divided by the (time of
## last spike - time of first spike).
meanfiringrate <- nspikes/ ( end - beg)
rates <- make.spikes.to.frate(spikes, time.interval=time.interval,
beg=beg, end=end)
## Parse the channel names to get the cell positions.
layout <- make.ejc.layout(channels)
## TODO; worry about multiple units recorded at the same location?
unit.offsets <- NULL #default value.
## check that the spikes are monotonic.
check.spikes.monotonic(spikes)
res <- list( channels=channels,
spikes=spikes, nspikes=nspikes, NCells=length(spikes),
meanfiringrate=meanfiringrate,
file=filename,
layout=layout,
rates=rates,
unit.offsets=unit.offsets,
rec.time=rec.time
)
class(res) <- "mm.s"
## feller.breaks = seq(from=0, to=500, by=50)
feller.breaks <- c(0, seq(35, by=70, length=9))
res$corr = corr.index(res, feller.breaks)
res
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.