Hello, I've been using the survival package to determine whether a gene's expression can be used in prediction of survival probabilities a while now using the CDGSR package to pull TCGA data. However, now that I've been thinking about using multiple genes' expression to stratify patients, I'm a bit stuck. I've written some code in R to do just that based on a paper I've read but can't remember the name anymore where the basic idea is to score each patient 1 or -1 based on whether they express a gene more than the mean expression of that gene in the population. Afterwards, the scores are added and each patient end up with a PIR score (I remember the original paper I've read used the PIR score but can't remember what it stands for now, I've just kept it PIR in case I run into that paper again), I then stratify the patients based on the median of PIR scores and perform the survival analysis. Would anyone recommend a better way of doing things? I'm pasting the code below - which is probably very crude as I'm not too experienced with R yet. Any suggestion/critique/idea would be greatly appreciated. library(cgdsr, survival, reshape, plyr) expression <- data.frame(rep(1, 563)) #Create a dataframe to put the date in for (i in 1:length(chr.list)) { #chr.list is a list of genes x <- getProfileData(mycgds, chr.list[i], "ov_tcga_mrna_median", "ov_tcga_all") if(length(table(x)) != 0 ) { #Check to make sure the gene is represented and expression values are obtained expression <- cbind(expression, x) } } #get the expression data from TCGA database d.frame <- na.omit(expression[, 2:length(expression)]) #Remove missing values d.frame$patient <- rownames(d.frame) # add a column with patient IDs d.frame.melt <- melt(d.frame, id = "patient") # melt the data d.frame.melt <- ddply(d.frame.melt, .(variable), transform, mean = mean(value)) # create a new column, with mean value for each gene in the group of patients d.frame.melt <- ddply(d.frame.melt, .(patient), transform, PIR = ifelse(value > mean, 1, -1)) #For each gene, if the expression is higher td.framen the mean, then the PIR score is 1, else it is -1 d.frame.PIR <- ddply(d.frame.melt, .(patient), colwise(sum, "PIR")) # Sum the PIR values d.frame.PIR <- ddply(d.frame.PIR, .(patient), transform, survival = ifelse(PIR > median(d.frame.PIR$PIR), 1, 0)) #stratify patients based on the PIR value clinical <- getClinicalData(mycgds, "ov_tcga_all") #Get clinical data clinical$patient <- rownames(clinical) d.frame.PIR <- join(d.frame.PIR, clinical, type = "inner") #merge clinical information with PIR stratified patient data frame d.frame.PIR$status <- d.frame.PIR$overall_survival_status=="DECEASED" d.frame.surv <- Surv(d.frame.PIR$overall_survival_months, d.frame.PIR$status) # run the survival analysis and draw a graph d.frame.diff <- survdiff(d.frame.surv ~ d.frame.PIR$survival) d.frame.fit <- survfit(d.frame.surv ~ d.frame.PIR$survival) p <- 1-pchisq(d.frame.diff$chisq, length(d.frame.diff$n) -1) plot(d.frame.fit, lty=1:2, mark.time=F) legend(10, .2, c("HOXB1-13 high expression", "HOXB1-13 low expression"), lty = 1:2) title(main=paste("Survival Probablity of TCGA Ovarian Patients based on HOXB Cluster (HOXB1-13) expression (p = ", round(p, digits = 4), ")", sep=""), xlab="Months", ylab="Probability") > sessionInfo() R version 2.14.2 (2012-02-29) Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) locale: [1] C/en_US.UTF-8/C/C/C/C attached base packages: [1] splines graphics grDevices utils datasets grid stats methods base other attached packages: [1] survival_2.36-12 reshape2_1.2.1 reshape_0.8.4 plyr_1.7.1 cgdsr_1.1.19 R.oo_1.8.3 R.methodsS3_1.2.2 RSQLite_0.11.1 [9] DBI_0.2-5 knitr_0.3 gplots_2.10.1 KernSmooth_2.23-7 caTools_1.12 bitops_1.0-4.1 gdata_2.8.2 gtools_2.6.2 [17] RColorBrewer_1.0-5 ggplot2_0.9.0 loaded via a namespace (and not attached): [1] MASS_7.3-17 Rcpp_0.9.10 codetools_0.2-8 colorspace_1.1-1 dichromat_1.2-4 digest_0.5.1 evaluate_0.4.1 formatR_0.3-4 highlight_0.3.1 [10] memoise_0.1 munsell_0.3 parser_0.0-14 proto_0.3-9.2 scales_0.2.0 stringr_0.6 tools_2.14.2 -- *Ahmet Z.* [[alternative HTML version deleted]]