I would like to calculate a Spearman correlation of gene expression using a big dataset containing around 15000 genes, to generate a figure with the correlation coefficients distribution (as appeared in the attached figure extracted from a paper from Derrien and coworkers, 2012, Genome Res, 22: 1775-1789). The idea is to show the correlation distribution among coding and non-coding genes in a particular type of cells.
I have previously used MAGNET (magnet.case.edu) to calculate the correlation coefficients in smaller datasets, but with this huge amount of genes the produced files are too big for the server, and I am crashing it all the time.
Is there any alternative to do it using another software?. Thanks in advance.
Hi,
If you use R, just import the expression data into R, ideally you should have a matrix of the normalised expression across samples. If the samples are the columns, you will need to transpose the matrix and then run the correlation. for example
if the matrix is called m with samples as columns and genes as rows, then the correlation matrix of all to all genes will be
corm = cor(t(m))
This will likely take about 3-5 minutes in R.
You can then plot the distribution using this command,
plot(density(corm))
This will take some time, as it will try to flatten your matrix into a big vector of length 15000x15000 and then will plot the density, these operations are mainly memory intensive operations as you need to store all this in the memory and calculate the density. If you have a big memory machine, these commands will run much faster there.
Hope this helps,
Best Regards,
Rafik
Hi Rafik
Many thanks for your suggestions.... With your indications and "a little help from my friends..:"... It's working. Thanks !!!
Cheers
i can add my two bits here as well :)
usually if the sample size is too big, i.e. the 15k X 15k matrix you are expecting. and you want to draw plots or calculate a statistic - then you can take a sample from the data say about 10,000, using the sample function in R. this will reduce your computational costs significantly and won't affect your end result i.e. plots of density.
furthermore, if the expression data is from RNA-Seq, then i would also advise taking a log of the data first i.e. log(data+1) - the plus one is there for avoiding taking logs of zeros. this is due to the nature of RNA-Seq data as when the mean expression value for a particular gene increases, you also have an increase in variance. taking the log will perform 'variance stabilization' of the data.
Thanks for the question and discussion. I would like to try something similar. I have expression values for mRNA and miRNA from sequencing of 9 tissue samples. Would using baseMeans from DESeq2 outfile make sense?
- Badges
- Science topic
- Similar topics
- Computer Science and Engineering
- Bioinformatics
- Bioinformatic Tools