Worked Analysis From Expression Profiles To Modules

From expression profiles to interpretable modules

A co-expression analysis is a sequence of choices, not a single calculation. Expression profiles are compared across samples, correlations are transformed into network weights, and the resulting network is used to group genes into modules. A module eigengene then summarizes the expression pattern of its genes, making it possible to examine associations with a measured trait. The result is useful only if the links in this chain are clear: a correlation is not a causal relationship, a weighted connection is not a physical interaction, and a module–trait association is not proof that the module drives the trait.

This tutorial works through a deliberately small, hypothetical data set so that the calculations can be checked. A real analysis would have many more genes and would require careful preprocessing, diagnostics, and stability assessment. The goal here is to make the end-to-end logic reproducible and to show exactly where assumptions enter.

Set up the profiles and assumptions

Suppose eight independent biological samples have been measured for four genes. Expression has already been transformed to a suitable scale, normalized across samples, and adjusted for known technical effects without using the trait. Each gene profile is centered across the eight samples and standardized to sample standard deviation 1. The continuous trait \(T\) is also centered and standardized. These steps put profiles on a common scale; they do not guarantee that the samples are independent or that important confounders have been removed.

To give exact, compact profiles, represent the eight sample values using six mutually orthogonal, centered contrast vectors \(q_1,\ldots,q_6\), each with sample standard deviation 1. Such vectors can be constructed from orthogonal contrasts among the eight samples. In the hypothetical data, define \(u=\sqrt{0.9}q_1+\sqrt{0.1}q_3\) and \(v=\sqrt{0.9}q_2+\sqrt{0.1}q_5\). Because the contrasts are orthogonal, both \(u\) and \(v\) have sample variance 1, and they are uncorrelated.

The four standardized gene profiles and the trait are:

Every profile has sample variance 1: for example, \(u+0.5q_4\) has variance \(1+0.25=1.25\), since \(u\) and \(q_4\) are orthogonal. The denominator standardizes it. This construction also makes the correlations transparent: \(x_1\) and \(x_2\) share \(u\), while \(x_3\) and \(x_4\) share \(v\); the two pairs use orthogonal contrasts. The samples are hypothetical, not observations from a study.

Calculate correlations and construct a weighted network

For standardized profiles, the sample correlation is the cross-product divided by \(n-1\). Here \(n=8\), and the profile variances are 1. For the first pair, the numerator corresponds to the covariance of \(u+0.5q_4\) and \(u-0.5q_4\), divided by their common variance \(1.25\). Orthogonality makes the cross-terms vanish:

By the same calculation, \(r_{34}=0.6\). Every cross-pair correlation is 0 because \(u,v,q_4,q_6\) are built from mutually orthogonal contrasts. Thus the correlation matrix has diagonal entries 1, within-pair entries 0.6, and cross-pair entries 0.

Choose an unsigned soft-thresholding rule with power \(\beta=2\): \(a_{ij}=|r_{ij}|^\beta\) for \(i\ne j\), and \(a_{ii}=0\). Setting the diagonal to zero means a gene is not treated as connected to itself. The within-pair weights are \(0.6^2=0.36\), and the cross-pair weights are \(0^2=0\). A weighted network retains graded connection strengths rather than choosing one correlation cutoff and turning every pair into a binary edge or non-edge.

The choice of unsigned adjacency treats positive and negative correlations of equal magnitude alike. A signed network instead encodes the sign in its weight, while a hard-threshold network would retain or discard each pair according to a cutoff. These alternatives answer somewhat different questions; the network definition must be reported rather than inferred from the word “co-expression.”

Identify modules from the network

In this four-gene example, the weighted network consists of two disconnected pairs: genes 1 and 2 are connected with weight 0.36, as are genes 3 and 4, and there are no cross-pair edges. A hierarchical clustering of genes using a network-derived dissimilarity, followed by a cut that separates the two pairs, therefore gives the illustrative modules \(M_1=\{1,2\}\) and \(M_2=\{3,4\}\). In a practical analysis, the particular clustering, dissimilarity, cut rule, and any merging decisions must be recorded. A module is not uniquely specified by the expression matrix alone.

This simple example makes the grouping unambiguous, but real networks usually have connections of varying strength rather than exact zeros. Choices such as the correlation definition, power, dissimilarity, and cut can change module membership. Sensitivity analysis and resampling are useful checks, but a single chosen network should not be presented as the only possible representation of the data.

Compute and interpret module eigengenes

A module eigengene is the first principal-component score of the samples in a module’s expression matrix. It is a one-dimensional summary chosen to capture the greatest possible sample variation among linear combinations of the module’s gene profiles. Its sign is arbitrary: multiplying all scores by \(-1\) leaves the principal-component direction unchanged, so the sign should be aligned consistently before comparing or interpreting eigengenes.

For \(M_1\), the two standardized profiles have the same shared component \(u\) and opposite \(q_4\) components. Their gene-by-gene correlation matrix is \(\left(\begin{smallmatrix}1&0.6\\0.6&1\end{smallmatrix}\right)\). Its leading eigenvector is \((1,1)/\sqrt{2}\), with eigenvalue 1.6; the other eigenvector \((1,-1)/\sqrt{2}\) has eigenvalue 0.4. Thus the first component combines the two profiles with equal positive weights. Their sum is proportional to \(u\), because the \(q_4\) components cancel. After standardizing the component score to sample standard deviation 1 and choosing its sign, the first eigengene is \(E_1=u\). The same reasoning gives \(E_2=v\) for \(M_2\).

Test module–trait associations

For a continuous trait and an eigengene, the sample Pearson correlation \(r\) can be tested against \(H_0:\rho=0\), where \(\rho\) is the population correlation. Under independent sampling from a bivariate normal population, the usual test statistic under \(H_0\) is \(t=r\sqrt{(n-2)/(1-r^2)}\), with \(df=n-2\). The test assesses evidence of a linear association; it does not establish causation. In this example there are eight independent samples, and we use the conventional two-sided test.

For \(E_1\) and \(T\), \(r=\sqrt{0.9}\). Substituting \(n=8\) gives:

The two-sided reference p-value is \(2P(t_6\geq 7.348)\), approximately 0.000325. For \(E_2\), the observed correlation is 0, so \(t=0\) and the two-sided p-value is 1. These p-values are calculated as if the module–trait tests were specified independently of the observed associations. Here the network and modules were formed without using the trait, an important safeguard against selecting modules precisely because they look associated.

Checks, alternatives, and limitations

A reproducible analysis should report the expression scale and preprocessing, the sample units, the correlation measure, the adjacency rule and power, the module-detection choices, the eigengene convention, the trait model, and the multiplicity adjustment. Diagnostics help reveal when the numerical pipeline may be producing an untrustworthy result:

• Inspect profiles and correlations. A single outlying sample can dominate Pearson correlations. Check sample-level quality and profile plots; consider whether a robust correlation or a rank-based measure is scientifically appropriate.
• Check technical and biological structure. Batch, cell composition, relatedness, or other shared structure can create correlations between genes that do not reflect the intended biological signal. A network built from confounded profiles can still look orderly.
• Check the module solution. Very small modules, one large catch-all module, or modules that change substantially under reasonable construction choices call for caution. A visually clean cluster is not by itself evidence of biological coherence.
• Check eigengene and test assumptions. Verify that eigengenes are calculated from the intended samples and that their sign is aligned consistently. The exact correlation t-test relies on independent sample units and a bivariate normal reference under the null; severe outliers, dependence, or strong departures from that model can make its p-value unreliable.
• Separate discovery from confirmation. If the trait is used to select, merge, or name modules and the same trait is then tested on those modules, nominal p-values can be optimistic. Use an independent cohort or a design that accounts for selection.

The main alternatives answer different questions. Signed adjacency emphasizes positively correlated genes, while unsigned adjacency also gives strong negative correlations high weights. A hard threshold is easy to interpret as an edge decision but can be sensitive to the selected cutoff; a soft-thresholded network retains graded weights but depends on the chosen power. Pearson correlation targets linear association and can be sensitive to extreme observations, whereas rank-based correlations target monotone association with different robustness properties. An eigengene compactly summarizes a module but can conceal subgroups or nonlinear patterns. These choices should follow the scientific aim and be checked for sensitivity, not selected only because they produce a preferred result.

A compact reproducible workflow

For a larger data set, the same logic can be implemented in the following order. The order helps keep network construction distinct from the trait tests and makes each choice auditable.