version 1.1.1

cran · Jan 14, 2011 · 094f868 · 094f868
1 parent b9dd288
commit 094f868
Show file tree

Hide file tree

Showing 35 changed files with 614 additions and 125 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -2,16 +2,16 @@ Package: survJamda
 Type: Package
 Title: Survival prediction by joint analysis of microarray gene
         expression data
-Version: 1.0.1
-Date: 2010-12-05
+Version: 1.1.1
+Date: 2011-01-14
 Author: Haleh Yasrebi
 Maintainer: Haleh Yasrebi<hyasrebi@yahoo.com>
 Description: Microarray gene expression data can be analyzed
         individually or jointly using merging methods or meta-analysis
-        for the survival prediction and risk assessment of patients.
-Depends: survival, survivalROC,ecodist,survJamda.data
+        to predict patients' survival and risk assessment.
+Depends: survival, survivalROC,ecodist,survcomp,survJamda.data
 License: GPL (>= 2)
 LazyLoad: yes
-Packaged: 2010-12-05 21:36:35 UTC; hyasrebi
+Packaged: 2011-01-27 11:47:32 UTC; hyasrebi
 Repository: CRAN
-Date/Publication: 2010-12-06 09:03:32
+Date/Publication: 2011-01-27 12:21:04
diff --git a/R/cal.cox.coef.R b/R/cal.cox.coef.R
@@ -0,0 +1,19 @@
+cal.cox.coef <-
+function (gnExpMat, survivaltime, censor){
+	require (survival)
+
+        cox.coef = NULL
+
+        max.col = ifelse (is.matrix(gnExpMat), ncol(gnExpMat), 1)
+
+        for (i in 1:max.col){
+                if(is.matrix(gnExpMat))
+                        var = gnExpMat[,i]
+                else
+                        var = gnExpMat
+                cox.t = coxph(Surv (survivaltime, censor)~var)          
+                cox.coef = c(cox.coef, cox.t$coef)
+        }
+  return (cox.coef)
+}
+
diff --git a/R/calPerformance.auc.plot.R b/R/calPerformance.auc.plot.R
@@ -4,16 +4,15 @@ function(lst, train.ind, test.ind, file.name,col, method)
 	train = lst$mat[train.ind,]
 	test = lst$mat[test.ind,]
 	options(warn = -1)
-        res = featureselection (train, lst$phyno$surv[train.ind], lst$phyno$censor[train.ind])
+        res = featureselection (train, lst$phyno$surv[train.ind], lst$phyno$censor[train.ind], method, 100)
 
-   	list.p = p.adjust (res$p, method = method)
+	my.func <- featureselection
+        list.p<- do.call(my.func, list(train, lst$phyno$surv[train.ind], lst$phyno$censor[train.ind], method, 100))
 
-	if (method == "none")
-		list.p = order(res$p)[1:100]
-	else
-		list.p = (list.p <= 0.05)
-	lp.train = res$coef[list.p]%*%t(train[,list.p])
-	lp = res$coef[list.p]%*%t(test[,list.p])
+   	cox.coef = cal.cox.coef (train, lst$phyno$surv[train.ind], lst$phyno$censor[train.ind])
+
+	lp.train = cox.coef[list.p]%*%t(train[,list.p])
+	lp = cox.coef[list.p]%*%t(test[,list.p])
 
 	plotROC(test.ind,lst$phyno$surv,lst$phyno$censor, lp, file.name,col)
 }

diff --git a/R/calPerformance.merge.indep.R b/R/calPerformance.merge.indep.R
@@ -1,28 +1,40 @@
 calPerformance.merge.indep <-
-function(lst, train.ind, test.ind, method,gn.nb)
+function(lst, train.ind, test.ind, method,gn.nb,perf.eval)
 {
 	train = lst$mat[train.ind,]
 	test = lst$mat[test.ind,]
 	options(warn = -1)
+
+	my.func <- featureselection
+	list.p<- do.call(my.func, list(train, lst$phyno$surv[train.ind], lst$phyno$censor[train.ind], method, gn.nb))
+       	cox.coef = cal.cox.coef (train, lst$phyno$surv[train.ind], lst$phyno$censor[train.ind])
 
-        res = featureselection (train, lst$phyno$surv[train.ind], lst$phyno$censor[train.ind])
-       	list.p = p.adjust (res$p, method = method)
+	lp.train = cox.coef[list.p]%*%t(train[,list.p])
+	lp = cox.coef[list.p]%*%t(test[,list.p])
 
-	if (method == "none")
-		list.p = order(res$p)[1:gn.nb]
-	else
-		list.p = (list.p <= 0.05)
-	lp.train = res$coef[list.p]%*%t(train[,list.p])
-	lp = res$coef[list.p]%*%t(test[,list.p])
+	lp.train = as.vector(lp.train)
+	lp = as.vector(lp)
 
-	roc.fit =survivalROC (Stime = lst$phyno$surv[test.ind], status =lst$phyno$censor[test.ind], marker=lp, predict.time = mean(lst$phyno$surv[test.ind]), span = 0.25*NROW(test)^(-0.20))
+	if(perf.eval == "auc"){
+		roc.fit =survivalROC (Stime = lst$phyno$surv[test.ind], status =lst$phyno$censor[test.ind], marker=lp, predict.time = mean(lst$phyno$surv[test.ind]), span = 0.25*NROW(test)^(-0.20))
 
-	sgn = ifelse (lp < median(lp.train),0, 1)
-	sgn = as.vector(sgn)
+		sgn = ifelse (lp < median(lp.train),0, 1)
+		sgn = as.vector(sgn)
 
-	cox.hr = coxph(Surv(lst$phyno$surv[test.ind],lst$phyno$censor[test.ind])~sgn)
+		cox.hr = coxph(Surv(lst$phyno$surv[test.ind],lst$phyno$censor[test.ind])~sgn)
 
-	cat ("AUC\tHR(CI)\t\tP-val\n")
-	cat (sprintf("%.2f",roc.fit$AUC), "\t", sprintf("%.2f",summary (cox.hr)[[6]][2]), "(", sprintf("%.2f",summary (cox.hr)[[7]][3]), "-", sprintf("%.2f",summary (cox.hr)[[7]][4]), ")\tp=", summary (cox.hr)[[6]][5],"\n",sep = "")
+		cat ("AUC\tHR(CI)\t\tP-val\n")
+		cat (sprintf("%.2f",roc.fit$AUC), "\t", sprintf("%.2f",summary (cox.hr)[[6]][2]), "(", sprintf("%.2f",summary (cox.hr)[[7]][3]), "-", sprintf("%.2f",summary (cox.hr)[[7]][4]), ")\tp=", summary (cox.hr)[[6]][5],"\n",sep = "")	}
+	else
+		if (perf.eval == "cindex"){
+			res = concordance.index(x=lp,surv.time = lst$phyno$surv[test.ind], surv.event =lst$phyno$censor[test.ind])
+			print(res)
+		}
+		else{
+			dt.tr = data.frame ("time"=lst$phyno$surv[train.ind], "event"=lst$phyno$censor[train.ind], "score" = lp.train)
+			dt.ts = data.frame ("time"=lst$phyno$surv[test.ind], "event"=lst$phyno$censor[test.ind], "score" = lp)
+			res = sbrier.score2proba(dt.tr,dt.ts)
+			print(res)
+		}
 }
 
diff --git a/R/calPerformance.meta.R b/R/calPerformance.meta.R
@@ -24,9 +24,11 @@ function(common.gene,zstat, i, j, geno.files, surv.data, method)
 		surv= surv[!is.na(surv)]
 
 		lp = comb[cm.gn]%*%t(test[,cm.gn])
+		lp = as.vector(lp)
 
 		sgn = ifelse(lp < median(lp), 0,1)
 		sgn = as.vector(sgn)
+
 		cox.hr = coxph(Surv(surv,censor)~sgn)
 
 		roc.fit =survivalROC (Stime = surv, status = censor, marker=lp, predict.time = mean(surv), span = 0.25*NROW(test)^(-0.20))

diff --git a/R/calPerformance.single.indep.R b/R/calPerformance.single.indep.R
@@ -1,29 +1,45 @@
 calPerformance.single.indep <-
-function(lst1, lst2, method,gn.nb)
+function(lst1, lst2, method,gn.nb,perf.eval)
 {
 	train.vec = lst1$phyno
 	test.vec = lst2$phyno
 	options(warn = -1)
 
-	lst = featureselection (lst1$mat, train.vec$surv,train.vec$censor)
+	#list.p = fs (lst1$mat, train.vec$surv, train.vec$censor, method, gn.nb)
+
+	my.func <- featureselection
+	list.p<- do.call(my.func, list(lst1$mat, train.vec$surv, train.vec$censor, method, gn.nb))
 
-	list.p = p.adjust (lst$p, method = method)
+	cox.coef = cal.cox.coef (lst1$mat, train.vec$surv, train.vec$censor)
 
-	if (method == "none")
-		list.p = order(lst$p)[1:gn.nb]
-	else
-		list.p = list.p <= 0.05
-	lp.train = lst$coef[list.p]%*%t(lst1$mat[,list.p])
-	lp = lst$coef[list.p]%*%t(lst2$mat[,list.p])
+	lp.train = cox.coef[list.p]%*%t(lst1$mat[,list.p])
+	lp = cox.coef[list.p]%*%t(lst2$mat[,list.p])
+
+	lp.train = as.vector(lp.train)
+	lp = as.vector(lp)
 
-	roc.fit =survivalROC (Stime = test.vec$surv, status = test.vec$censor, marker=lp, predict.time = mean(test.vec$surv), span = 0.25*NROW(lst2$mat)^(-0.20))
+	if(perf.eval == "auc"){
+		roc.fit =survivalROC (Stime = test.vec$surv, status = test.vec$censor, marker=lp, predict.time = mean(test.vec$surv), span = 0.25*NROW(lst2$mat)^(-0.20))
 
-	sgn = ifelse (lp < median(lp.train),0, 1)
-	sgn = as.vector(sgn)
+		sgn = ifelse (lp < median(lp.train),0, 1)
+		sgn = as.vector(sgn)
 
-	cox.hr = coxph(Surv(test.vec$surv, test.vec$censor)~sgn)
+		cox.hr = coxph(Surv(test.vec$surv, test.vec$censor)~sgn)
 
-	cat ("AUC\tHR(CI)\t\t\tP-val\n")
-        cat (sprintf("%.2f",roc.fit$AUC), "\t", sprintf("%.2f",summary (cox.hr)[[6]][2]), "(", sprintf("%.2f",summary (cox.hr)[[7]][3]), "-", sprintf("%.2f",summary (cox.hr)[[7]][4]),")\t\tp=", summary (cox.hr)[[6]][5],"\n",sep = "")
+		cat ("AUC\tHR(CI)\t\t\tP-val\n")
+        	cat (sprintf("%.2f",roc.fit$AUC), "\t", sprintf("%.2f",summary (cox.hr)[[6]][2]), "(", sprintf("%.2f",summary (cox.hr)[[7]][3]), "-", sprintf("%.2f",summary (cox.hr)[[7]][4]),")\t\tp=", summary (cox.hr)[[6]][5],"\n",sep = "")
+	}
+	else
+		if (perf.eval == "cindex"){
+			res = concordance.index(x=lp,surv.time = test.vec$surv, surv.event =test.vec$censor)
+			print(res)
+		}
+		else{
+			dt.tr = data.frame ("time"=train.vec$surv, "event"=train.vec$censor, "score" = lp.train)
+			dt.ts = data.frame ("time"=test.vec$surv, "event"=test.vec$censor, "score" = lp)
+			res = sbrier.score2proba(dt.tr,dt.ts)
+			print(res)
+		}
+
 }
 
diff --git a/R/cross.val.combat.R b/R/cross.val.combat.R
@@ -43,26 +43,22 @@ function (x,y,censor,batchID, method,gn.nb,plot.roc, ngroup, iter)
 	test.adj = compute.combat("test", "testSample")
 	test.adj = t(test.adj)
 
-	lst = featureselection(train.adj, y[-groups[[j]]],censor[-groups[[j]]])
+	my.func <- featureselection
+	p.list<- do.call(my.func, list(x[-groups[[j]],], y[-groups[[j]]],censor[-groups[[j]]], method, gn.nb))
 
-	p.list <- p.adjust(lst$p,method=method)
-	if (method == "none")
-		p.list = order(p.list)[1:gn.nb]
-	else{
-		p.list = (p.list<= .05)	
-		gn.nb = sum(p.list)
-		cat ("Selected genes nb: ", gn.nb, "\n")
-	}
+	cox.coef = cal.cox.coef (x[-groups[[j]],], y[-groups[[j]]],censor[-groups[[j]]])
 
-  	lp.train = lst$coef[p.list]%*%t(train.adj)[p.list,]
+ 	lp.train = cox.coef[p.list]%*%t(train.adj)[p.list,]
+	lp.train = as.vector(lp.train)
 
 	if (is.vector(test.adj[,p.list]) && length(test.adj[,p.list]) ==length(lst$coef[p.list]))
 		m = test.adj[,p.list]
 	else
 		m = t(test.adj[,p.list])
 
-	lp = lst$coef[p.list]%*%m
-
+	lp = cox.coef[p.list]%*%m
+	lp = as.vector(lp)
+
 	predict.time =  mean(y[groups[[j]]][censor[groups[[j]]]==1])
 
 	roc.fit =survivalROC (Stime = y[groups[[j]]], status = censor[groups[[j]]], marker =lp, predict.time = predict.time, span = 0.25*NROW(test.adj)^(-0.20))

diff --git a/R/cross.val.surv.R b/R/cross.val.surv.R
@@ -27,40 +27,38 @@ function (x, y, censor, ngroup, iter, method, zscore,gn.nb,gn.nb.display,plot.ro
     	sign = NULL	
 	all.fp = NULL
 	all.tp = NULL
-	 options(warn=-1) 
-
+	options(warn=-1) 
+#######################################
    for (j in 1:ngroup) {
 	if (zscore){
 		x[-groups[[j]], ] = scale (t(scale(t(x[-groups[[j]], ]))))
 		x[groups[[j]], ] = scale (t(scale(t(x[groups[[j]], ]))))
 	}
-	lst = featureselection(x[-groups[[j]],], y[-groups[[j]]],censor[-groups[[j]]])
-	p.list <- p.adjust(lst$p,method=method)
-	if (method == "none")
-		p.list = order(p.list)[1:gn.nb]
-	else{
-		p.list = (p.list<= .05)	
-		gn.nb = sum(p.list)
-		cat ("Selected genes nb: ", gn.nb, "\n")
-	}
 
 	if (gn.nb.display >= 1){
 		if (is.matrix(x))
 			print(sort(colnames(x)[p.list][1:gn.nb.display]))
 	}
-
-	lp.train = lst$coef[p.list]%*%t(x[-groups[[j]],p.list])
+
+	my.func <- featureselection
+	p.list<- do.call(my.func, list(x[-groups[[j]],], y[-groups[[j]]],censor[-groups[[j]]], method, gn.nb))
+
+		cox.coef = cal.cox.coef (x[-groups[[j]],], y[-groups[[j]]],censor[-groups[[j]]])
+
+	lp.train = cox.coef[p.list]%*%t(x[-groups[[j]],p.list])
+	lp.train = as.vector(lp.train)
 
 	if (is.vector(x[groups[[j]],p.list]) && length(x[groups[[j]],p.list]) == length(lst$coef[p.list]))
 		m = x[groups[[j]],p.list]
 	else
 		m = t(x[groups[[j]],p.list])
 
-	lp = lst$coef[p.list]%*%m
+	lp = cox.coef[p.list]%*%m
+	lp = as.vector(lp)
 
 	predict.time =  mean(y[groups[[j]]][censor[groups[[j]]]==1])
 
-	roc.fit =survivalROC (Stime = y[groups[[j]]], status = censor[groups[[j]]], marker =lp, predict.time = predict.time, span = 0.25*NROW(x[groups[[j]],order(lst$p)[1:gn.nb]])^(-0.20))
+	roc.fit =survivalROC (Stime = y[groups[[j]]], status = censor[groups[[j]]], marker =lp, predict.time = predict.time, span = 0.25*NROW(x[groups[[j]],])^(-0.20))
 	pred.fit[j] = roc.fit$AUC
 
 	all.fp = rbind(all.fp,roc.fit$FP)
@@ -81,4 +79,3 @@ function (x, y, censor, ngroup, iter, method, zscore,gn.nb,gn.nb.display,plot.ro
 	val = c(mean(pred.fit), summary (cox.hr)[[6]][2])
     return(val)
 }
-
diff --git a/R/eval.merge.simulate.R b/R/eval.merge.simulate.R
@@ -4,7 +4,6 @@ function(d1,d2,tot.genes, gene.nb, zscore)
 	mat = rbind(d1$ds1,d2$ds1)
 
 	cat("\nMerged data set\n")
-	iter.crossval(mat, c(d1$T,d2$T), c(d1$censor,d2$censor),ngroup = 10,zscore =1,
- gn.nb =gene.nb,gn.nb.display = 0)
+	iter.crossval(mat, c(d1$T,d2$T), c(d1$censor,d2$censor),ngroup = 10,zscore =1,gn.nb =gene.nb,gn.nb.display = 0)
 }
 
diff --git a/R/eval.subset.R b/R/eval.subset.R
@@ -12,26 +12,21 @@ function (x, y, censor,iter, method, gn.nb, train.nb)
 	train.ind = lst.samples$train.ind
 	test.ind = lst.samples$test.ind
 
-	lst = featureselection(x[train.ind,], y[train.ind],censor[train.ind])
-
-	p.list <- p.adjust(lst$p,method=method)
-	if (method == "none")
-		p.list = order(p.list)[1:gn.nb]
-	else{
-		p.list = (p.list<= .05)	
-		gn.nb = sum(p.list)
-		cat ("Nb selected genes: ", gn.nb, "\n")
-	}
-
+	my.func <- featureselection
+	p.list<- do.call(my.func, list(x[train.ind,], y[train.ind],censor[train.ind], method, gn.nb))
+
+	cox.coef = cal.cox.coef(x[train.ind,], y[train.ind],censor[train.ind])
 
-	lp.train = lst$coef[p.list]%*%t(x[train.ind,p.list])
+	lp.train = cox.coef[p.list]%*%t(x[train.ind,p.list])
+	lp.train = as.vector(lp.train)
 
 	if (is.vector(x[test.ind,p.list]) && length(x[test.ind,p.list]) == length(lst$coef[p.list]))
 		m = x[test.ind,p.list]
 	else
 		m = t(x[test.ind,p.list])
 
-	lp = lst$coef[p.list]%*%m
+	lp = cox.coef[p.list]%*%m
+	lp = as.vector(lp)
 
 	roc.fit =survivalROC (Stime = as.vector(y[test.ind]), status = as.vector(censor[test.ind]), marker=lp, predict.time = mean(y[test.ind]), span = 0.25*NROW(x[test.ind,])^(-0.20))
 

diff --git a/R/featureselection.R b/R/featureselection.R
@@ -1,5 +1,5 @@
 featureselection <-
-function (gnExpMat, survivaltime, censor){
+function (gnExpMat, survivaltime, censor, method = "none",gn.nb){
 	require (survival)
 
         ploglik = NULL
@@ -18,6 +18,14 @@ function (gnExpMat, survivaltime, censor){
                 ploglik = c (ploglik, pv)
                 cox.coef = c(cox.coef, cox.t$coef)
         }
-  return (list(coef = cox.coef, p = ploglik))
+	ploglik <- p.adjust(ploglik,method=method)
+	if (method == "none")
+		ploglik = order(ploglik)[1:gn.nb]
+	else{
+		ploglik = (ploglik<= .05)	
+		gn.nb = sum(ploglik)
+		cat ("Selected genes nb: ", gn.nb, "\n")
+	}
+  return (ploglik)
 }
 
diff --git a/R/iter.crossval.R b/R/iter.crossval.R
@@ -17,14 +17,14 @@ function (data, surv, censor,ngroup=10, plot.roc = 0, method = "none", zscore=0,
 	cat ("Iteration\tAUC\tHR(CI)\t\tP-val\n")
 	for (i in 1:ngroup){
 		new.lst = cross.val.surv(data, surv, censor,ngroup, i, method, zscore,gn.nb,gn.nb.display,plot.roc)
-res = rbind (res, new.lst)
+		res = rbind (res, new.lst)
 	}
 
-if(ngroup != length(surv)){
-cat ("Avg AUC+/-SD\tHR(CI)\n")
-if (plot.roc)
-legend (0.55,0.1, legend = paste("AUC+/-SD =", sprintf("%.2f",as.numeric(mean(res[,1], na.rm = TRUE))), "+/-", sprintf("%.2f",sd (res[,1],na.rm = TRUE)), sep = " "), bty = "n")
-cat (sprintf("%.2f",as.numeric(mean(res[,1], na.rm = TRUE))),  "+/-", sprintf("%.2f",sd (res[,1],na.rm = TRUE)), "\t", gm(res[,2]), "(", sprintf("%.2f",ci.gm(res[,2])[1]), "-", sprintf("%.2f",ci.gm(res[,2])[2]), ")\n", sep = "")
-}
+	if(ngroup != length(surv)){
+		cat ("Avg AUC+/-SD\tHR(CI)\n")
+		if (plot.roc)
+			legend (0.55,0.1, legend = paste("AUC+/-SD =", sprintf("%.2f",as.numeric(mean(res[,1], na.rm = TRUE))), "+/-", sprintf("%.2f",sd (res[,1],na.rm = TRUE)), sep = " "), bty = "n")
+		cat (sprintf("%.2f",as.numeric(mean(res[,1], na.rm = TRUE))),  "+/-", sprintf("%.2f",sd (res[,1],na.rm = TRUE)), "\t", gm(res[,2]), "(", sprintf("%.2f",ci.gm(res[,2])[1]), "-", sprintf("%.2f",ci.gm(res[,2])[2]), ")\n", sep = "")
+	}
 }
 
diff --git a/R/main.merge.indep.valid.R b/R/main.merge.indep.valid.R
@@ -1,5 +1,5 @@
 main.merge.indep.valid <-
-function(geno.files,surv.data,gn.nb=100,method = "none", normalization= "zscore1")
+function(geno.files,surv.data,gn.nb=100,method = "none", normalization= "zscore1",perf.eval = "auc")
 {
 	require(survival)
         require(survivalROC)
@@ -28,9 +28,9 @@ function(geno.files,surv.data,gn.nb=100,method = "none", normalization= "zscore1
 			lst = prep(common.gene,geno.files,surv.data,x,y)
 
 		if (normalization == "zscore1" || normalization == "combat")
-			splitMerged.indep (geno.files,lst, x, y, method,gn.nb)
+			splitMerged.indep (geno.files,lst, x, y, method,gn.nb,perf.eval)
 		else
-			splitZscore2.merge.indep (common.gene,geno.files,surv.data,lst, x, y, method,gn.nb)	
+			splitZscore2.merge.indep (common.gene,geno.files,surv.data,lst, x, y, method,gn.nb,perf.eval)	
 	}
 }