Hello all,

I'm using a relatively new package for time series analysis called SSClust
by Christian Davis and Wenxuan Davis (2006) (can be found at
genemerge.bioteam.com/SSClust.html).

The script gives the error:

> source("SSClust.R")
Error in eval.with.vis(expr, envir, enclos) :
        subscript out of bounds
> traceback()
3: eval.with.vis(expr, envir, enclos)
2: eval.with.vis(ei, envir)
1: source("SSClust.R")

I have included the entire script below.

My data is correctly formatted, and as far as I know, I have all relevant
packages in place to run the propgram.

Can anyone give me a push in the right direction to fix this error?

With thanks

Maurice

##############################################################
# This program implements the SSC algorithm described in     #
# Ma, Castillo-Davis, Zhong, and Liu (2005) Nuc. Acids. Res. #
# "Curve Clustering to Discover Patterns in Time Course Gene #
# Expression Data"                                           #
#                                                            #
# Author: Wenxuan Zhong                                      #
# Email:  wenxuan@stat.harvard.edu                           #
##############################################################

rm(list=ls())
#___________________________________________________________________________________
#parameter specification

# "nchain" is how many chains you want to use when performing the RCEM
# "nclust" is the number of clusters (increment this by one until BIC goes
up)
# "threshold" is the p-value cut-off for the RCEM algorithm
# replace "data.txt" with your data input file with the following
tab-delimited
# format:
#
# genename     time1  time2  time3 ...
# geneA 0.54   0.98   0.22  ...
# geneB 0.11   0.76   0.004 ...
# ...

nchain = 2
nclust = 5
my.data = read.table("/data/Copy of R_data.txt",header=T, row.names=1
,sep="\t")
threshold = 0.1

#end parameter specification
#___________________________________________________________________________________



#_____________________
my.name=rownames(my.data)
my.sample=as.matrix(my.data)[,1:10]
n=dim(my.sample)[1]
m=dim(my.sample)[2]
tt=1:m
library(gss)
source("R_code/ssanova1.R")
source("R_code/mcem.R")
like.rec=bic.rec=NULL
my.clust=list()
clust.centers=list()
for(chain in 1:nchain){
clust=kmeans(my.sample,nclust)$cluster
zeta=NULL
varht=NULL
clust.initial=matrix(nrow=nclust,ncol=m)

for(kk in 1:nclust){
         label=which(clust==kk)
         temp=NULL
         temp.pred=NULL
         temp=ssanova1(clusti~tm,
                random=~1|my.geneid,id.basis=1:m,alpha=1,ext=0,
              data=data.frame(clusti=as.numeric(as.vector(t(my.sample
[label,]))),
                              tm=rep(tt,length(label)),
                              my.geneid=as.factor
(sort(rep(1:length(label),m)))
                             )
              )
         temp.pred=predict(temp,data.frame(tm=tt,offset=rep(0,length(tt))))
         clust.initial[kk,]=temp.pred
         my.ran=mkran(~1|my.geneid,data.frame(my.geneid=as.factor(sort(rep(1:length(label),m))))
)
         random.zb=my.ran$z%*%temp$b
         temp.pred=predict(temp,data.frame(tm=tt))
         zeta[kk]=temp$zeta
         varht[kk]=temp$varht
         clust.initial[kk,]=temp.pred
         rm(temp)
print(kk)   ######
 }

clust.sim=NULL
clust.sim=em.clust(my.sample,tt,clust,nclust,clust.initial
,zeta,varht,thres=threshold,iter.max=20)#
clust.centers[[chain]]=clust.sim$mean.curve
my.zeta<-my.varht<-gamma.prior<-NULL
my.zeta=clust.sim$var.ratio
my.varht=clust.sim$variance
gamma.prior=clust.sim$pp
like.rec[chain]=clust.sim$likelihood
bic.rec[chain]=clust.sim$bic
clust.temp=NULL
for(i in 1:n)
    clust.temp[i]=post(my.sample[i,],nclust,gamma.prior,clust.centers
[[chain]],my.zeta,my.varht)$clust
my.clust[[chain]]=clust.temp
}

####prepare output

output.clust=NULL


output.clust=my.clust[[which(bic.rec==min(bic.rec))[1]]]
output.center=matrix(nrow=nclust,ncol=m)
output.center=clust.centers[[which(bic.rec==min(bic.rec))[1]]]
output.bic=min(bic.rec)

dump(c("output.clust","output.center","output.bic"),"output.Rdata")


source("R_code/show.R")
gene.name.clust(my.name,nclust,output.clust)
write(my.name,"all_genes.txt")
postscript("raw_curves.ps")
mclust.plot(my.sample,nclust,output.clust)
dev.off()


x.select=my.sample
nclass=nclust
cl=output.clust
pp=dim(x.select)[2]
tt=1:pp
UCL=matrix(nrow=nclass,ncol=length(tt))
LCL=matrix(nrow=nclass,ncol=length(tt))
curv=matrix(nrow=nclass,ncol=length(tt))
library(gss)


for(i in 1:nclass)
  {
    class=x.select[cl==i,]
    nn=dim(class)[1]
    y=as.vector(t(class))
    my.time=rep(tt,nn)
    geneid=as.factor(sort(rep(tt,nn)))
    temp=ssanova1(y~my.time,random=~1|geneid,id.basis=tt,alpha=1,ext=0)
    temp.pred=predict(temp,data.frame(my.time=tt),se.fit=T)
    UCL[i,]=temp.pred$fit+1.96*temp.pred$se.fit
    LCL[i,]=temp.pred$fit-1.96*temp.pred$se.fit
    curv[i,]=temp.pred$fit
}

postscript("mean_curves.ps")
par(mfrow=c(2,2))
for(i in 1:nclass)
  {
   plot(tt,curv[i,],type="l",ylab=paste("cluster",i,sep="",collapse=NULL))
   lines(tt,LCL[i,],col=2)
   lines(tt,UCL[i,],col=2)

        }

dev.off()

##### print cluster number and BIC score ######
print(paste("Clusters =", nclust, sep=" "), quote = F);
print(paste("BIC Score =", round(bic.rec,2), sep=" "), quote = F);

	[[alternative HTML version deleted]]

