3
\$\begingroup\$

I have written a code for a project which aims at finding the Cox regression coefficients in a mixed model for microarray data. The study was carried out on the Affymetrix Hgu133a platform. In the following code, I have first attempted to collapsing the probes to genes with the WGCNA package, and then carry out cox mixed effects models in a loop using the batch date as a random effects variable.I will be grateful if anyone could just check if the code is right, because it takes a very long time in giving results. I want to know not just about speeding up the code , but also if the code is correct.

Just for clarity:

  • probes = affymetrix probes
  • genes = genes which are mapped to the affymetrix probes
  • tdmfs = time variable
  • edmfs = event variable
  • chrdate = date of scanning the affymetrix data (proxy for batch effects)
  • totcomb = data containing the expression data with the expression from probes in columns, also contains the time, event and chrdate variable.

This is a different arrangement from the usual microarray data, where the expression data are usually arranged in rows.

This is a part of a much larger project where the same code is run in different random models. But the general format is the same, so if this gives a correct result, the Cox regression results for the others will also be correct, and any way this code snippet can be made to run quicker, will make the entire project quicker.

#### The expression data and the other data are separated,
#### since only the expression data will be used for collapsing the probes to genes
zz=length(colnames(totcomb))
data=totcomb[,c(1:22283)] #### expression data since affymetrixhgu133a has 22283 probes and they are in the first 22283 columns of the data frame
demo=totcomb[,c(22284:zz)]


##### For the collapsing the probes are supposed to be in rows and the samples in columns, so the data has to be transposed
data=t(data)

#### Only that data representing the probes having a mapping to a Entrez Gene ID are kept. 
data=data[c(row.names(data) %in% probes),]
#### The probes are then collapsed to genes. The gene names are represented in the row names
datacol=(collapseRows(data,genes,probes))[[1]]
#### The number of genes finally the data is collapsed to 
len=length(row.names(datacol))
#### The collapsed  expression data is to be transposed again so that the genes are in the columns which is needed for the cox regression that is to follow
datacol=t(datacol)
#### The expression data for the genes are then merged again with the data containing the estimated molecular subtypes, AURKA score, scan date and follow up and event times
datacol=merge(datacol, demo, by="row.names")
row.names(datacol)=datacol$Row.names
datacol$Row.names=NULL 

####Make an empty matrix to store results
resmat=matrix(nrow=len, ncol=4)

####Model with the scan date serving as a proxy for batch and controlled as a random effects variable in a mixed effects Cox model. 
#### The Cox regression coefficients along with the z score and the AIC of the model are calculated for each gene through a loop and the values entered into the empty matrix (column 1 contains the coefficients, column 2 the zscore for the genes and column 3 the AIC )

####Model (i): for just the batch effects modeled as a random variable 
for (i in 1:len) {
         coxcomb=coxme(Surv(tdmfs,edmfs)~datacol[,c(i)]+(1|chrdate), data=datacol)
                 resmat[c(i),c(1)]=fixef(coxcomb)
         resmat[c(i),c(2)]=resmat[c(i),c(1)]/sqrt(diag(vcov(coxcomb)))
         resmat[c(i),c(3)]=-2*(coxcomb$loglik[2])+2*coxcomb$df[1]
}
#### The coefficients are then mapped to their Entrez IDs
resmat[,c(4)]=names(datacol)[1:len]
\$\endgroup\$
  • 1
    \$\begingroup\$ 1. Variables like "coxcombunivlistrandombatch" make your code rather difficult to read. 2. Did you try out your code? What happens? I'd recommend to choose a (much) smaller sub sample to validate your model first, if time is an issue. \$\endgroup\$ – Vincent Oct 7 '14 at 5:19
  • \$\begingroup\$ I tried out the code, and it gives a result, but it takes very long -about 30 minutes for the above code. Also sorry about the long names- as I mentioned in the main question, this is a part of a much larger project, where 36 random models are compared and takes more than 20 hours to run. The names are such that I can differentiate between the models. \$\endgroup\$ – Nilotpal Oct 7 '14 at 8:27
  • \$\begingroup\$ Just checked cran.r-project.org/web/packages/coxme/index.html. R provides great functionality in terms of ease of use and straight forward implementation of model evaluation. That comes at the price of speed. Looking at the documentation of the website it seems to me that this a rather computationally intensive procedure, so I really would not be surprise about the time it takes to compute a large data set multiple times. \$\endgroup\$ – Vincent Oct 7 '14 at 9:05
  • \$\begingroup\$ Thanks to all who have made my initial post more readable \$\endgroup\$ – Nilotpal Oct 7 '14 at 9:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.