Commit 4146ad9e authored by damien's avatar damien
Browse files

Added script to compare 1-point and n-point.

parent 6d23722a
# jsonlite is required to read the 1-point POP
tryCatch(library(jsonlite), error=function(e) install.packages('jsonlite'), finally=library(jsonlite))
readOnePoint <- function(path) {
extractJson <- function(csv) {
ret <- read.csv(csv, sep=';', quote="'")
ret$Prob <- sapply(as.character(ret$Prob), jsonlite::fromJSON)
names(ret$Prob) <- ret$Id
as.data.frame(ret)
}
allfiles <- list.files(path, full.names=T)
markers <- sub('.*pedigree-and-probabilities\\.(.*)\\.csv', '\\1', allfiles)
all1p <- lapply(allfiles, extractJson)
names(all1p) <- markers
all1p
}
# from http://stackoverflow.com/a/17289991
read.tcsv = function(file, header=TRUE, sep=",", ...) {
n = max(count.fields(file, sep=sep), na.rm=TRUE)
x = readLines(file)
.splitvar = function(x, sep, n) {
var = unlist(strsplit(x, split=sep))
length(var) = n
return(var)
}
x = do.call(cbind, lapply(x, .splitvar, sep=sep, n=n))
x = apply(x, 1, paste, collapse=sep)
out = read.csv(text=x, sep=sep, header=header, ...)
return(out)
}
readNPoint <- function(path) {
chr <- list.dirs(path, recursive=F, full.names=F)
chrdir <- paste(path, chr, sep='/')
gen <- unique(list.dirs(chrdir, recursive=F, full.names=F))
ret <- list()
for (g in gen) {
for (ch in chr) {
ret[[g]][[ch]] <- lapply(list.files(paste(path, ch, g, sep='/'), full.names=T), function(f) {
x <- read.tcsv(f, sep=';', header=T)
#print(typeof(x))
#n <- x$markers
#rownames(x) <- n
#x$markers <- NULL
#as.numeric(x[-1,])
x
})
}
}
ret
}
onePointByGenAndNum <- function(p.data, gen.name, ind.num) {
i <- which(p.data[[1]]$Gen == gen.name)[ind.num]
tmp <- lapply(p.data, function(row) row[row$Gen == gen.name, 'Prob'][i])
allgeno <- sort(unique(do.call(c, lapply(tmp, function(e) names(e[[1]])))))
ret <- list()
allmark <- sort(names(p.data))
ret <- sapply(allgeno, function(g) sapply(allmark, function(m) { w <- p.data[[m]]$Prob[[i]]; print(w); print(typeof(w)); ifelse(is.null(w[[g]]), 0, w[[g]]) }))
print(typeof(ret))
ret <- as.data.frame(ret)
rownames(ret) <- allmark
ret
}
nPointOnMarkersByGenAndNum <- function(n.data, gen.name, ind.num) {
ret <- lapply(n.data[[gen.name]], function(np) { i <- np[[ind.num]]; i[i$markers != '',] })
ret <- do.call(rbind, ret)
rownames(ret) <- ret$markers
ret$markers <- NULL
ret$locus <- NULL
col.idx <- sort.int(colnames(ret), index.return=T)$ix
ret[,col.idx]
}
read.data <- function(path.base) {
list(op=readOnePoint(paste(path.base, '.1-point', sep='')), np=readNPoint(paste(path.base, '.n-point', sep='')))
}
compare.ind <- function(data, gen, ind) {
xo <- onePointByGenAndNum(data$op, gen, ind)
xn <- nPointOnMarkersByGenAndNum(data$np, gen, ind)
n <- ncol(xn)
r <- floor(sqrt(n))
while (n %% r) {
r <- r - 1
}
par(mfrow=c(r, n / r))
for (a in colnames(x)) {
plot(xo[rownames(x), a], col='blue', type='p', ylim=c(0, 1), xlab=a)
lines(xn[,a], col='red', type='l')
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment