x=read.table("",head=T,row.names=1,sep="\t")

(if column 1 is your gene column)

cor(x,method="pearson")

will give you the correlation matrix

If you want the p-values rather than using apply you can try the "psych" package, it requires "mnormt".

library("psych")

corr.test(x)$p (the spelling corr.test is the correct one for the function provided by psych) will give you a matrix with the correlation p values

What you can next do is to cluster them by correlation coefficient, either in R, or using external tools like cluster 3.0, gives nice pairwise (symmetrical) dendograms

Dear Mohamed,

I think the WGCNA (https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/) is suitable for your aim.

Calculating all the pairwise correlation for 60k genes would be computationally intractable usually. I would remove uninformative genes first, such as 1) the ones that are not expressed or very lowly expressed in most samples or 2) the ones whose expression level is almost identical or displays small variation throughout samples. This initial filtering helps you not only to reduce the computational overhead but also to get more clearly distinct group of genes of interest for downstream analysis. Then you can calculate the pairwise correlation with the default R function, 'cor'.

'cor' only gives you the correlation coefficients, not P-values. You can use 'cor.test' to get both.

To make correlation matrix you can use 'cor.test' on one-transcript-at-a-time or on everything at once. The slow and intuitive way with a "for-loop". Or faster, but less intuitive way using "apply".

The result from 'cor.test' is a list. You can get the correlation coefficient and p-values like this:

cor.test(x, y, method="pearson")$p.value

cor.test(x, y, method="pearson")$est

x=read.table("",head=T,row.names=1,sep="\t")

(if column 1 is your gene column)

cor(x,method="pearson")

will give you the correlation matrix

If you want the p-values rather than using apply you can try the "psych" package, it requires "mnormt".

library("psych")

corr.test(x)$p (the spelling corr.test is the correct one for the function provided by psych) will give you a matrix with the correlation p values

What you can next do is to cluster them by correlation coefficient, either in R, or using external tools like cluster 3.0, gives nice pairwise (symmetrical) dendograms

I think you can use this, but in my opinion you should absolutely calibrate the colors and your threshold values for "pheatmap" because otherwise you'll have a hard time comparing several heatmaps, or if your values are very close, you may have 2 or 3 different color intensities for a variation of 0.1 which is low (a value of 0.8 can be blue and a value of 0.9 can be red) :

### Bioconductor Packages

source("http://bioconductor.org/biocLite.R")

biocLite()

biocValid()

biocLite(c("pheatmap"))

### Correlation

#DATA = data.frame with gene in row and sample in colum

cormM = cor(DATA, method = "pearson")

### Heatmap

## Set Value for legend

heatmap.MIN.Val = max(0, (trunc(min(cormM) * 10) - 1) / 10)

heatmap.MAX.Val = 1

heatmap.RANGE.Val = (heatmap.MAX.Val - heatmap.MIN.Val) / 100

## Set color

# If you want to make it automaticaly, leave "annotation" and "annotation_colors" when you call pheatmap function

grouplabelcol = rainbow(ncol(DATA), s=1, alpha = 1)

# OR WITH YOUR OWN COLORS

grouplabelcol = c("red", "blue", "green", "violet")

names(grouplabelcol)=c(colnames(DATA))

anno_colors

RNA-seq data usually follow a poisson distribution which prevents to use Pearson. Make sure that your data follow a normal distribution and that the variance are equal. You can use shapiro-wilk test, anderson-darling test and the "fat pencil" test which is less stringent. Otherwise I would go for a Spearman correlation.

If you want to make a correlation table, have a look at the function "pairs" in R

Edit: Also try to compare log with RPKM (FPKM if paired) and rlog to transform your data.

Just adding on the previous answers, loading the WGCNA package is a good idea, because it has a much faster "cor" function that the basic one.

If you go with the cor.test() option I would also adjust my p-values for multiple comparisons:

adjusted = p.adjust(p_values_from_cor.test, "BH") # BH correction

I agree with Lim's comment that you should filter your data first. I have previously used the WGCNA package and filtered the genes based on variance (i.e. kept the most variable genes - upper 75 percentile):

countMatrix = read.csv("count_matrix.txt", header=TRUE, sep="\t", check.names=FALSE, row.names=1)

#Filter out lowly expressed genes, in this case keep genes expressed more than 5 in at least 3 samples. But could be too strict

use = rowSums(countMatrix > 5) >=3

countMatrixFiltered = countMatrix[ use, ]

#log2 transform (x+1 to avoid zeroes)

log2countMatrixFiltered = log2(countMatrixFiltered+1)

#Extract the genes with the highest variance (upper 75 percentile).

library(genefilter)

rv = rowVars(log2countMatrixFiltered)

q75 = quantile(rowVars(log2countMatrixFiltered), .75)

q75log2countMatrixFiltered = log2countMatrixFiltered[rv > q75, ]

I guess your data must be really big. Maybe you want the pairwise correlation of ~20k genes (or even more)? If that so, I would use R which handles matrices really well (and easy to use too).

Software for that task would be challenging. So far I've used BioVinci for calculating and visualizing pairwise correlation. It is not as good as coding, but for visualization, it did well for me (you can download it at https://vinci.bioturing.com)

Dear authors ,how about if I want to pick out the gene list that is most correlated with a specific gene , say PD-1,by the Spearman correlation analysis?(

Pearson correlation) Great thanks!

This could be useful https://reneshbedre.github.io/blog/corr.html