#'Plot heatmap visualization for time-series data
#'@description Plot heatmap visualization for time-series data
#'@usage PlotHeatMap2(mSetObj=NA, imgName, format="png", dpi=72, width=NA, smplDist="pearson", clstDist="average", colors="bwm", viewOpt="overview", hiRes=FALSE, sortInx = 1, useSigFeature, drawBorder)
#'@param mSetObj Input the name of the created mSetObj (see InitDataObjects)
#'@param imgName Input a name for the plot
#'@param format Select the image format, "png", or "pdf".
#'@param dpi Input the dpi. If the image format is "pdf", users need not define the dpi. For "png" images,
#'the default dpi is 72. It is suggested that for high-resolution images, select a dpi of 300.
#'@param width Input the width, there are 2 default widths, the first, width = NULL, is 10.5.
#'The second default is width = 0, where the width is 7.2. Otherwise users can input their own width.
#'@param smplDist Select distance measure: euclidean, pearson, or minkowski
#'@param clstDist Select clustering algorithm: ward, average, complete, single
#'@param colors Select heatmap colors: bwm, gray
#'@param viewOpt Select overview or detailed view: overview or detail
#'@param hiRes Select high-resolution or not: logical, default set to F
#'@param sortInx Sort by index
#'@param useSigFeature Use significant features only: F or T (default false)
#'@param drawBorder Show cell borders: F or T (default F)
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
#'@export
#'
PlotHeatMap2<-function(mSetObj=NA, imgName, format="png", dpi=72, width=NA, smplDist="pearson",
clstDist="average", colors="bwm", viewOpt="overview", hiRes=FALSE, sortInx = 1,
useSigFeature, drawBorder){
mSetObj <- .get.mSet(mSetObj);
if(mSet$dataSet$facA.type==TRUE){
facA <- as.factor(as.numeric(levels(mSet$dataSet$facA))[mSet$dataSet$facA]);
}else{
facA <- mSetObj$dataSet$facA;
}
if(mSet$dataSet$facB.type==TRUE){
facB <- as.factor(as.numeric(levels(mSet$dataSet$facB))[mSet$dataSet$facB]);
}else{
facB <- mSetObj$dataSet$facB;
}
if(sortInx == 1){
ordInx <- order(facA, facB);
}else{
ordInx <- order(facB, facA);
}
new.facA <- facA[ordInx];
new.facB <- facB[ordInx];
# set up data set. note, need to transpose the data for two way plotting
data <- mSetObj$dataSet$norm[ordInx, ];
# use features from ANOVA2
if(useSigFeature){
hits <- colnames(data) %in% rownames(mSetObj$analSet$aov2$sig.mat);
data <- mSetObj$dataSet$norm[ordInx, hits];
}
hc.dat <- as.matrix(data);
colnames(hc.dat) <- substr(colnames(data), 1, 18) # some names are too long
# set up parameter for heatmap
if(colors=="gbr"){
colors <- colorRampPalette(c("green", "black", "red"), space="rgb")(256);
}else if(colors == "heat"){
colors <- heat.colors(256);
}else if(colors == "topo"){
colors <- topo.colors(256);
}else if(colors == "gray"){
colors <- colorRampPalette(c("grey90", "grey10"), space="rgb")(256);
}else{
colors <- rev(colorRampPalette(RColorBrewer::brewer.pal(10, "RdBu"))(256));
}
if(drawBorder){
border.col<-"grey60";
}else{
border.col <- NA;
}
imgName = paste(imgName, "dpi", dpi, ".", format, sep="");
if(viewOpt == "overview"){
if(is.na(width)){
w <- 9;
}else if(width == 0){
w <- 7.2;
}else{
w <- 7.2;
}
mSetObj$imgSet$htmaptwo <- imgName;
h <- w;
}else{
if(is.na(width)){
minW <- 650;
myW <- nrow(hc.dat)*11 + 150;
if(myW < minW){
myW <- minW;
}
w <- round(myW/72,2);
}else if(width == 0){
w <- 7.2;
}else{
w <- 7.2;
}
if(ncol(hc.dat) >100){
myH <- ncol(hc.dat)*12 + 120;
}else if(ncol(hc.dat) > 50){
myH <- ncol(hc.dat)*12 + 60;
}else{
myH <- ncol(hc.dat)*12 + 20;
}
mSetObj$imgSet$htmaptwo <- imgName;
h <- round(myH/72,2);
}
if(format=="pdf"){
pdf(file = imgName, width=w, height=h, bg="white", onefile=FALSE);
}else{
Cairo::Cairo(file = imgName, unit="in", dpi=dpi, width=w, height=h, type=format, bg="white");
}
annotation <- data.frame(new.facB, new.facA);
colnames(annotation) <- c(mSetObj$dataSet$facB.lbl, mSetObj$dataSet$facA.lbl);
rownames(annotation) <-rownames(hc.dat);
pheatmap::pheatmap(t(hc.dat),
annotation=annotation,
fontsize=8, fontsize_row=8,
clustering_distance_rows = smplDist,
clustering_distance_cols = smplDist,
clustering_method = clstDist,
border_color = border.col,
cluster_rows = T,
cluster_cols = F,
scale = 'row',
color = colors);
dev.off();
mSetObj$analSet$htmap2 <- list(dist.par=smplDist, clust.par=clstDist);
return(.set.mSet(mSetObj));
}
#'Perform ASCA
#'@description The ASCA algorithm was adapted from the ASCA-genes method
#'(analysis of variance (ANOVA) simultaneous component analysis)
#'by Maria Jose Nueda (mj.nueda@ua.es) and Ana Conesa (aconesa@ivia.es)
#'@param mSetObj Input the name of the created mSetObj (see InitDataObjects)
#'@param a specify the number of components for facA
#'@param b specify the number of components for facB
#'@param x specify the number of components for interaction AB
#'@param res specify the number of model residuals
#' type is string, indicating the type of analysis
#' "abc" separately
#' "aab" facA joins with AB
#' "bab" facB joins with AB
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
#'@export
#'
Perform.ASCA <- function(mSetObj=NA, a=1, b=2, x=2, res=2){
mSetObj <- .get.mSet(mSetObj);
X <- mSetObj$dataSet$norm;
Fac = c(a, b, x, res);
# Designa (p x I)
Designa <- model.matrix(~ mSetObj$dataSet$facA-1);
# Designb (p x J)
Designb <- model.matrix(~ mSetObj$dataSet$facB-1);
n<-ncol(X);
p<-nrow(X);
I<-ncol(Designa);
J<-ncol(Designb);
Faca=Fac[1]; # number components Model a
Facb=Fac[2]; # number components Model b
Facab=Fac[3]; # number components Model ab to not consider interations set to 0
Facres=Fac[4]; # number components Residues
# Calculate Overall Mean
offset<-apply(X,2,mean);
# remove the overall mean from the matrix
Xoff <- X-(cbind(matrix(1,nrow=p,ncol=1))%*%rbind(offset));
Model.a<-ASCAfun1(Xoff,Designa,Faca);
Model.b<-ASCAfun1(Xoff,Designb,Facb);
if (Facab != 0 ) {
Model.ab<-ASCAfun2(Xoff,Designa,Designb,Facab);
}
# Collecting models
models <- ls(pattern="Model");
output <- vector(mode="list");
Xres <- Xoff;
for (i in 1: length(models)) {
mymodel <- get(models[i], envir=environment());
output <- c(output, list(mymodel));
Xres <- Xres - mymodel$X;
rm(mymodel);
gc();
}
names(output) <- models
# Model.res
Model.res <- ASCAfun.res(Xres,Facres);
LIST<-list(Model.res);
names(LIST)<-c("Model.res");
mSetObj$analSet$asca <- list(
facNum = Fac,
Xoff = Xoff,
models = c(output,LIST)
);
return(.set.mSet(mSetObj));
}
#'Calculate the Important Variable Cutoff
#'@description This function calculates the all important features based on
#'a specfic cutoff.
#'@usage CalculateImpVarCutoff(mSetObj, spe.thresh, lev.thresh)
#'@param mSetObj Input the name of the created mSetObj (see InitDataObjects)
#'@param spe.thresh alpha threshold, less is better, default less than 5 percentile based chi-square
#'note: spe and leverage are vectors, not a single value, but a list to store the result
#'note: the last model is Model.res, no spe
#'Calculate leverage cutoff based on permutation
#'Calculate the reference distribution of leverages
#'note: leverage.perm is a list with each member in a 3 column matrix
#'@param lev.thresh leverage threshold, the higher better, default more than 95 percentile of permuted leverage
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
#'@export
#'
CalculateImpVarCutoff <- function(mSetObj=NA, spe.thresh = 0.05, lev.thresh = 0.95){
mSetObj <- .get.mSet(mSetObj);
asca.models <- mSetObj$analSet$asca$models;
spe.lims <- lev.lims <- numeric(3);
md.nms <- names(asca.models)[1:3]; # note, last model is Model.res, ignore
names(spe.lims) <- names(lev.lims) <- md.nms;
for (nm in md.nms) { # a, ab, b
model <- asca.models[[nm]];
# get SPE cutoff based on Chiq distribution
m <- mean(model$SPE);
v <- var(model$SPE);
g <- v/(2*m);
h <- 2*m*m/v;
lim <- g*qchisq(1-spe.thresh, df=h);
spe.lims[nm] <- lim;
}
mSetObj$analSet$asca$SPE.cutoff <- spe.lims;
if(is.null(mSetObj$analSet$asca$leverage.perm)){
mSetLev <<- mSetObj
# lev.perm is a list with each 3 col matrix (lvA, lvV, lvAB)
res.perm <- Perform.permutation(20, Get.asca.leverage);
# perm.num may be adjusted on public server
perm.num <- res.perm$perm.num;
lev.perm <- res.perm$perm.res;
rm(mSetLev, pos = ".GlobalEnv")
# convert to a list with 3 members each is a permutation of all variable
# for a single factor
rowNum <- length(lev.perm);
colNum <- ncol (mSetObj$dataSet$norm); # each col is variable
lvA.mat <- lvB.mat <- lvAB.mat <- matrix(0, nrow = rowNum, ncol=colNum);
for(m in 1:length(lev.perm)){
mat <- lev.perm[[m]];
lvA.mat[m,] <- mat[,1]; # facA
lvB.mat[m,] <- mat[,2]; # facB
lvAB.mat[m,] <- mat[,3]; # facAB
}
perm.lv <- list("Model.a"=lvA.mat,
"Model.b"=lvB.mat,
"Model.ab"=lvAB.mat);
mSetObj$analSet$asca$leverage.perm <- perm.lv;
rm(lev.perm);
gc();
}
for (nm in md.nms){
lv.mat <- mSetObj$analSet$asca$leverage.perm[[nm]];
# get the quantile for each var
quant1 <- apply(lv.mat, 2, quantile, probs = lev.thresh);
# get the quantile for each model
lim <- quantile(quant1, probs = lev.thresh);
lev.lims[nm] <- lim;
}
mSetObj$analSet$asca$leverage.cutoff <- lev.lims;
# now get all significant and outlier for each factor based on the threshold
sig.list <- out.list <- list();
for (nm in md.nms){
model <- asca.models[[nm]];
lv <- model$leverage;
spe <- model$SPE;
lv.cutoff <- mSetObj$analSet$asca$leverage.cutoff[nm];
spe.cutoff <- mSetObj$analSet$asca$SPE.cutoff[nm];
lvInx <- lv >= lv.cutoff;
speInx <- spe <= spe.cutoff;
sigInx <- lvInx & speInx;
outInx <- spe > spe.cutoff;
sig.mat <- cbind(lv[sigInx], spe[sigInx]);
colnames(sig.mat) <- c("Leverage", "SPE");
rownames(sig.mat) <- colnames(mSetObj$dataSet$norm)[sigInx];
# order by leverage
ordInx <- order(sig.mat[,1], decreasing=TRUE);
sig.mat <- sig.mat[ordInx,,drop=F];
out.mat <- cbind(lv[outInx], spe[outInx]);
colnames(out.mat) <- c("Leverage", "SPE");
rownames(out.mat) <- colnames(mSetObj$dataSet$norm)[outInx];
# order by SPE
ordInx <- order(out.mat[,2], decreasing=TRUE);
out.mat <- out.mat[ordInx,,drop=F];
# must use double [[, to use dynamical name and assign arbitury list element
sig.list[[nm]] <- sig.mat;
out.list[[nm]]<- out.mat;
nm <- gsub("\\.", "_", nm);
write.csv(sig.mat, file=paste("Sig_features_", nm, ".csv", sep=""));
write.csv(out.mat, file=paste("Outliers_", nm, ".csv", sep=""));
}
mSetObj$analSet$asca$sig.list <- sig.list;
mSetObj$analSet$asca$out.list <- out.list;
return(.set.mSet(mSetObj));
}
#'Function to perform ASCA
#'@description Perform ASCA
#'@param X Numeric, number of compounds
#'@param Design Number of levels in the factor
#'@param Fac Numeric, the factor
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
ASCAfun1<-function (X, Design, Fac) {
n <- ncol(X) # number of genes
I <- ncol(Design) # number of levels in the factor
NK<-NULL
XK<-matrix(NA,nrow=I,ncol=n)
for (i in 1:I) {
sub<-X[Design[,i]==1,]
NK[i]<-nrow(sub)
XK[i,]<-apply(sub,2,mean)
}
NK<-sqrt(NK)
# Weigh the data of the Submodel with the corresponding number of measurement occasions
XKw<- NK*XK
PCA<-PCA.GENES(XKw)
scw<-PCA$scores[,1:Fac]
ld<-PCA$loadings[,1:Fac]
ssq<-PCA$var.exp
if(Fac==1) {
scw<-as.matrix(scw)
ld<-as.matrix(ld)
}
# Re-weigth the scores
sc<-scw/NK
XKrec<-sc%*%t(ld)
Xa<-NULL
TPa<-NULL
for (i in 1:nrow(X)){
position<-which(Design[i,]==1)
Xa<-rbind(Xa,XK[position,])
TPa<-rbind(TPa,XKrec[position,])
}
Ea<-Xa-TPa
# leverage & SPE
leverage<-apply(ld^2,1,sum)
SPE<-apply(Ea^2,2,sum)
output<-list(XK,sc,ld,ssq,Xa,TPa,Ea,leverage,SPE, Fac)
names(output)<-c("data","scores","loadings","var.exp","X","TP","E","leverage","SPE", "facNum");
output
}
#'Function to perform ASCA
#'@description Perform ASCA
#'@param X Numeric, number of compounds
#'@param Desa Number of levels in the factor TIME
#'@param Desb Number of levels in the other factor
#'@param Fac Numeric, the factor
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
ASCAfun2<-function (X, Desa, Desb, Fac) {
n <- ncol(X) # number of genes
I <- ncol(Desa) # number of levels in the factor TIME
J <- ncol(Desb) # number of levels in the other factor
XK1<-matrix(NA,nrow=I,ncol=n);
for (i in 1:I) {
sub<-X[Desa[,i]==1,]
XK1[i,]<-apply(sub,2,mean)
}
XK2<-matrix(NA,nrow=J,ncol=n);
for (j in 1:J) {
sub<-X[Desb[,j]==1,]
XK2[j,]<-apply(sub,2,mean)
}
NK<-matrix(NA,nrow=I,ncol=J)
XK<-matrix(NA,nrow=I*J,ncol=n)
nms.I <- colnames(Desa);
nms.J <- colnames(Desb);
row.nm <- vector(mode="character", length=I*J);
k=1
for (j in 1:J){
for (i in 1:I){
sub<-X[(Desa[,i]+Desb[,j])==2,]
NK[i,j]<-sqrt(nrow(sub))
XK[k,]<-apply(sub,2,mean)-XK1[i,]-XK2[j,];
row.nm[k] <- paste(nms.I[i], nms.J[j], sep=".");
k=k+1
}
}
XKw<-XK*(as.numeric(NK))
rownames(XKw) <- row.nm;
PCA<-PCA.GENES(XKw)
scw<-PCA$scores[,1:Fac]
ld<-PCA$loadings[,1:Fac]
ssq<-PCA$var.exp
if(Fac==1) {
scw<-as.matrix(scw)
ld<-as.matrix(ld)
}
# Re-weigth the scores
sc<-scw/(as.numeric(NK))
XKrec<-sc%*%t(ld)
Xab<-NULL
TPab<-NULL
for (i in 1:nrow(X)){
position1<-which(Desa[i,]==1)
position2<-which(Desb[i,]==1)
Xab<-rbind(Xab,XK[I*(position2-1)+position1,])
TPab<-rbind(TPab,XKrec[I*(position2-1)+position1,])
}
Eab<-Xab-TPab
leverage<-apply(ld^2,1,sum)
SPE<-apply(Eab^2,2,sum)
output<-list(XK,sc,ld,ssq,Xab,TPab,Eab,leverage,SPE, Fac)
names(output)<-c("data","scores","loadings","var.exp","X","TP","E","leverage","SPE", "facNum");
output
}
#'Function to perform ASCA
#'@description Perform ASCA
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'@param X Input list of compounds
#'@param Fac Numeric
#'McGill University, Canada
#'License: GNU GPL (>= 2)
ASCAfun.res <- function(X, Fac){
PCA<-PCA.GENES(X);
sc<-PCA$scores[,1:Fac];
ld<-PCA$loadings[,1:Fac];
ssq<-PCA$var.exp;
if(Fac==1) {
sc<-as.matrix(sc)
ld<-as.matrix(ld)
}
TPres<-sc%*%t(ld)
if(Fac==0){
sc=0
ld=0
TPres<-matrix(0,nrow(X),ncol(X))
}
Eres<-X-TPres
output<-list(sc,ld,ssq,X,TPres,Eres, Fac)
names(output)<-c("scores","loadings","var.exp","X","TP","E", "facNum");
output
}
#'Perform ASCA model validation by permutation
#'@description Perform ASCA model validation by permutation
#'we use Manly's unrestricted permutation of observations
#'which esentially permutes the data over all cells in the
#'designed experiment, then calculates the score for
#'each main factor or interaction components.
#'This will get the null distribution for all effects in one go
#'@usage Perform.ASCA.permute(mSetObj=NA, perm.num)
#'@param mSetObj Input name of the created mSet Object
#'@param perm.num Select the number of permutations, default is 20
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
#'@export
#'
Perform.ASCA.permute<-function(mSetObj=NA, perm.num=20){
# since there are three factors a, b, ab, it is easier
# to permute the data, and let factors fixed, we can pre-calculate
# the needed information for each factor to save computation time
mSetObj <- .get.mSet(mSetObj);
facA <- mSetObj$dataSet$facA;
facB <- mSetObj$dataSet$facB;
facAB <- as.factor(paste(facA, facB, sep=":"));
lvAB <- levels(facAB);
lvAB.mat <- do.call(rbind, strsplit(lvAB, ":"))
lvAB.mat <- cbind(lvAB,lvAB.mat);
# factor size is needed for each iteration
facA.size <- getFactorSize(facA);
facB.size <- getFactorSize(facB);
facAB.size <- getFactorSize(facAB);
# record these information
mSetObj$analSet$asca$perm.info <- list(
facA.size = facA.size,
facB.size = facB.size,
facAB = facAB,
facAB.size = facAB.size,
lvAB.mat = lvAB.mat
);
mSetPerm <<- mSetObj
perm.orig <- Get.asca.tss(dummy, perm=F);
res.perm <- Perform.permutation(perm.num, Get.asca.tss);
rm(mSetPerm, pos = ".GlobalEnv");
gc(); # garbage collection
# perm.num may be adjusted on public server
perm.num <- res.perm$perm.num;
perm.res <- res.perm$perm.res;
# convert to matrix
perm.res <- do.call(rbind, perm.res);
perm.res <- rbind(perm.orig, perm.res);
colnames(perm.res) <- c(mSetObj$dataSet$facA.lbl, mSetObj$dataSet$facB.lbl, "Interaction");
# calculate the significant p value as the proportion of sampled permutations better than or equal to original one
# note, the precision is determined by the permutation number i.e. for 100 time, no better than original
# p value is < 0.01, we can not say it is zero
better.mat <- sweep(perm.res[-1,], 2, perm.res[1,]); # subtract perm.res[1,] for each permutated rows
better.hits <- apply(better.mat>=0, 2, sum);
#p.vec <- better.hits/perm.num;
p.res <- vector(mode="character", length=3);
p.res[better.hits == 0] <- paste("p < ", 1/perm.num, " (1/", perm.num, ")", sep="");
p.res[better.hits > 0] <- paste("p = ", signif(better.hits[better.hits > 0]/perm.num, digits=5), " (", better.hits[better.hits > 0], "/", perm.num, ")", sep="");
## added for test
write.csv(perm.res, file="perm.res.csv");
mSetObj$analSet$asca$perm.p <- p.res;
mSetObj$analSet$asca$perm.mat <- perm.res;
return(.set.mSet(mSetObj));
}
getFactorSize <- function(fac){
lvs <- levels(fac);
size.vec <- numeric(length(lvs));
for(i in length(lvs)){
size.vec[i] <- sum(fac == lvs[i]);
}
size.vec;
}
#'Function for ASCA permutation
#'@description Dummy is used only for the purpose to maintain lapply API
#'this is used for permutation on ANOVA paritions,
#'not on the SCA/PCA part, so the number of selected components is
#'not applicable in this step
#'@param dummy Dummy variable
#'@param perm Logical, TRUE by default
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
#'@export
#'
Get.asca.tss <- function(dummy, perm=T){
X <- mSetPerm$analSet$asca$Xoff;
facA <- mSetPerm$dataSet$facA;
facB <- mSetPerm$dataSet$facB;
facAB <- mSetPerm$analSet$asca$perm.info$facAB;
facA.size <- mSetPerm$analSet$asca$perm.info$facA.size;
facB.size <- mSetPerm$analSet$asca$perm.info$facB.size;
facAB.size <- mSetPerm$analSet$asca$perm.info$facAB.size;
lvAB.mat <- mSetPerm$analSet$asca$perm.info$lvAB.mat;
if(perm){
# permute factor is faster than permute data matrix?
ordInx <- order(runif(length(mSetPerm$dataSet$facA)));
facA <- facA[ordInx];
facB <- facB[ordInx];
facAB <- facAB[ordInx];
}
# BHAN: because of the mean & sd functions have been changed in latest version(R > 3.0.2)
# to calculate column means; replace mean --> colMeans
mnA <- by(X, facA, colMeans);
mnB <- by(X, facB, colMeans);
mnAB <- by(X, facAB, colMeans);
# mnAB should subtract A and B effect
for(i in 1:nrow(lvAB.mat)){
mnAB[[lvAB.mat[i,1]]] <- mnAB[[lvAB.mat[i,1]]] - mnA[[lvAB.mat[i,2]]] - mnB[[lvAB.mat[i,3]]];
}
dist.a <- sum(unlist(lapply(mnA, sumSquare), use.names=F)*(facA.size));
dist.b <- sum(unlist(lapply(mnB, sumSquare), use.names=F)*(facB.size));
dist.ab <- sum(unlist(lapply(mnAB, sumSquare), use.names=F)*(facAB.size));
# return all the distance
c(dist.a, dist.b, dist.ab);
}
# euclidean distance to zero
sumSquare <- function(x){
sum(x*x);
}
Get.asca.leverage <- function(dummy){
X <- mSetLev$analSet$asca$Xoff;
Fac <-mSetLev$analSet$asca$facNum;
Faca=Fac[1]; # number components Model a
Facb=Fac[2]; # number components Model b
Facab=Fac[3]; # number components Model ab to not consider interations set to 0
# permute facA and facB
ordInx <- order(runif(length(mSetLev$dataSet$facA)));
facA <- mSetLev$dataSet$facA[ordInx];
facB <- mSetLev$dataSet$facB[ordInx];
# Designa (p x I)
Desa <- model.matrix(~facA-1);
# Designb (p x J)
Desb <- model.matrix(~facB-1);
n<-ncol(X);
p<-nrow(X);
I <- ncol(Desa)
J <- ncol(Desb)
################################
######### factor A #############
################################
NK1<-numeric(I);
XK1<-matrix(NA,nrow=I,ncol=n);
for (i in 1:I) {
sub<-X[Desa[,i]==1,]
NK1[i]<-nrow(sub)
XK1[i,]<-apply(sub,2,mean)
}
NK1<-sqrt(NK1);
XKw1<- NK1*XK1;
lvA <- Get.Leverage(XKw1, Faca);
# factor B
NK2<-numeric(J);
XK2<-matrix(NA,nrow=J,ncol=n);
for (i in 1:J) {
sub<-X[Desb[,i]==1,]
NK2[i]<-nrow(sub)
XK2[i,]<-apply(sub,2,mean)
}
NK2<-sqrt(NK2);
XKw2<- NK2*XK2;
lvB <- Get.Leverage(XKw2, Facb);
# interaction AB
if (Facab != 0 ) {
NK3<-matrix(NA,nrow=I,ncol=J)
XK3<-matrix(NA,nrow=I*J,ncol=n);
k=1
for (j in 1:J){
for (i in 1:I){
sub<-X[(Desa[,i]+Desb[,j])==2,]
NK3[i,j]<-sqrt(nrow(sub));
XK3[k,]<-apply(sub,2,mean)-XK1[i,]-XK2[j,];
k=k+1
}
}
XKw3<-XK3*(as.numeric(NK3));
lvAB <- Get.Leverage(XKw3, Facab);
}else{
lvAB <- 0;
}
# return all the distance, note, each is a vector
cbind(lvA, lvB, lvAB);
}
#'Fast leverage calculation for permutation purpose
#'@description note, the leverage combines all components
#'the importance feature is for the factor not per components
#'@param XKw Features
#'@param Fac Factor
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
Get.Leverage <- function(XKw, Fac){
PCA <- PCA.GENES(XKw);
ld <- PCA$loadings[,1:Fac];
if(Fac==1) {
ld<-as.matrix(ld);
}
# leverage
apply(ld^2,1,sum);
}
#'Obtain principal components into a matrix that has more variables than individuals
#'@description X is a matrix that has as columns the compounds that were considered as variables in the PCA analysis.
#'First we center the matrix by columns (Xoff) and then we obtain the eigenvalues and the
#'eigenvectors of the matrix Xoff%*%t(Xoff) and we
#'use the equivalences between the loadings and scores to obtain the solution
#'@param X Input matrix that has as columns the compounds that were considered as variables in the PCA analysis
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
PCA.GENES<-function(X){
n<-ncol(X)
p<-nrow(X)
offset<-apply(X,2,mean)
Xoff<-X-(cbind(matrix(1,p,1))%*%rbind(offset));
Xoff <- data.matrix(Xoff);
#eigen command sorts the eigenvalues in decreasing orden.
eigen<-eigen(Xoff%*%t(Xoff)/(p-1))
var<-cbind(eigen$values/sum(eigen$values),cumsum(eigen$values/sum(eigen$values)))
loadings2<-eigen$vectors
scores2<-t(Xoff)%*%loadings2
normas2 <- sqrt(apply(scores2^2,2,sum))
scores1 <- loadings2%*%diag(normas2);
rownames(scores1) <- rownames(X);
loadings1 <- scores2%*%diag(1/normas2)
output <- list(eigen,var,scores1,loadings1)
names(output) <- c("eigen","var.exp","scores","loadings")
output
}
#'Plot scree plots for each model in ASCA
#'@description Plot scree plots for each model in ASCA
#'@usage PlotModelScree(mSetObj, imgName, format="png", dpi=72, width=NA)
#'@param mSetObj Input name of the created mSet Object.
#'@param imgName Input a name for the plot
#'@param format Select the image format, "png", or "pdf".
#'@param dpi Input the dpi. If the image format is "pdf", users need not define the dpi. For "png" images,
#'the default dpi is 72. It is suggested that for high-resolution images, select a dpi of 300.
#'@param width Input the width, there are 2 default widths, the first, width = NULL, is 10.5.
#'The second default is width = 0, where the width is 7.2. Otherwise users can input their own width.
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
#'@export
#'
PlotModelScree <- function(mSetObj=NA, imgName, format="png", dpi=72, width=NA){
mSetObj <- .get.mSet(mSetObj);
imgName = paste(imgName, "dpi", dpi, ".", format, sep="");
if(is.na(width)){
w <- 9;
}else if(width == 0){
w <- 7;
}else{
w <- width;
}
h <- w;
mSetObj$imgSet$asca.scree <- imgName;
Cairo::Cairo(file = imgName, unit="in", dpi=dpi, type=format, width=w, height=h, bg="white");
models <- mSetObj$analSet$asca$models;
# note four plots, model a, b, ab and res
par(mfrow=c(2,2),oma=c(0,0,3,0), cex=1.0)
for ( i in 1: length (models)) {
pc.var <- models[[i]]$var.exp[,1];
if(length(pc.var) > 8){
pc.var <- pc.var[1:8];
}
## Correct Typo: Expalined --> ExPlained by B. Han (17 Sep 2013)
plot(pc.var, type="b", main=paste(names(models)[[i]]),
xlab="Component", ylab="Explained variability", axes=F);
axis(2);
axis(1, at=0:length(pc.var));
box();
}
title("Scree plots of each model", outer=TRUE)
dev.off();
return(.set.mSet(mSetObj));
}
#'Plot score plots of each ASCA model for component 1 against time
#'@description Plot score plots of each ASCA model for component 1 against time
#'@usage PlotASCAModel(mSetObj=NA, imgName, format="png", dpi=72, width=NA, type, colorBW=FALSE)
#'@param mSetObj Input the name of the created mSetObj (see InitDataObjects)
#'@param imgName Input a name for the ASCA score plot
#'@param format Select the image format, "png", or "pdf".
#'@param dpi Input the dpi. If the image format is "pdf", users need not define the dpi. For "png" images,
#'the default dpi is 72. It is suggested that for high-resolution images, select a dpi of 300
#'@param width Input the width, there are 2 default widths, the first, width = NULL, is 10.5.
#'The second default is width = 0, where the width is 7.2. Otherwise users can input their own width.
#'@param type select model a or b
#'@param colorBW Logical, use black/white coloring (T) or not (F)
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
#'@export
#'
PlotASCAModel<-function(mSetObj=NA, imgName, format="png", dpi=72, width=NA, type, colorBW=FALSE){
mSetObj <- .get.mSet(mSetObj);
if(type == "a"){
md <- mSetObj$analSet$asca$models$Model.a;
lbls <- as.character(levels(mSetObj$dataSet$facA));
fac.lbl <- mSetObj$dataSet$facA.lbl;
}else{
md <- mSetObj$analSet$asca$models$Model.b;
lbls <- as.character(levels(mSetObj$dataSet$facB));
fac.lbl <- mSetObj$dataSet$facB.lbl;
}
pcNum <- md$facNum;
# plot at most top 3 PCs
if(pcNum > 3){
pcNum <- 3;
}
imgName = paste(imgName, "dpi", dpi, ".", format, sep="");
h.adj <- ifelse(md$facNum > 1, 5/6, 1)
if(is.na(width)){
w <- ifelse(md$facNum > 1, 6, 5);
}else if(width == 0){
w <- ifelse(md$facNum > 1, 6, 5);
}else{
w <- width;
}
h <- w*h.adj;
if(type == "a"){
mSetObj$imgSet$asca.modelA <- imgName;
}else{
mSetObj$imgSet$asca.modelB <- imgName;
}
Cairo::Cairo(file = imgName, unit="in", dpi=dpi, width=w, height=h, type=format, bg="white");
par(mfrow = c(1, pcNum), cex=1.0);
for(j in 1:pcNum){
## add 'xlab=fac.lbl' & replace ylab as "Scores Component #1 (XX% of variation explained)"
## by B. Han (17 Sep 2013)
# plot(md$scores[,j], type="b", ylab="Concentration / Intensity", col=cols[j], pch=19,
# main=paste(fac.lbl, ", component", j, sep=""), axes=F);
# BHan: modified for Black/White color
# cols[j] --> color
colortype <- ifelse(colorBW, "black", (j + 1));
plot(md$scores[,j], type="b", ylab=paste("Scores (",round(md$var.exp[j,1]*100,2),"% of variation explained)"),
col=colortype, pch=19,
main=paste(fac.lbl, ", component", j, sep=""), axes=F,
xlab=fac.lbl);
axis(2);
axis(1, label=lbls, at=1:length(lbls));
box();
}
dev.off();
return(.set.mSet(mSetObj));
}
#'Plot ASCA interaction plots
#'@description Plot ASCA interaction plots
#'@usage PlotInteraction(mSetObj=NA, imgName, format="png", dpi=72, colorBW=FALSE, width=NA)
#'@param mSetObj Input name of the created mSet Object
#'@param imgName Input a name for the plot
#'@param format Select the image format, "png", or "pdf".
#'@param dpi Input the dpi. If the image format is "pdf", users need not define the dpi. For "png" images,
#'the default dpi is 72. It is suggested that for high-resolution images, select a dpi of 300.
#'@param colorBW Logical, use black and white (TRUE) or colors (FALSE)
#'@param width Input the width, there are 2 default widths, the first, width = NULL, is 10.5.
#'The second default is width = 0, where the width is 7.2. Otherwise users can input their own width.
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
#'@export
#'
PlotInteraction <- function(mSetObj=NA, imgName, format="png", dpi=72, colorBW=FALSE, width=NA){
mSetObj <- .get.mSet(mSetObj);
md <- mSetObj$analSet$asca$models$Model.ab;
ab.lbls <- as.character(levels(mSetObj$dataSet$facA));
ba.lbls <- as.character(levels(mSetObj$dataSet$facB));
pcNum <- md$facNum;
if(pcNum > 3){
pcNum <- 3;
}
imgName = paste(imgName, "dpi", dpi, ".", format, sep="");
if(is.na(width)){
w <- ifelse(md$facNum > 1, 8, 5)
}else if(width == 0){
w <- ifelse(md$facNum > 1, 7, 4.6);
}else{
w <- width;
}
h <- 8;
mSetObj$imgSet$asca.modelAB <- imgName;
Cairo::Cairo(file = imgName, unit="in", dpi=dpi, width=w, height=9, type=format, bg="white");
lmat<-matrix(1:(4*pcNum), nrow=4, byrow=F);
lwid<-rep(4.0, pcNum);
lhei<-rep(c(4.0, 0.4), 2);
layout(lmat, widths = lwid, heights = lhei, respect = FALSE)
# BHan: modified for black/white color; change line type and symbol
numscores <- ifelse( (length(ba.lbls) > length(ab.lbls)), length(ba.lbls), length(ab.lbls) );
if( colorBW ) {
cols <- "black";
linestyle<- (1:numscores) + 1;
pchstyle <- (1:numscores) + 4;
} else {
cols <- (1:numscores) + 1;
linestyle <- 1;
pchstyle <- 19;
}
for(j in 1:pcNum){ # plot layout column wise
scores <- md$scores[,j];
#b y & a x
md.scores <- matrix(scores, nrow=length(levels(mSetObj$dataSet$facB)), byrow=T);
# note, matplot plot each col, need to transpose
par(mar = c(3,4,3,2), cex=1.0);
## replace ylab as "Scores (XX% of variation explained)" by B. Han (17 Sep 2013)
# matplot(t(md.scores), type="b", pch=19, lty=1, axes=F, col = cols,
# ylab="Concentration / Intensity", main=paste("Interaction, component", j, sep=""));
matplot(t(md.scores), type="b", pch=pchstyle, lty=linestyle, axes=F, col = cols,
ylab=paste("Scores (",round(md$var.exp[j,1]*100,2),"% of variation explained)"),
main=paste("Interaction, component", j, sep=""));
axis(1, label=ab.lbls, at=1:length(ab.lbls));
axis(2);
box();
par(mar = c(0,0,0,0), cex=1.0);
plot.new();
# legend("center", horiz=T, legend = as.character(ba.lbls), pch=pchstyle, col=(1:length(ba.lbls))+1, lty=1, bty="n");
legend("center", horiz=T, legend = as.character(ba.lbls), pch=pchstyle, col=cols, lty=linestyle, bty="n");
#b x & a y
op <- par(mar = c(3,4,4,2), cex=1.0);
# cols <- (1:ncol(md.scores)) + 1; # duplicated
## replace ylab as "Scores (XX% of variation explained)" by B. Han (17 Sep 2013)
# matplot(md.scores, type="b", pch=19, lty=1, col= cols, axes=F,
# ylab="Concentration / Intensity", main=paste("Interaction, component", j, sep=""));
matplot(md.scores, type="b", pch=pchstyle, lty=linestyle, col= cols, axes=F,
ylab=paste("Scores (",round(md$var.exp[j,1]*100,2),"% of variation explained)"),
main=paste("Interaction, component", j, sep=""));
axis(1, label=ba.lbls, at=1:length(ba.lbls));
axis(2);
box();
op <- par(mar = c(0,0,0,0), cex=1.0);
plot.new();
# legend("center", horiz=T, legend = as.character(ab.lbls), pch=pchstyle, col=(1:length(ab.lbls))+1, lty=1, bty="n");
legend("center", horiz=T, legend = as.character(ab.lbls), pch=pchstyle, col=cols, lty=linestyle, bty="n");
}
dev.off();
return(.set.mSet(mSetObj));
}
#'Plot the important variables for each factor
#'@description Plot the important variables for each factor
#'@usage PlotAscaImpVar(mSetObj=NA, imgName, format, dpi, width=NA, type)
#'@param mSetObj Input name of the created mSet Object
#'@param imgName Input a name for the plot
#'@param format Select the image format, "png", or "pdf".
#'@param dpi Input the dpi. If the image format is "pdf", users need not define the dpi. For "png" images,
#'the default dpi is 72. It is suggested that for high-resolution images, select a dpi of 300.
#'@param width Input the width, there are 2 default widths, the first, width = NULL, is 10.5.
#'The second default is width = 0, where the width is 7.2. Otherwise users can input their own width.
#'@param type select model a, b, or ab
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
#'@export
#'
PlotAscaImpVar <- function(mSetObj=NA, imgName, format, dpi, width=NA, type){
mSetObj <- .get.mSet(mSetObj);
if(type == "a"){
lvg <- mSetObj$analSet$asca$models$Model.a$leverage;
spe <- mSetObj$analSet$asca$models$Model.a$SPE;
lv.cutoff <- mSetObj$analSet$asca$leverage.cutoff["Model.a"];
spe.cutoff <- mSetObj$analSet$asca$SPE.cutoff["Model.a"];
lbl <- mSetObj$dataSet$facA.lbl;
}else if(type == "b"){
lvg <- mSetObj$analSet$asca$models$Model.b$leverage;
spe <- mSetObj$analSet$asca$models$Model.b$SPE;
lv.cutoff <- mSetObj$analSet$asca$leverage.cutoff["Model.b"];
spe.cutoff <- mSetObj$analSet$asca$SPE.cutoff["Model.b"];
lbl <- mSetObj$dataSet$facB.lbl;
}else{
lvg <- mSetObj$analSet$asca$models$Model.ab$leverage;
spe <- mSetObj$analSet$asca$models$Model.ab$SPE;
lv.cutoff<- mSetObj$analSet$asca$leverage.cutoff["Model.ab"];
spe.cutoff <- mSetObj$analSet$asca$SPE.cutoff["Model.ab"];
lbl <- "Interaction";
}
imgName = paste(imgName, "dpi", dpi, ".", format, sep="");
if(is.na(width)){
w <- 7;
}else if(width == 0){
w <- 7;
}else{
w <- width;
}
h <- w*6/7;
if(type == "a"){
mSetObj$imgSet$asca.impA<-imgName;
}else if(type == "b"){
mSetObj$imgSet$asca.impB<-imgName;
}else{
mSetObj$imgSet$asca.impAB<-imgName;
}
Cairo::Cairo(file = imgName, unit="in", dpi=dpi, width=w, height=h, type=format, bg="white");
PlotSigVar(lvg, spe, lv.cutoff, spe.cutoff,lbl);
dev.off();
return(.set.mSet(mSetObj));
}
#'Supporting function for plotting important variables for each factor
#'@description Supporting function for plotting important variables for each factor
#'note, by control xpd to plot legend outside the plotting area
#'without using layout
#'@param x Input the X variable
#'@param y Input the Y variable
#'@param xline Input the xline
#'@param yline Input the yline
#'@param title Input the title
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
PlotSigVar <- function(x, y, xline, yline, title){
par(mar=c(5,4,3,8), xpd=F);
plot(x, y, xlab="Leverage", ylab="SPE", main=title);
axis.lims <- par("usr"); # x1, x2, y1 ,y2
bad.col <- rgb(0, 0, 1, 0.2);
polygon(c(axis.lims[1], axis.lims[1], axis.lims[2], axis.lims[2]), c(yline, axis.lims[4], axis.lims[4], yline),
col = bad.col, border = NA);
good.col <- rgb(1, 0, 0, 0.2);
polygon(c(xline, xline, axis.lims[2], axis.lims[2]), c(axis.lims[3], yline, yline, axis.lims[3]),
col = good.col, border = NA);
abline(v = xline, col="red");
abline(h = yline, col="red");
# outside x2, and middle y
lgd.x <- axis.lims[2];
lgd.y <- (axis.lims[4] - axis.lims[3])/2;
par(xpd=T);
legend(lgd.x, lgd.y, legend = c("Well-modelled", "Outliers"),
fill=c(good.col, bad.col), bty="n");
box();
}
#'Plot ASCA permutation
#'@description Plot plsda classification performance using different components
#'@usage PlotASCA.Permutation(mSetObj=NA, imgName, format="png", dpi=72, width=NA)
#'@param mSetObj Input name of the created mSet Object
#'@param imgName Input a name for the plot
#'@param format Select the image format, "png", or "pdf".
#'@param dpi Input the dpi. If the image format is "pdf", users need not define the dpi. For "png" images,
#'the default dpi is 72. It is suggested that for high-resolution images, select a dpi of 300.
#'@param width Input the width, there are 2 default widths, the first, width = NULL, is 10.5.
#'The second default is width = 0, where the width is 7.2. Otherwise users can input their own width.
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
#'@export
#'
PlotASCA.Permutation <- function(mSetObj=NA, imgName, format="png", dpi=72, width=NA){
mSetObj <- .get.mSet(mSetObj);
perm.mat<-mSetObj$analSet$asca$perm.mat;
perm.p<-mSetObj$analSet$asca$perm.p;
imgName = paste(imgName, "dpi", dpi, ".", format, sep="");
if(is.na(width)){
w <- 9;
}else if(width == 0){
w <- 9;
}else{
w <- width;
}
h <-2*w;
mSetObj$imgSet$asca.perm <- imgName;
Cairo::Cairo(file = imgName, unit="in", dpi=dpi, width=w, height=h, type=format, bg="white");
par(mfrow=c(3,1), cex=1.0);
nms <- colnames(perm.mat);
for(i in 1:3){
limX <-range(perm.mat[,i])*c(0.9, 1.1);
hst <- hist(perm.mat[-1,i], breaks = "FD", freq=T, xlim = limX, axes=F, bty="l",
ylab="Frequency", xlab= 'Permutation test statistics', col="gray", main=nms[i]);
ticks <- heckbert(limX[1], limX[2], 8);
axis(1, at=ticks);
# add the indicator using original label
h <- max(hst$counts)
arrows(perm.mat[1,i], h/5, perm.mat[1,i], 0, col="red", lwd=2);
text(perm.mat[1,i], h/2, paste('Observed \n statistic \n', perm.p[i]), xpd=T);
}
dev.off();
return(.set.mSet(mSetObj));
}
##############################################
##############################################
########## Utilities for web-server ##########
##############################################
##############################################
#'Table of features well modelled by ASCA
#'@description Table of features well modelled by ASCA
#'@param mSetObj Input the name of the created mSetObj (see InitDataObjects)
#'@param nm Input the name of the well modelled features
#'@author Jeff Xia \email{jeff.xia@mcgill.ca}
#'McGill University, Canada
#'License: GNU GPL (>= 2)
#'@export
GetSigTable.ASCA <- function(mSetObj=NA, nm){
mSetObj <- .get.mSet(mSetObj);
if(nm == "Model.a"){
nmLbl <- paste("main effect", mSetObj$dataSet$facA.lbl);
}else if(nm == "Model.b"){
nmLbl <- paste("main effect", mSetObj$dataSet$facB.lbl);
}else{
nmLbl <- paste("interaction effect between", mSetObj$dataSet$facA.lbl, "and", mSetObj$dataSet$facB.lbl);
}
GetSigTable(mSetObj$analSet$asca$sig.list[[nm]], paste("ASCA. The table shows features that are well modelled by ", nmLbl, ".", sep=""), mSetObj$dataSet$type);
#return(.set.mSet(mSetObj));
}
GetAscaSigMat <- function(mSetObj=NA, type){
mSetObj <- .get.mSet(mSetObj);
if(type == "sigA"){
sig.mat <- CleanNumber(mSetObj$analSet$asca$sig.list[["Model.a"]])
}else if(type == "outA"){
sig.mat <- CleanNumber(mSetObj$analSet$asca$out.list[["Model.a"]])
}else if(type == "sigB"){
sig.mat <- CleanNumber(mSetObj$analSet$asca$sig.list[["Model.b"]]);
}else if(type == "outB"){
sig.mat <- CleanNumber(mSetObj$analSet$asca$out.list[["Model.b"]]);
}else if(type == "sigAB"){
sig.mat <- CleanNumber(mSetObj$analSet$asca$sig.list[["Model.ab"]]);
}else if(type == "outAB"){
sig.mat <- CleanNumber(mSetObj$analSet$asca$out.list[["Model.ab"]])
}else{
print(paste("Unknown data type: ", type));
return(0);
}
fileNm <- paste("asca_",type, ".csv", sep="");
write.csv(signif(sig.mat,5), file=fileNm);
mSetObj$analSet$asca$sig.nm <- fileNm;
if(.on.public.web){
.set.mSet(mSetObj);
return(sig.mat);
}else{
return(.set.mSet(mSetObj));
}
}
GetAscaSigRowNames <- function(mSetObj=NA, type){
mSetObj <- .get.mSet(mSetObj);
if(type == "sigA"){
return(rownames(mSetObj$analSet$asca$sig.list[["Model.a"]]))
}
if(type == "outA"){
return(rownames(mSetObj$analSet$asca$out.list[["Model.a"]]))
}
if(type == "sigB"){
return(rownames(mSetObj$analSet$asca$sig.list[["Model.b"]]))
}
if(type == "outB"){
return(rownames(mSetObj$analSet$asca$out.list[["Model.b"]]))
}
if(type == "sigAB"){
return(rownames(mSetObj$analSet$asca$sig.list[["Model.ab"]]))
}
if(type == "outAB"){
return(rownames(mSetObj$analSet$asca$out.list[["Model.ab"]]))
}
return(0);
}
GetAscaSigColNames <- function(type){
return(c("Leverage", "SPE"));
}
GetAscaSigFileName <- function(mSetObj=NA){
mSetObj <- .get.mSet(mSetObj);
mSetObj$analSet$asca$sig.nm
}
#'Heckbert algorithm
#'@description function to calculate tick mark based on Heckbert algorithm
#'available in the "labeling" package implemented by Justin Talbot
#'adapted from the imagemap package
#'Heckbert's labeling algorithm
#'Heckbert, P. S. (1990) Nice numbers for graph labels, Graphics Gems I, Academic Press Professional, Inc.
#'@param dmin Heckbert
#'@param dmax Heckbert
#'@param m Heckbert
#'@author Justin Talbot \email{jtalbot@@stanford.edu}
heckbert <- function(dmin, dmax, m){
range <- .heckbert.nicenum((dmax-dmin), FALSE)
lstep <- .heckbert.nicenum(range/(m-1), TRUE)
lmin <- floor(dmin/lstep)*lstep
lmax <- ceiling(dmax/lstep)*lstep
seq(lmin, lmax, by=lstep)
}
.heckbert.nicenum <- function(x, round){
e <- floor(log10(x))
f <- x / (10^e)
if(round)
{
if(f < 1.5) nf <- 1
else if(f < 3) nf <- 2
else if(f < 7) nf <- 5
else nf <- 10
}
else
{
if(f <= 1) nf <- 1
else if(f <= 2) nf <- 2
else if(f <= 5) nf <- 5
else nf <- 10
}
nf * (10^e)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.