How to Run Correlation Analysis on Your Own RNA-Seq Data

Oak
rna-seq correlation coexpression multi-omics bioinformatics
A shallow reef at dusk filled with rows of glowing anemone stalks, some pulsing green in unison and connected by faint glowing threads.

By the end of this guide you will have a gene-by-sample expression matrix on a variance-stabilized scale, a sample-sample correlation matrix you can use as a quality check, a filtered gene-gene co-expression matrix with the dominant technical structure removed, and a set of correlations between your transcripts and external measurements such as plasma proteins, blood chemistry, or continuous glucose readings. You need FASTQ files (or an already-quantified count matrix), a reference transcriptome, R 4.3 or later with DESeq2, limma, WGCNA, and ppcor, and roughly 16 GB of RAM. You also need a sample sheet that records, for every library, the collection date and time, the library preparation protocol, the sequencing run, the RNA integrity number (RIN), and whatever phenotypic covariates you have. That sample sheet matters more than any of the code below, because almost every spurious correlation in RNA-seq traces back to a variable that was not written down.

1. Quantify consistently, and never mix protocols inside one correlation matrix

Correlation is a comparison, so every sample in the matrix must have been measured the same way. Quantify all libraries against the same reference and the same tool version. We use Salmon for personal-scale work because it is fast enough to re-run the whole archive when you change references, and its bias models matter for correlation specifically: GC and positional bias vary between runs and will show up as run-level structure in your gene-gene correlations.

salmon index -t gencode.v44.transcripts.fa.gz -i salmon_idx_v44 --gencode -k 31

for s in $(cut -f1 samples.tsv | tail -n +2); do
  salmon quant -i salmon_idx_v44 -l A \
    -1 fastq/${s}_R1.fastq.gz -2 fastq/${s}_R2.fastq.gz \
    --validateMappings --gcBias --seqBias \
    --numGibbsSamples 20 --threads 16 \
    -o quant/${s}
done

The --numGibbsSamples 20 flag gives you posterior samples of each transcript’s abundance, which is the cheapest way to find out whether a gene’s estimate is stable enough to correlate at all. Genes with high inferential variance (paralog families, short transcripts, anything in an immunoglobulin locus) produce correlations that are an artifact of read assignment ambiguity.

Do not pool libraries built with different chemistries. Ribosomal-RNA-depleted and polyA-selected libraries from the same RNA sample measure overlapping but distinct transcript populations, with the depleted libraries capturing non-polyadenylated and intronic signal that polyA selection removes, so protocol becomes the largest source of variance in any joint matrix.1 If your longitudinal archive contains both, analyze them as separate matrices and compare the results, rather than concatenating them and hoping a batch correction absorbs the difference.

Import to gene level with tximport, scaling by effective transcript length so that gene-level counts remain comparable when isoform usage shifts:

library(tximport); library(readr)
tx2gene <- read_tsv("tx2gene.gencode.v44.tsv")
files <- file.path("quant", samples$id, "quant.sf"); names(files) <- samples$id
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

2. Pick the transform before you pick the coefficient

The single most common error in RNA-seq correlation analysis is computing Pearson correlation on raw TPM or FPKM values. Count data is heteroscedastic: the variance of a gene’s estimate grows with its mean, so a handful of high-expressing samples dominate the covariance, and a gene expressed at 5 counts carries the same nominal weight as one at 5,000 while having ten times the relative noise. The fix is a variance-stabilizing transform.

library(DESeq2)
dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ 1)
keep <- rowSums(counts(dds) >= 10) >= ceiling(0.5 * ncol(dds))
dds <- dds[keep, ]
vsd <- vst(dds, blind = TRUE)
E <- assay(vsd)   # genes x samples, roughly homoscedastic, log2-like scale

Use blind = TRUE when the matrix is going into an exploratory correlation analysis, so the dispersion estimates are not informed by a design you are about to test. If you have fewer than about 15 samples, rlog() is the better-behaved choice; above that, vst() is much faster and nearly equivalent.

Two further notes on scale. First, TPM is compositional: it sums to a fixed total, so an increase in one large gene mathematically forces decreases elsewhere, producing negative correlations that reflect arithmetic rather than biology. Median-of-ratios normalization inside DESeq2 mitigates this but does not eliminate it, which is why a few strongly induced genes (hemoglobin in a hemolyzed draw, interferon genes during a viral infection) can shift the entire correlation matrix. Second, if you only have a TPM table and no counts, log2(TPM + 1) followed by quantile normalization is a serviceable fallback, and Spearman correlation is safer than Pearson on that scale.

3. Use sample-sample correlation as your first quality check

Before asking anything about genes, ask whether your samples make sense. Compute correlations between samples on the 2,000 most variable genes, which is where structure lives, and inspect the matrix as a heatmap ordered by collection date.

vg <- head(order(rowVars(E), decreasing = TRUE), 2000)
S <- cor(E[vg, ], method = "pearson")
pheatmap::pheatmap(S, annotation_col = samples[, c("protocol","run","RIN","date")])

For repeated sampling of one person from the same tissue, adjacent-in-time pairs typically sit above 0.95 on this scale, and anything below about 0.90 deserves investigation before it enters downstream analysis. When you see blocks, read them against the annotation columns. Blocks aligned with sequencing run are a batch effect. Blocks aligned with RIN are degradation, which selectively depletes long transcripts and 5’ ends and will create a spurious correlation module made entirely of transcript length. Blocks aligned with nothing in your sample sheet mean a variable you failed to record, and in whole blood the usual suspect is cell composition: the first principal component of a blood RNA-seq series is frequently the neutrophil-to-lymphocyte ratio, so a paired complete blood count with differential is worth having for every draw.

4. Filter aggressively, because low counts manufacture correlation

The filter in step 2 (at least 10 counts in at least half the samples) is not cosmetic. At low counts the Poisson component of the noise dominates, and two genes that share nothing biologically will still co-vary because both track library depth and both spend part of their range near zero. The result is a long tail of high-coefficient, low-information gene pairs, and in a co-expression network those pairs form dense modules of poorly expressed genes that resist interpretation.

We also drop genes with unstable quantification, using the Gibbs posteriors from step 1:

# inferential relative variance per transcript, aggregated to gene level
library(fishpond)
se <- tximeta::tximeta(coldata)             # imports inferential replicates
se <- fishpond::computeInfRV(se)
drop <- rownames(se)[mcols(se)$meanInfRV > 0.5]
E <- E[!rownames(E) %in% drop, ]

A practical answer to the read-depth question, since it determines what you can correlate: gene-level differential expression and correlation on moderately expressed genes is comfortable at 20 to 30 million paired reads per library, and correlations among the top few thousand expressed genes are already stable around 10 million. Transcription factors, cytokines, and other low-abundance transcripts need considerably more depth, and no amount of statistical care recovers a correlation for a gene sitting at 3 counts. Isoform-level correlation is a different regime again and benefits from 50 million reads or more with long fragments.

5. Choose the coefficient deliberately, and know what ANOVA does instead

On variance-stabilized data, Pearson correlation is the default we use: it is the maximum-likelihood estimator under approximate joint normality, it is fast enough for a 15,000-gene square matrix, and it feeds directly into the network and factor methods below. Spearman’s rank correlation is the right choice when you have not variance-stabilized, when a single timepoint is an outlier you do not want to delete, or when you expect monotone but strongly nonlinear relationships, such as a saturating response. The cost is power: on clean data Spearman typically needs roughly ten percent more samples for the same detection ability, and ties among low counts distort ranks. Biweight midcorrelation, implemented as WGCNA::bicor(x, y, maxPOutliers = 0.1), is a good middle ground and what we use for network construction.

To answer a question that comes up constantly: ANOVA is not a correlation test. Analysis of variance partitions variance across discrete groups and asks whether group means differ, so it is the tool for “does expression differ between my fasted and fed draws.” Correlation asks whether two continuous variables co-vary monotonically or linearly. The two are connected through the linear model, where a one-way ANOVA on two groups is equivalent to a t-test and to the test on a single binary regression coefficient, and the F-test on a simple regression of y on continuous x is algebraically the same test as Pearson’s r. So if your predictor is continuous, use correlation or regression. If it is categorical with three or more levels, use ANOVA or its count-based analogue, the likelihood ratio test in DESeq2. If you have repeated measurements on one person and want both, a linear mixed model with a subject or batch random effect is the correct frame.

6. Remove the dominant confounders before computing gene-gene correlations

Any variable that affects thousands of genes at once, such as sequencing batch, RIN, or cell composition, induces correlation between every pair of genes it touches. This is the most consequential decision in the whole pipeline, and we would rather over-correct a known technical factor than interpret a module that turns out to be a batch.

library(limma)
design <- model.matrix(~ 1, data = samples)              # or ~ condition to preserve
covars <- model.matrix(~ RIN + neut_frac + lymph_frac, data = samples)[, -1]
E_adj <- removeBatchEffect(E, batch = samples$run, covariates = covars, design = design)

For structure you cannot name, estimate surrogate variables rather than guess:

library(sva)
mod  <- model.matrix(~ condition, data = samples)
mod0 <- model.matrix(~ 1, data = samples)
n_sv <- num.sv(E, mod, method = "be")
svobj <- sva(E, mod, mod0, n.sv = n_sv)

Include the surrogate variables as covariates in removeBatchEffect, but check first that they do not correlate with the biological variable you care about. If SV1 tracks your condition of interest, you will delete the signal along with the noise. An alternative that avoids the deletion problem entirely is partial correlation: keep the full matrix and condition on the covariates pairwise.

library(ppcor)
pcor.test(E["ENSG00000115263", ], glucose_auc, cbind(samples$RIN, samples$neut_frac),
          method = "pearson")

7. Build co-expression modules instead of reading 100 million pairs

A filtered matrix of 15,000 genes has about 112 million pairs, which is too many to interpret and too many to correct for individually. Reduce to modules. WGCNA clusters genes by topological overlap, which is more stable than raw correlation because it asks whether two genes share neighbors rather than only whether they track each other.

library(WGCNA); enableWGCNAThreads(16)
datExpr <- t(E_adj)
sft <- pickSoftThreshold(datExpr, powerVector = 1:20, networkType = "signed")
net <- blockwiseModules(datExpr, power = sft$powerEstimate,
                        networkType = "signed", corType = "bicor",
                        maxPOutliers = 0.10, minModuleSize = 30,
                        mergeCutHeight = 0.25, deepSplit = 2,
                        numericLabels = TRUE, maxBlockSize = 20000)
MEs <- moduleEigengenes(datExpr, net$colors)$eigengenes

Use networkType = "signed" so that anti-correlated genes do not land in the same module, which matters when you later want to interpret a module eigengene’s direction. The eigengenes, one value per module per sample, are the objects you correlate against phenotype in the next step, and reducing 15,000 tests to 30 makes the multiple-testing problem tractable.

Two cautions. WGCNA on fewer than about 15 samples produces unstable modules, and pairwise correlation of any kind captures only two-way dependence. Methods built on higher-order dependence structures can recover relationships that a correlation matrix misses, particularly conditional and multi-gene interactions.2 If your modules look uninformative and you have depth, that is the direction to explore rather than tuning mergeCutHeight indefinitely.

8. Correlate expression against your other measurements

This is where a personal molecular profile earns its keep, because transcript abundance in isolation is hard to interpret while a transcript that tracks a plasma protein or a glucose excursion is a lead. The simple version correlates module eigengenes against each phenotype and controls the false discovery rate:

pheno <- samples[, c("hba1c","crp","apob","glucose_cv","tg_hdl")]
R <- cor(MEs, pheno, use = "pairwise.complete.obs", method = "spearman")
P <- WGCNA::corPvalueStudent(R, nSamples = nrow(MEs))
Q <- matrix(p.adjust(P, "BH"), nrow = nrow(P), dimnames = dimnames(P))

Use Spearman here, because clinical analytes are frequently skewed and a single high CRP value from an incidental infection will otherwise drive the result. For continuous glucose data, decide what you are correlating against before you compute anything: mean glucose, coefficient of variation, time in range, and postprandial area under the curve are different quantities with different noise properties, and the CGM window should be aligned to the blood draw rather than to the calendar day. If you have dense timepoints, lagged cross-correlation is informative, since a transcriptional response to a metabolic state need not be simultaneous.

ccf(as.numeric(MEs$ME3), as.numeric(glucose_daily_cv), lag.max = 7, plot = TRUE)

Interpret lagged results cautiously. Both series are autocorrelated in time, and standard confidence bands assume independence, so they are far too narrow. Test significance by circular block permutation of one series instead of trusting the default bands.

When you have several omics layers measured on the same samples, pairwise correlation across layers scales badly and wastes the shared structure. Multivariate integration methods designed for multi-block data, such as multiple co-inertia analysis, find axes of common variation across transcriptomic, proteomic, and metabolomic blocks and are better suited to this problem than a large cross-correlation matrix.3 Studies that pair transcriptomic and plasma proteomic profiling routinely find that integrated axes carry signal absent from any single layer, which is the practical argument for measuring more than one.4 If you would rather run a pipeline than assemble the pieces, workflow implementations exist that wire together expression, correlation, and enrichment steps for multi-omics input.5

Everything in this section is measurement and hypothesis generation. A correlation between a module eigengene and HbA1c in your own data describes covariation in one person’s samples, and it does not establish causation, diagnose anything, or indicate what to do next. Take any finding that touches on health to a clinician who can put it in the context of your history and examination.

Common problems

Correlations near 1.0 across most gene pairs mean an uncorrected global factor is still present, usually library depth or degradation. Recompute the sample-sample correlation matrix and check the annotation alignment before adjusting anything else.

A module made entirely of long or short transcripts is a degradation or 3’ bias artifact. Regress transcript length against module membership to confirm.

Correlations that vanish after you add one more sample were never real. With 8 timepoints, the 95 percent confidence interval on an observed r of 0.7 runs roughly from 0.06 to 0.94, so treat any single-sample-sensitive result as noise and prioritize adding timepoints over adding statistical sophistication.

Strong negative correlations between highly expressed genes are often compositional. Verify on counts with median-of-ratios normalization rather than on TPM.

Repeated samples from one person are not independent observations. A naive Pearson p-value overstates significance when adjacent timepoints are similar, so use block bootstrap or permutation that preserves the temporal ordering.

Finally, a correlation is not a mechanism. Even well-controlled co-expression networks mix direct regulation, shared upstream drivers, and cell-composition shifts, which is why biomarker work treats correlation as a screening step feeding orthogonal validation.6

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Li Chen, Ruirui Yang, Tony Kwan, et al. Paired rRNA-depleted and polyA-selected RNA sequencing data and supporting multi-omics data from human T cells. Scientific Data, 2020. https://doi.org/10.1038/s41597-020-00719-4 ↩

  2. Hui Tang, Tao Zeng, Luonan Chen. High-Order Correlation Integration for Single-Cell or Bulk RNA-seq Data Analysis. Frontiers in Genetics, 2019. https://doi.org/10.3389/fgene.2019.00371 ↩

  3. Chen Meng, Bernhard Kuster, Aedín C Culhane, et al. A multivariate approach to the integration of multi-omics datasets. BMC Bioinformatics, 2014. https://doi.org/10.1186/1471-2105-15-162 ↩

  4. Xuan Lai, Shengyuan Xu, Fuchun Zhang, et al. Plasma multi-omics profiling reveals unique precursors of post-COVID cognitive decline during clinical recovery in older patients. Journal of Pharmaceutical Analysis, 2026. https://doi.org/10.1016/j.jpha.2026.101778 ↩

  5. Bianka Alexandra Pasat, Eleftherios Pilalis, Katarzyna Mnich, et al. MultiOmicsIntegrator: a nextflow pipeline for integrated omics analyses. Bioinformatics Advances, 2024. https://doi.org/10.1093/bioadv/vbae175 ↩

  6. Upasna Srivastava, Swarna Kanchan, Minu Kesheri, et al. Integrative omics approaches for identification of biomarkers. Integrative Omics, 2024. https://doi.org/10.1016/b978-0-443-16092-9.00010-2 ↩