Hi , I am trying to perform a big number of case control status permutations, this is my code:
Y<-matrix(0,1285,100000)
Y<-apply(Y,c(1,2), function(x) sample(c(0,1),1))
NAMES <- names(df[,51:7298])
res<-as.data.frame(matrix(nrow=100000,ncol=length(NAMES),dimnames=NULL))
for (s in 1:100000) {
for(i in 51:7298)
{
fit<-glm(Y[,s]~df$cov1+df[,i]+df$cov2+ df$cov3,data=df,family=binomial())
sfit <-summary(fit)
j=i-50
res[s,j]<-as.numeric(sfit$coefficients[3,4])
colnames(res)<-NAMES
}}
It's too slow under R, will it be better to use C++ or Python ? I have a beginner level in both Python and C++language, any script suggestions are welcome !
Thank you
1 answer
I would rather invest into optimising code in a language you know well rather than trying out several ones at which you are a beginner. Check out parallelization in R, be it using foreach or functions like mclapply from the parallel package. Also, critically review your code towards redundant operations that might slow down execution. What exactly is this code doing and why? What is the input data format and what is the final output?
Your code is actually fine, it is simply a lot of data. I suggest you replace the outer loop by mclapply or any of the available parallelization commands, get a HPC node with a lot of memory and many CPUs and then wait till everything is finished. I do not know about memory usage but if you use, say a 72-core node if available, this should speed up things dramatically with mclapply.
Log in to answer this question.
-
Bumphunter function in minfi
written by eebloom •
-
Cochran–Armitage test results
written by Eliza •
-
CNV annotation using GAIA
written by melissachua90 •
-
Tutorial: Gene Set Enrichment Analysis
written by Hamid Ghaedi
-
Tutorial: Methylation Analysis Tutorial in R_part2
written by Hamid Ghaedi
-
Compute the correlations between submatrices
written by pablo
-
edgeR DE; choosing right coefficient from design matrix with many covariates
written by Cece •
-
Tutorial: Network plot from expression data in R using igraph
written by Kevin Blighe
-
WGCNA on gene abundance data
written by cara78 •
-
Tutorial: Survival analysis of TCGA patients integrating gene expression (RNASeq) data
written by TriS
If time really matters, C++ is the best.
Yours is an embarrassingly parallel problem. Learning to recognize such cases and how to handle them will go a long way in helping you deal with heavy computations. As already mentioned, check the various R packages available for this. Also this is where you'd make use of a compute cluster if you've access to one. Recoding in a fast language (C, C++...but not python) should be the absolute last resort (e.g. when all you can use is a laptop) or done as an exercise.
Thank you very much Jean-Karim, The problem that I am failing even on a compute cluster. The loop I wrote take a lot of memory exceeding the cluster performance (Or may be I am not using the best way ... ) That's way I am lokking how to mek better this or to write it in anothjer language ... C++ still not easy for me as a beginner
Start by ensuring that the outer loop sends jobs to different nodes of your cluster. Rewriting without parallelization may not solve your memory issue if your cluster nodes don't have the required RAM. The way your code is written you keep data for each run in memory so your memory use keeps on growing. Your res matrix is about 5.5 GB, your Y matrix is 1 GB so this piece of code needs at least 7 GB without including the data in df and space required to compute the regression model. Also stick with matrices if you don't need the flexibility/features of data frames. Accessing elements of a data frame can be several orders of magnitude slower than with a matrix as you can see with this:
On my machine, this shows the data frame to be ~50 times slower than the matrix.