What methods exist to calculate RNA expression profile similarity

Question

Some of the work in our lab requires a comparison of a strain across several experimental conditions. We are looking to identify most similar experimental conditions based on the gene transcription response similarity from the cell.

While we could easily invent and create home-grown methods to do it, their implementation and testing are a laborious project in themselves and are outside the scope of our current work.

Are there any methods for RNA expression profile similarity computation that has been already published? If yes, what is your experience using them?

Answer 1

The suggestion to use proportionality is probably the correct one if you are interested in similar patterns between samples. But not if you are interested in absolute differences.

For example: the following two samples are similar in pattern, but not similar on absolute levels:

        Sample 1    Sample 2    Sample 3
Gene A    10          100           80
Gene B     8           80          100  
Gene C    12          120          120

Samples 1 and 2 have a perfect proportionality (phi is 0) and also a perfect correlation (as a side note any pair with a perfect proportionality will always have a perfect correlation). However, in terms of logfold changes, samples 2 and 3 are more like each other.

Of course in real life you would never see a comparison like the sample 1 - sample 2 one because normalisation would have removed the scale difference. This was exactly the point brought up by the proportionality paper. But normalisation methods don't normally guarantee that the sum of expression for each sample is identical, and such differences can still occur.

An alternative that might suit more in the second case is either the euclidean distance between the samples, or the euclidean distance on the first two components of a principle component or multi-dimensional scaling. The latter is effectively using mean logFC between samples.

Assuming that x is a matrix containing normalised, log transformed expression values, you could use R and limma to calculate distance in multi-dimensional scaled space as follows:

library(limma)
mds <- plotMDS(x, plot=FALSE)
mds <- data.frame(mds$x, mds$y)
distances <- dist(mds)

Answer 2

There's a new paper that's just appeared on the subject of "proportionality", including a method by which RNA expression might be compared.

This is a new concept for me, and the article is not easy enough for me to read that I can write up a quick summary; the authors don't seem to devote a section in the paper to defining "proportionality". However, here's an interesting chunk from the article:

We graphed the network of relationships between these mRNAs (S5 Fig.), an approach similar to gene co-expression network [12] or weighted gene co-expression analysis [13] but founded on proportionality and therefore valid for relative data. The network revealed one cluster of 96, and many other smaller clusters of mRNAs behaving proportionally across conditions.
...
We are also keen to raise awareness that correlation (and other statistical methods that assume measurements come from real coordinate space) should not be applied to relative abundances. This is highly relevant to gene coexpression networks [12]. Correlation is at the heart of methods like Weighted Gene Co-expression Network Analysis [13] and heatmap visualization [14]. These methods are potentially misleading if applied to relative data.

Answer 3

If you have a lot of samples with very different environmental conditions, then a Weighted Gene Correlation Network Analysis (WGCNA) might be appropriate.

This type of analysis looks for genes that trace similar (or opposing) expression patterns throughout different conditions (e.g. high-medium-medium-low-absent-high would be highly negatively correlated with low-medium-medium-high-very high-low).

That particular paper introduces the concept of "modules", which are groups of genes that share similar expression patterns. Functions are available for plotting how the expression of canonical module members changes throughout the different conditions, and for identifying which module (or modules) a particular gene is likely to be a member of.

WGCNA works best when there's a lot of different expression changes in the different conditions, which sounds like it would fit well with your project. However, it concentrates more on the genes, rather than the conditions (which seems like it would be less useful for you).

Answer 4

The euclidean distance is probably the simplest, both conceptually and in terms of implementation. It's by no means an elegant solution, and it may not perform well in certain circumstances.

Euclidean distance is easiest to conceptualize as the distance between two points in a two dimensional space.

Y
^
|
|
|
|            * p = (3, 3)
|
|
|
|        * q = (2, 1)
|
----------------------------> X

In this example, the distance between the two points is

d(p, q) = sqrt( (p_x-q_x)^2 + (p_y-q_y)^2 )
        = sqrt( (2-3)^2 + (1-3)^2 )
        = sqrt(5)
        ≈ 2.24

For a gene expression profile with two genes, this is exactly how the euclidean distance would be calculated, using expression values from one gene as the X axis and expression values from the other gene as the Y axis. Realistically, though, gene expression profiles typically contain thousands or tens of thousands of genes, so instead we use the generalization of the distance calculation for N dimensions.

d(p, q) = sqrt( (p_1-q_1)^2 + (p_2-q_2)^2 + ... + (p_N-q_N)^2 )

Packages for R and Python make these types of calculations trivial once you have the data loaded into the correct data structure. See Ian's answer for some example R code.