Nothing
.First.lib <- function(lib, pkg) {
library.dynam("prepRISA", pkg, lib)
}
################################
# prepRISA - 9 March 2007
#
# Jean Thioulouse - Labo Biometrie et Biologie Evolutive, UMR CNRS 5558
# Universite Lyon 1, Batiment Mendel, 43 Boulevard du 11 Novembre 1918
# 69622 Villeurbanne Cedex - France http://pbil.univ-lyon1.fr/JTHome/
# Tel : (33) 4 72 43 27 56 Fax : (33) 4 72 43 13 88
# msn: thioulouse@hotmail.com # jabber: jt69 # .mac: j.thioulouse #
#
################################
"prepRISA" <- function ()
{
#
# user dialogs
#
cat("Data file (no def.) ? ")
filename <- readLines(n=1)
if (filename == "") stop("Filename needed")
cat("Number of groups (def. = 1) ? ")
ngr <- as.integer(readLines(n=1))
if (is.na(ngr)) ngr <- 1
cat("Sum (def. = 1) or Max (2) of peaks ? ")
ioptsm <- as.integer(readLines(n=1))
if (is.na(ioptsm)) ioptsm <- 1
cat("Row percentage tranformation (def. = yes = 1) ? ")
ioptpt <- as.integer(readLines(n=1))
if (is.na(ioptpt)) ioptpt <- 1
cat("Max. number of peaks (def. = All = 0) ? ")
maxnp <- as.integer(readLines(n=1))
if (is.na(maxnp)) maxnp <- 0
cat("Start (def. = mini = 0) ? ")
start <- as.integer(readLines(n=1))
if (is.na(start)) start <- 0
#
# reading input file
#
# tab <- read.table(filename)
tab <- read.table(filename, colClasses=c("factor", "numeric", "numeric", "numeric"))
names(tab) <- c("nech","np","pb","hp")
#
# computing min/max plus a few declarations
#
if (start != 0) mini <- start else mini <- min(tab$pb)
maxi <- max(tab$pb)
claw <- vector(mode="numeric",length=ngr)
ncla <- vector(mode="numeric",length=ngr)
maxl <- rep(maxi, ngr-1)
lim <- vector(mode="numeric",length=ngr+1)
nech <- rep(c(1:nlevels(tab$nech)),table(tab$nech)[unique(tab$nech)])
nlig <- nrow(tab)
pb <- tab$pb
hp <- tab$hp
cat("--- Min. length =", mini, "--- Max. length =", maxi,"---\n")
#
# follow up on user dialogs
#
for (i in 1:ngr) {
cat("Group number",i,":\n")
cat("Class width (def. = 2) ? ")
claw[i] <- as.integer(readLines(n=1))
if (is.na(claw[i])) claw[i] <- 2
}
if (ngr>1) {
for (i in 1:(ngr-1)) {
cat("Group number",i,":\n")
cat("Upper limit ? ")
maxl[i] <- as.integer(readLines(n=1))
}
}
for (i in 1:ngr) {
cat("Group number",i,"class width:", claw[i],"\n")
}
cat("Class limits:\n")
cat("Mini =",mini,"\n")
if (ngr>1) {
for (i in 1:(ngr-1)) {
cat(maxl[i],"\n")
}
}
cat("Maxi =",maxi,"\n")
#
# computing the table of class limits
#
lim[1] <- mini
if (ngr>1) lim[2:ngr] <- maxl
lim[ngr+1] <- maxi
nbech <- length(levels(tab$nech))
cat("Number of samples =",nbech,"\n")
npi <- by(tab$np, tab$nech, max)
for (i in 1:nbech) {
cat("Sample number :",i, "number of peaks =", npi[i],"\n")
}
tablim1 <- mini
i <- 1
while (i <= ngr) {
while (tablim1 <= lim[i+1]) {
ncla[i] <- ncla[i] + 1
tablim1 <- tablim1 + claw[i]
}
i <- i + 1
}
for (i in 1:ngr) {
cat("Group number : ",i,"number of classes =",ncla[i],"\n")
}
nclatot=sum(ncla)
cat("Total number of classes =", nclatot, "\n")
maxnlim <- max(ncla)
cat("Max number of classes =", maxnlim, "\n")
tablim <- matrix(0, nrow=maxnlim+1, ncol=ngr)
tablim[1,] <- lim[1:ngr]
for (i in 1:ngr) {
for (j in 2:(ncla[i]+1)) {
tablim[j, i] <- tablim[j-1, i] + claw[i]
}
tablim[ncla[i]+1, i] <- lim[i+1]
}
#
# Call of the C function to compute the distribution of peaks
#
classe <- function()
.C("repClass",
as.integer(ioptsm),
as.integer(ngr),
as.integer(nbech),
as.integer(nlig),
as.integer(nclatot),
as.integer(maxnlim),
as.integer(pb),
as.integer(ncla),
as.integer(t(tablim)),
as.double(hp),
as.integer(nech),
tabh = double(nbech*nclatot),
PACKAGE = "prepRISA")$tabh
#
# recovering C function output
#
tabh <- matrix(classe(), nrow=nclatot, ncol=nbech)
tabh <- t(tabh)
#
# truncate the number of peaks
#
if ( (maxnp > 0) && ((maxnp < nclatot)) ) {
for (j in 1:nbech) {
listpics <- sort(tabh[j,], decreasing=TRUE)
tabh[j, tabh[j,] < listpics[maxnp]] <- 0
}
}
#
# setting dataframe row names
#
row.names(tabh) <- unique(tab$nech)
tabh <- as.data.frame(tabh)
names(tabh) <- as.character(tablim[1:nclatot])
#
# remove columns of 0
#
tabhsz <- tabh[,apply(tabh,2,sum)!=0]
#
# row percentages
#
if (ioptpt == 1) {
tabh <- tabh/apply(tabh,1,sum)
tabhsz <- tabhsz/apply(tabhsz,1,sum)
}
#
# build output file names
#
extTot <- "_Tot"
extSZ <- "_NoZ"
if (ioptsm == 1) {
extTot <- paste(extTot,"S",sep="")
extSZ <- paste(extSZ,"S",sep="")
} else {
extTot <- paste(extTot,"M",sep="")
extSZ <- paste(extSZ,"M",sep="")
}
if (ioptpt == 1) {
extTot <- paste(extTot,"P",sep="")
extSZ <- paste(extSZ,"P",sep="")
}
#
# take care of filename extension(s)
#
spl1 <- strsplit(filename, ".", fixed=TRUE)[[1]]
if (length(spl1) > 1) {
spl2 <- spl1[2:length(spl1)]
TotOF <- paste(spl1[1],extTot,".",paste(spl2,collapse="."), sep="")
NZOF <- paste(spl1[1],extSZ,".",paste(spl2,collapse="."),sep="")
} else {
TotOF <- paste(spl1[1],extTot,sep="")
NZOF <- paste(spl1[1],extSZ,sep="")
}
#
# write output files
#
write.table(tabh, file=TotOF, quote=FALSE)
cat("Full table :",TotOF,"Number of rows =", nrow(tabh), "Number of columns =",ncol(tabh),"\n")
write.table(tabhsz, file=NZOF, quote=FALSE)
cat("Table without columns of 0 :",NZOF,"Number of rows =", nrow(tabhsz), "Number of columns =",ncol(tabhsz),"\n")
invisible(as.data.frame(tabhsz))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.