How to get Cox p-value for 20,000 genes?

If you run the following code, you will have a data frame real.dat which has 1063 samples for 20531 genes. There are 2 extra columns named time and event where time is the survival time and event is death in case of 1 and 0 in case of censored.

lung.dat <- read.table("genomicMatrix_lung")

# For clinical data, get only rows which do not have NA in column "X_EVENT"
lung.no.na.dat <- lung.clin.dat[!is.na(lung.clin.dat$X_EVENT), ] # Getting the transpose of main lung cancer data ge <- t(lung.dat) # Getting a vector of all the id's in the clinical data frame without any 'NA' values keep <- lung.no.na.dat$sampleID

# getting only the samples(persons) for which we have a value rather than 'NA' values
real.dat <- ge[ge[, 1] %in% keep, ]

# adding the 2 columns from clinical data to gene expression data
keep_again <- real.dat[, 1]
temp_df <- lung.no.na.dat[lung.no.na.dat$sampleID %in% keep_again, ] # naming the columns into our gene expression data col_names <- ge[1, ] colnames(real.dat) <- col_names dd <- temp_df[, c('X_TIME_TO_EVENT', 'X_EVENT')] real.dat <- cbind(real.dat, dd) # renaming the 2 new added columns colnames(real.dat)[colnames(real.dat) == 'X_TIME_TO_EVENT'] <- 'time' colnames(real.dat)[colnames(real.dat) == 'X_EVENT'] <- 'event'  I want to get the univariate Cox regression p-value for each gene in the above data frame. Now, when I try to run the coxph function from the survival package even for one gene, it shows the following error - > coxph(Surv(time, event) ~ HIF3A, real.dat) Error in fitter(X, Y, strats, offset, init, control, weights = weights, : NA/NaN/Inf in foreign function call (arg 6) In addition: Warning message: In fitter(X, Y, strats, offset, init, control, weights = weights, : Ran out of iterations and did not converge What am I doing wrong here? You can download the data from here. 1 Answer I think the problem is that, for some reason, the genes in your data are not coded as numeric variables. For the HIF3A, for example: > str(real.dat$HIF3A)
Factor w/ 1059 levels "0.0000","0.3468",..: 309 327 294 108 611 458 347 346 576 417 ...
- attr(*, "names")= chr [1:1063] "V2" "V3" "V4" "V5" ...


It should be numeric for the Cox model. The error you get is because coxph runs into trouble with fitting a model with a factor with 1059 levels. So yeah, check your data before fitting models.

NB: when converting factors to numeric, you should use as.numeric(as.character( )), as you can see from this:

> real.dat$HIF3A[1:5] V2 V3 V4 V5 V6 4.6261 4.7796 4.5313 2.2999 6.6477 1059 Levels: 0.0000 0.3468 0.4397 0.5484 0.5720 0.6396 0.9456 0.9635 0.9674 1.0976 1.1084 ... 9.9993 > as.numeric(real.dat$HIF3A)[1:5]
[1] 309 327 294 108 611
> as.numeric(as.character(real.dat\$HIF3A))[1:5]
[1] 4.6261 4.7796 4.5313 2.2999 6.6477

• Hi Theodor, many many thanks for taking the time to look at the data. That was indeed the problem here. Also, many thanks for suggesting the solution :). – nafizh Apr 20 '16 at 6:10