.Rhistory

signature %in% master[,2]
library(survival)
### The File takes 218 samples ###############################################################################
#### This file does Prefiltering also based on those genes which discriminate PvsZ and ARE releated to survival
##### This script CLUSTERS Periphery Cells ON THE PvsZ signature and PFS survival  ###########################
library('affy')
library('xlsx')
library('limma')
library('survival')
rm(list =ls())
## load data ##
load("/home/abidata/Dinis/VIP/Main.Data/ExpressionConsole.normalised.RData")
## The phenoData
tab <- read.xlsx(file  = '/home/abidata/Dinis/VIP/Sample.Info/15-07-14 VIP Sample Collection PvsZ DEG.xlsx', sheetIndex =1)
list.patients <- tab[,2]
### Only those patients who have a NON NA is PFS
pheno  <- pData(eset.ec.norm)
pheno.ready <- pheno[!is.na(pheno$PFS),]
list.patients.final <- list.patients[(list.patients %in% pheno.ready[,3])]
## Include ONLY those features that have correspoding annotation ~ 27148 out of 70000
## Include ONLY those Samples which HAVE PFS and PvsZ Annotation
exprs <- eset.ec.norm[featureNames(eset.ec.norm) %in% rownames(anno.GENENAME),sampleNames(eset.ec.norm) %in% list.patients.final]
#############################################################################################################################
######## SIGNATURE WAS IDENTIFIED WITH THOSE FEATURES WHICH DISTINGUISHED PvsZ Clustering AND HAD SURvival Information ######
#############################################################################################################################
## LIMMA PREFILTERING #################################################
pheno <-pData(exprs)
pheno$Case <- as.factor(as.matrix(pheno$Case))
pheno$Topo <- as.factor(as.matrix(pheno$Topo))
pheno$Class <- as.factor(as.matrix(pheno$Class))
mm <- model.matrix( ~ 0 + Topo + Case  + Class , pheno)
cc <- colnames(mm)
cc <-  gsub("_","", cc)
cc <- gsub("-","", cc)
colnames(mm) <- cc
cont.matrix <- makeContrasts(ZvsP = Topocenter - Topoperiphery, levels= mm)
fit <- lmFit(exprs, design = mm)
fit2 <- contrasts.fit(fit, cont.matrix)
fit3 <- eBayes(fit2)
DEG_table <- topTable(fit3, adjust="BH", coef='ZvsP', number = Inf)
list <-  (DEG_table)[DEG_table$adj.P.Val < 0.05,1]
########### Fitting Univariate Cox's Models ######################################################
Xsig.ready <- as.data.frame(t(exprs(exprs[featureNames(exprs) %in% list,])))
######## Getting survival times (PFS) and status #####################
Time <- (as.numeric(as.matrix(pheno$PFS)))
status <- rep(1, dim(pheno)[1]) # 1 because the relapse always occurs
surv.obj <- Surv(Time, status)
##### FITTING A UNIVARIATE COX REGRESSION MODEL ################################
pvalue.sig <- c(0)
pvalue.adj <- c(0)
for ( i in 1:ncol(Xsig.ready)){
q <- unlist(summary(coxph(surv.obj ~ Xsig.ready[,i], data = Xsig.ready)))
pvalue.sig[i] <- q$logtest.pvalue
}
pvalue.adj <- p.adjust(pvalue.sig, method = "fdr")
list.final <- colnames(Xsig.ready)[(pvalue.adj < 0.2)]
signature <- list.final
signature %in% master[,2]
master <-  read.xlsx(file = '/home/abidata/Dinis/VIP/Main.Data/Master_File_FINAL.xlsx', sheetIndex =1)
signature %in% master[,2]
signature
sum((signature %in% master) +0)
master[,2]
##########################################
exprs.final <- exprs[featureNames(exprs) %in% signature,]
annot <- anno.GENENAME[anno.GENENAME[,1] %in% signature,]
featureNames(exprs.final) ==  annot[,1]
###################################################
#### Prepare the data set for clustering ##########
exprs.final
pData(exprs.final)
table(pData(exprs.final)$Case)
exprs.last.final <- exprs.final[,pData(exprs.final)$Topo == 'periphery']
exprs.last.final
exprs.last.final$Topo
table(exprs.last.final$Case)
table(pData(exprs.final)$Case)
cbind(table(pData(exprs.final)$Case),table(exprs.last.final$Case))
hist(exprs.last.final$Case)
hist(table(exprs.last.final$Case))
table(exprs.last.final$Case)
sum((table(exprs.last.final$Case) ==1)+0)
sum((table(exprs.last.final$Case) ==2)+0)
sum((table(exprs.last.final$Case) ==2)+0)
sum((table(exprs.last.final$Case) ==3)+0)
sum((table(exprs.last.final$Case) ==1)+0)
sum((table(exprs.last.final$Case) ==2)+0)
sum((table(exprs.last.final$Case) ==3)+0)
sum((table(exprs.last.final$Case) ==4)+0)
sum((table(exprs.last.final$Case) ==5)+0)
sum((table(exprs.last.final$Case) ==6)+0)
length(table(exprs.last.final$Case))
dim(Y)
library(MASS)
library(mixtools)
library(matrixcalc)
library(stats)
library(affy)
library(MASS)
library(mixtools)
library(matrixcalc)
library(stats)
library(Runuran)
library(truncnorm)
library(Matrix)
library(MCMCpack)
library(psych)
library(VGAM)
library(MixSim)
library(statmod)
library(flexclust)
library(mixAK)
library(mclust)
library(monomvn)
library(cluster)
library(flexmix)
library(survival)
library(utils)
library(rms)
library(pec)
library(ipred)
library(verification)
library(Hmisc)
library(glmpath)
library(glmnet)
library(gplots)
library(doMC)
library(sparcl)
library(NMF)
library(mcfa)
library(kernlab)
library(class)
library(reshape)
library(impute)
library(GGally)
library(vargplvm)
library(vargplvm)
library(lungExpression)
rm(list =ls())
load("/home/bit/ashar/ownCloud/DPMM_RESULTS/ONE_VIEW/Verhark/Final/GeneLevel/GeneLevel.RData")
p1
table(c.final)
logrank <- survdiff(surv.ob ~ c.final)
setwd('/home/bit/ashar/Dropbox/Code/DPmixturemodel/DPplusAFT')
library(MASS)
library(mixtools)
library(matrixcalc)
library(stats)
library(Runuran)
library(truncnorm)
library(Matrix)
library(MCMCpack)
library(psych)
library(VGAM)
library(MixSim)
library(statmod)
library(flexclust)
library(mixAK)
library(mclust)
library(monomvn)
library(cluster)
library(flexmix)
library(survival)
library(utils)
library(rms)
library(pec)
library(ipred)
library(verification)
library(Hmisc)
library(glmpath)
library(glmnet)
library(gplots)
library(doMC)
library(sparcl)
library(NMF)
library(mcfa)
library(kernlab)
library(class)
library(reshape)
library(impute)
library(GGally)
library(xlsx)
setwd('/home/bit/ashar/Dropbox/Code/DPmixturemodel/DPplusAFT')
library(MASS)
library(mixtools)
library(matrixcalc)
library(stats)
library(Runuran)
library(truncnorm)
library(Matrix)
library(MCMCpack)
library(psych)
library(VGAM)
library(MixSim)
library(statmod)
library(flexclust)
library(mixAK)
library(mclust)
library(monomvn)
library(cluster)
library(flexmix)
library(survival)
library(utils)
library(rms)
library(pec)
library(ipred)
library(verification)
library(Hmisc)
library(glmpath)
library(glmnet)
library(gplots)
library(doMC)
library(sparcl)
library(NMF)
library(mcfa)
library(kernlab)
library(class)
library(reshape)
library(impute)
library(GGally)
library(xlsx)
logrank <- survdiff(surv.ob ~ c.final)
logrank
exp(time[c==1])
mean(exp(time[c==1]))
mean(exp(time[c==2]))
mean(exp(time[c==3]))
mean(exp(time[c==4]))
table(c.final.new)
mean(exp(time.new[c.new.final==4]))
mean(exp(time.new[c.final.new==4]))
mean(exp(time.new[c.final.new==1]))
mean(exp(time.new[c.final.new==3]))
mean(exp(time.new[c.final.new==4]))
mean(exp(time.new[c.final.new==2]))
Nps
rank <- matrix(0, nrow = N, ncol =Nps)
for (j in 1:Nps){
Y.scaled <- matrix(0, nrow = N, ncol =D)
for ( v in 1:2){
clust <- which(c.list[[j]] == v)
Y.scaled[clust,1:D] <- scale(Y[clust,1:D], center = TRUE, scale = TRUE)
}
for ( i in 1:N){
rank[i,j] <- dMVN(as.vector(t(Y[i,1:D])), mean = mu.list[[j]][c.final[i],1:D], Q= S.list[[j]][c.final[i],1:D,1:D], log = TRUE) - dMVN(as.vector(t(Y[i,1:D])), mean = mu.list[[j]][1,1:D], Q= S.list[[j]][1,1:D,1:D], log = TRUE) +  dnorm(x = That[i], mean = beta0.list[[j]][c.final[i]] + betahat.list[[j]][c.final[i],1:D] %*% as.vector(t(Y.scaled[i,1:D])), sd = sqrt(sigma2.list[[j]][c.final[i]]), log =TRUE) -  dnorm(x = That[i], mean = beta0.list[[j]][1] + betahat.list[[j]][1,1:D] %*% as.vector(t(Y.scaled[i,1:D])), sd = sqrt(sigma2.list[[j]][1]), log =TRUE)
}
}
avg.rank <- apply(rank,1,mean)
order.zo <- range01(avg.rank)
order.train <- sort(order.zo,index.return = TRUE, decreasing = TRUE)
Y.order <- Y[order.train$ix,]
c.final.order <- c.final[order.train$ix]
c.final.order
=c("violet","blue","green","red")[c.final.order]
c("violet","blue","green","red")
c.final.order
c("violet","blue","green","red")[c.final.order]
hmcols<-colorRampPalette(c("green","black","red"))(256)
heatmap.2(x = t(Y.order), Rowv =TRUE ,Colv = FALSE, dendrogram = c("row"), scale ="row",ColSideColors =c("green","red","blue","orange")[c.final.order], labRow = colnames(Y.vijver.train), labCol = NA, main = ' \n Training Set \n 58 Gene DPMM signature', col =hmcols, cexRow =0.40, trace ="none", key.title = "Color Key")
heatmap.2(x = t(Y.order), Rowv =TRUE ,Colv = FALSE, dendrogram = c("row"), scale ="row",ColSideColors =c("green","red","blue","orange")[c.final.order], labRow = colnames(Y), labCol = NA, main = ' \n Training Set \n 49 Gene DPMM signature', col =hmcols, cexRow =0.40, trace ="none", key.title = "Color Key")
c.final.order
which(c.final.order==1)
which(c.final.order==3)
which(c.final.order==4)
which(c.final.order==1)
c(which(c.final.order==1),which(c.final.order==3),which(c.final.order==4),which(c.final.order==2))
order.2 <- c(which(c.final.order==1),which(c.final.order==3),which(c.final.order==4),which(c.final.order==2))
order.2 <- c(which(c.final.order==1),which(c.final.order==3),which(c.final.order==4),which(c.final.order==2))
Y.order.2 <- Y.order[order.2,]
c.final.order.2 <- c.final.order[order.2]
c.final.order.2
colnames(Y.order)
colnames(Y.order.2)
hmcols<-colorRampPalette(c("green","black","red"))(256)
heatmap.2(x = t(Y.order.2), Rowv =TRUE ,Colv = FALSE, dendrogram = c("row"), scale ="row",ColSideColors =c("green","red","blue","orange")[c.final.order.2], labRow = colnames(Y.order.2), labCol = NA, main = ' \n Training Set \n 49 Gene DPMM signature', col =hmcols, cexRow =0.40, trace ="none", key.title = "Color Key")
legend("topright", legend = c("Best","Good Moderate","Bad moderate","Worst"),fill = c("Green","Blue","Orange","Red"), cex = 0.4)
p5
surv.ob <- Surv(exp(time),censoring)
Classes <- c.final
Classes <- as.factor(Classes)
surv.fit <- survfit(surv.ob ~ Classes)
p5 <- ggsurv(surv.fit, main = " DPMM \n Kaplan Meier Estimators", surv.col = c("green","red","blue","orange"), cens.col ="Blue")
p5
D
pdf('DPMMSignature.pdf')
hmcols<-colorRampPalette(c("green","black","red"))(256)
heatmap.2(x = t(Y.order.2), Rowv =TRUE ,Colv = FALSE, dendrogram = c("row"), scale ="row",ColSideColors =c("green","red","blue","orange")[c.final.order.2], labRow = colnames(Y.order.2), labCol = NA, main = ' \n Training Set \n 47 Gene DPMM signature', col =hmcols, cexRow =0.40, trace ="none", key.title = "Color Key")
legend("topright", legend = c("Best","Good Moderate","Bad moderate","Worst"),fill = c("Green","Blue","Orange","Red"), cex = 0.4)
p5
dev.off()
pc <- prcomp(Y)
pc.pred <- predict(pc,newdata = Y)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], colour= c("green","red","blue","orange")[as.factor(c.final)])) + ggtitle(" DPMM Clustering") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour = "Classes")
p1
pc <- prcomp(Y)
pc.pred <- predict(pc,newdata = Y)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2]])) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour = "Classes")
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2]]) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour = "Classes")
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2]])) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour = "Classes")
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], colour= c("green","red","blue","orange")[as.factor(c.final)])) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour = "Classes")
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2]) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour = "Classes")
)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2])) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour = "Classes")
p1
pc <- prcomp(Y)
pc.pred <- predict(pc,newdata = Y)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2])) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(aes(colour=c("green","red","blue","orange")[c.final.order.2]),shape=19) + labs(y = "PC1", x = "PC2", colour = "Classes")
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2])) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(aes(colour=c("green","red","blue","orange")[c.final.order.2]),shape=19) + labs(y = "PC1", x = "PC2")
p1
c("green","red","blue","orange")[c.final]
pc <- prcomp(Y)
pc.pred <- predict(pc,newdata = Y)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2])) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(aes(colour=c("green","red","blue","orange")[c.final]),shape=19) + labs(y = "PC1", x = "PC2", colour =  c("Worst","Best","Moderate Good","Moderate Bad"))
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2])) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour =  c("Worst","Best","Moderate Good","Moderate Bad"))
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], fill = c("green","red","blue","orange")[c.final] )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour =  c("Worst","Best","Moderate Good","Moderate Bad"))
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], fill = c.final )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour =  c("Worst","Best","Moderate Good","Moderate Bad"))
p1
c.final
as.factor(c.final)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], fill = as.factor(c.final) )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2", colour =  c("Worst","Best","Moderate Good","Moderate Bad"))
p1
colour=c.final
c.final
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2] )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(aes(colour=c.final),shape=19) + labs(y = "PC1", x = "PC2", colour =  c("Worst","Best","Moderate Good","Moderate Bad"))
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2] )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(aes(colour=as.factor(c.final)),shape=19) + labs(y = "PC1", x = "PC2", colour =  c("Worst","Best","Moderate Good","Moderate Bad"))
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2] )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(aes(colour=c("green","red","blue","orange")[as.factor(c.final)]),shape=19) + labs(y = "PC1", x = "PC2", colour =  c("Worst","Best","Moderate Good","Moderate Bad"))
p1
cbPalette <- c("green","red","blue","orange")
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2] )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(aes(colour=as.factor(c.final)),shape=19) + labs(y = "PC1", x = "PC2", colour =  c("Worst","Best","Moderate Good","Moderate Bad"))
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2] )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(aes(colour=as.factor(c.final)),shape=19) + labs(y = "PC1", x = "PC2", colour =  c("Worst","Best","Moderate Good","Moderate Bad")) + scale_fill_manual(values=cbPalette)
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2] )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(aes(colour=as.factor(c.final)),shape=19) + labs(y = "PC1", x = "PC2") + scale_fill_manual(values=cbPalette)
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2] )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_fill_manual(values=cbPalette)
p1
c.final
ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], fill= c.final )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak")
geom_point(shape=19)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], fill= c.final )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2")
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], fill= as.factor(c.final) )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_fill_manual(breaks = c("1", "2", "3","4"),values=cbPalette)
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], fill= as.factor(c.final) )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_fill_manual(breaks = c("1", "2", "3","4"),values=cbPalette)
p1
cbPalette <- c("green","red","blue","orange")
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], fill= as.factor(c.final) )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_fill_manual(breaks = c("1", "2", "3","4"),values=cbPalette)
p1
scale_fill_manual(breaks = c("1", "2", "3","4"),values=  c("green","red","blue","orange"))
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], fill= as.factor(c.final) )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_fill_manual(breaks = c("1", "2", "3","4"),values=  c("green","red","blue","orange"))
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], color = as.factor(c.final) )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_color_manual(breaks = c("1", "2", "3","4"),values=  c("green","red","blue","orange"))
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], color = as.factor(c.final)] )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_color_manual(breaks = c("1", "2", "3","4"),values=  c("green","red","blue","orange"))
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], color = as.factor(c.final) )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_color_manual(breaks = c("1", "2", "3","4"),values=  c("green","red","blue","orange"))
p1
pc$sd[1]
(pc$sd[1]+ pc$sd[2])/sum(pc$sd)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], color = as.factor(c.final) )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak \n 18% Variance Explained") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_color_manual(breaks = c("1", "2", "3","4"),values=  c("green","red","blue","orange"))
p1
pdf('DPMMSignature.pdf')
hmcols<-colorRampPalette(c("green","black","red"))(256)
heatmap.2(x = t(Y.order.2), Rowv =TRUE ,Colv = FALSE, dendrogram = c("row"), scale ="row",ColSideColors =c("green","red","blue","orange")[c.final.order.2], labRow = colnames(Y.order.2), labCol = NA, main = ' \n Training Set \n 47 Gene DPMM signature', col =hmcols, cexRow =0.40, trace ="none", key.title = "Color Key")
legend("topright", legend = c("Best","Good Moderate","Bad moderate","Worst"),fill = c("Green","Blue","Orange","Red"), cex = 0.4)
p5
p1
dev.off()
lograkn
logrank
p5 <- ggsurv(surv.fit, main = " DPMM \n Kaplan Meier Estimators \n pvalue 5e -05", surv.col = c("green","red","blue","orange"), cens.col ="Blue")
pdf('DPMMSignature.pdf')
hmcols<-colorRampPalette(c("green","black","red"))(256)
heatmap.2(x = t(Y.order.2), Rowv =TRUE ,Colv = FALSE, dendrogram = c("row"), scale ="row",ColSideColors =c("green","red","blue","orange")[c.final.order.2], labRow = colnames(Y.order.2), labCol = NA, main = ' \n Training Set \n 47 Gene DPMM signature', col =hmcols, cexRow =0.40, trace ="none", key.title = "Color Key")
legend("topright", legend = c("Best","Good Moderate","Bad moderate","Worst"),fill = c("Green","Blue","Orange","Red"), cex = 0.4)
p5
p1
dev.off()
rank <- matrix(0, nrow = N, ncol =Nps)
for (j in 1:Nps){
Y.scaled <- matrix(0, nrow = N, ncol =D)
for ( v in 1:4){
clust <- which(c.list[[j]] == v)
Y.scaled[clust,1:D] <- scale(Y[clust,1:D], center = TRUE, scale = TRUE)
}
for ( i in 1:N){
rank[i,j] <- dMVN(as.vector(t(Y[i,1:D])), mean = mu.list[[j]][c.final[i],1:D], Q= S.list[[j]][c.final[i],1:D,1:D], log = TRUE) - dMVN(as.vector(t(Y[i,1:D])), mean = mu.list[[j]][1,1:D], Q= S.list[[j]][1,1:D,1:D], log = TRUE) +  dnorm(x = That[i], mean = beta0.list[[j]][c.final[i]] + betahat.list[[j]][c.final[i],1:D] %*% as.vector(t(Y.scaled[i,1:D])), sd = sqrt(sigma2.list[[j]][c.final[i]]), log =TRUE) -  dnorm(x = That[i], mean = beta0.list[[j]][1] + betahat.list[[j]][1,1:D] %*% as.vector(t(Y.scaled[i,1:D])), sd = sqrt(sigma2.list[[j]][1]), log =TRUE)
}
}
avg.rank <- apply(rank,1,mean)
order.zo <- range01(avg.rank)
order.train <- sort(order.zo,index.return = TRUE, decreasing = TRUE)
Y.order <- Y[order.train$ix,]
c.final.order <- c.final[order.train$ix]
######## Reordering Again ################
order.2 <- c(which(c.final.order==1),which(c.final.order==3),which(c.final.order==4),which(c.final.order==2))
Y.order.2 <- Y.order[order.2,]
c.final.order.2 <- c.final.order[order.2]
surv.ob <- Surv(exp(time),censoring)
Classes <- c.final
Classes <- as.factor(Classes)
surv.fit <- survfit(surv.ob ~ Classes)
p5 <- ggsurv(surv.fit, main = " DPMM \n Kaplan Meier Estimators \n pvalue 5e -05", surv.col = c("green","red","blue","orange"), cens.col ="Blue")
############ Generating some Plots ##########################
pc <- prcomp(Y)
pc.pred <- predict(pc,newdata = Y)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], color = as.factor(c.final) )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak \n 18% Variance Explained") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_color_manual(breaks = c("1", "2", "3","4"),values=  c("green","red","blue","orange"))
pdf('DPMMSignature.pdf')
hmcols<-colorRampPalette(c("green","black","red"))(256)
heatmap.2(x = t(Y.order.2), Rowv =TRUE ,Colv = FALSE, dendrogram = c("row"), scale ="row",ColSideColors =c("green","red","blue","orange")[c.final.order.2], labRow = colnames(Y.order.2), labCol = NA, main = ' \n Training Set \n 47 Gene DPMM signature', col =hmcols, cexRow =0.40, trace ="none", key.title = "Color Key")
legend("topright", legend = c("Best","Good Moderate","Bad moderate","Worst"),fill = c("Green","Blue","Orange","Red"), cex = 0.4)
p5
p1
dev.off()
c.list.new.save
c.final.new
c.final.new
That
That.new
time.new
rank <- matrix(0, nrow = N, ncol =Nps)
for (j in 1:Nps){
Y.scaled <- matrix(0, nrow = N, ncol =D)
for ( v in 1:4){
clust <- which(c.new.list.save[[j]] == v)
Y.scaled[clust,1:D] <- scale(Y.new[clust,1:D], center = TRUE, scale = TRUE)
}
for ( i in 1:N){
rank[i,j] <- dMVN(as.vector(t(Y.new[i,1:D])), mean = mu.list[[j]][c.final.new[i],1:D], Q= S.list[[j]][c.final.new[i],1:D,1:D], log = TRUE) - dMVN(as.vector(t(Y.new[i,1:D])), mean = mu.list[[j]][[1,1:D], Q= S.list[[j]][1,1:D,1:D], log = TRUE) +  dnorm(x = time.new[i], mean = beta0.list[[j]][c.final.new[i]] + betahat.list[[j]][c.final.new[i],1:D] %*% as.vector(t(Y.scaled[i,1:D])), sd = sqrt(sigma2.list[[j]][c.final.new[i]]), log =TRUE) -  dnorm(x = time.new[i], mean = beta0.list[[j]][1] + betahat.list[[j]][1,1:D] %*% as.vector(t(Y.scaled[i,1:D])), sd = sqrt(sigma2.list[[j]][1]), log =TRUE)
}
}
c.final.new[i]
dMVN(as.vector(t(Y.new[i,1:D])), mean = mu.list[[j]][[1,1:D], Q= S.list[[j]][1,1:D,1:D], log = TRUE)
S.list[[j]][1,1:D,1:D]
mu.list[[j]][[1,1:D]
]
rank <- matrix(0, nrow = N, ncol =Nps)
for (j in 1:Nps){
Y.scaled <- matrix(0, nrow = N, ncol =D)
for ( v in 1:4){
clust <- which(c.new.list.save[[j]] == v)
Y.scaled[clust,1:D] <- scale(Y.new[clust,1:D], center = TRUE, scale = TRUE)
}
for ( i in 1:N){
rank[i,j] <- dMVN(as.vector(t(Y.new[i,1:D])), mean = mu.list[[j]][c.final.new[i],1:D], Q= S.list[[j]][c.final.new[i],1:D,1:D], log = TRUE) - dMVN(as.vector(t(Y.new[i,1:D])), mean = mu.list[[j]][1,1:D], Q= S.list[[j]][1,1:D,1:D], log = TRUE) +  dnorm(x = time.new[i], mean = beta0.list[[j]][c.final.new[i]] + betahat.list[[j]][c.final.new[i],1:D] %*% as.vector(t(Y.scaled[i,1:D])), sd = sqrt(sigma2.list[[j]][c.final.new[i]]), log =TRUE) -  dnorm(x = time.new[i], mean = beta0.list[[j]][1] + betahat.list[[j]][1,1:D] %*% as.vector(t(Y.scaled[i,1:D])), sd = sqrt(sigma2.list[[j]][1]), log =TRUE)
}
}
avg.rank <- apply(rank,1,mean)
order.zo <- range01(avg.rank)
order.train <- sort(order.zo,index.return = TRUE, decreasing = TRUE)
Y.order.new <- Y.new[order.train$ix,]
c.final.order.new <- c.final.new[order.train$ix]
order.2.new <- c(which(c.final.order.new==1),which(c.final.order.new==3),which(c.final.order.new==4),which(c.final.order.new==2))
Y.order.2.new <- Y.order.new[order.2.new,]
c.final.order.2.new <- c.final.order.new[order.2.new]
surv.ob <- Surv(exp(time.new),censoring.new)
Classes <- c.final.new
Classes <- as.factor(Classes)
surv.fit <- survfit(surv.ob ~ Classes)
p5 <- ggsurv(surv.fit, main = " DPMM \n Kaplan Meier Estimators \n pvalue 3e -02", surv.col = c("green","red","blue","orange"), cens.col ="Blue")
p5
############ Generating some Plots ##########################
pc <- prcomp(Y.new)
pc.pred <- predict(pc,newdata = Y.new)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], color = as.factor(c.final) )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak \n 18% Variance Explained") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_color_manual(breaks = c("1", "2", "3","4"),values=  c("green","red","blue","orange"))
p1
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], color = as.factor(c.final.new) )) + ggtitle(" DPMM Clustering \n 47 Gene Signature Verhaak \n 18% Variance Explained") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_color_manual(breaks = c("1", "2", "3","4"),values=  c("green","red","blue","orange"))
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], color = as.factor(c.final.new) )) + ggtitle(" DPMM Clustering\n Testing Data \n 47 Gene Signature Verhaak \n 18% Variance Explained") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_color_manual(breaks = c("1", "2", "3","4"),values=  c("green","red","blue","orange"))
p1
pc$sd[1] + pc$sd[2]/(sum(pc$sd))
(pc$sd[1] + pc$sd[2])/(sum(pc$sd))
pc <- prcomp(Y.new)
pc.pred <- predict(pc,newdata = Y.new)
p1 <- ggplot(as.data.frame(pc.pred), aes(x=pc.pred[,1], y= pc.pred[,2], color = as.factor(c.final.new) )) + ggtitle(" DPMM Clustering\n Testing Data \n 47 Gene Signature Verhaak \n 18% Variance Explained") + geom_point(shape=19) + labs(y = "PC1", x = "PC2") + scale_color_manual(breaks = c("1", "2", "3","4"),values=  c("green","red","blue","orange"))
(pc$sd[1] + pc$sd[2])/(sum(pc$sd))
hmcols<-colorRampPalette(c("green","black","red"))(256)
heatmap.2(x = t(Y.order.2.new), Rowv =TRUE ,Colv = FALSE, dendrogram = c("row"), scale ="row",ColSideColors =c("green","red","blue","orange")[c.final.order.2.new], labRow = colnames(Y.order.2.new), labCol = NA, main = ' \n Testing Set \n 47 Gene DPMM signature', col =hmcols, cexRow =0.40, trace ="none", key.title = "Color Key")
pdf('TestingDPMMSignature.pdf')
hmcols<-colorRampPalette(c("green","black","red"))(256)
heatmap.2(x = t(Y.order.2.new), Rowv =TRUE ,Colv = FALSE, dendrogram = c("row"), scale ="row",ColSideColors =c("green","red","blue","orange")[c.final.order.2.new], labRow = colnames(Y.order.2.new), labCol = NA, main = ' \n Testing Set \n 47 Gene DPMM signature', col =hmcols, cexRow =0.40, trace ="none", key.title = "Color Key")
legend("topright", legend = c("Best","Good Moderate","Bad moderate","Worst"),fill = c("Green","Blue","Orange","Red"), cex = 0.4)
p5
p1
dev.off()
final.betahat
heatmapdata <- as.data.frame(final.betahat)
dim(heatmapdata)
colnames(heatmapdata) <- signature
signature
colnames(heatmapdata) <- colnames(Y)
heatmap.2(t(as.matrix(heatmapdata)),dendrogram="none", col =cm.colors(180), margins=c(6,10), main = "Posterior prob. \n for Selection \n ", cexCol = 0.85, cexRow = 0.7, Rowv = FALSE)
hmcols<-colorRampPalette(c("white","black"))(256)
heatmap.2(t(as.matrix(heatmapdata)),dendrogram="none", col =hmcols, margins=c(6,10), main = "Posterior prob. \n for Selection \n ", cexCol = 0.85, cexRow = 0.7, Rowv = FALSE)
pdf("Heatmap_FeatureSelection_VerhaakDataSet.pdf")
hmcols<-colorRampPalette(c("white","black"))(256)
heatmap.2(t(as.matrix(heatmapdata)),dendrogram="none", col =hmcols, margins=c(6,10), main = "Posterior prob. \n for Selection \n ", cexCol = 0.85, cexRow = 0.7, Rowv = FALSE, labCol = c("Best","Worst","Moderate_Good","Moderate_Bad"))
dev.off()
heatmap.2(t(as.matrix(heatmapdata)),dendrogram="none", col =hmcols, margins=c(6,10), main = "Posterior prob. \n for Selection \n ", cexCol = 0.85, cexRow = 0.7, Rowv = FALSE, labCol = c("Best","Worst","Moderate_Good","Moderate_Bad"))
pdf("Heatmap_FeatureSelection_VerhaakDataSet.pdf")
hmcols<-colorRampPalette(c("white","black"))(128)
heatmap.2(t(as.matrix(heatmapdata)),dendrogram="none", col =hmcols, margins=c(6,10), main = "Posterior prob. \n for Selection \n ", cexCol = 0.85, cexRow = 0.7, Rowv = FALSE, labCol = c("Best","Worst","Moderate_Good","Moderate_Bad"))
dev.off()
dim(signature.vk)
length(signature.vk)
N
N.new