Biological network construction and interpretation

Practical 11: Genomics, systems biology and multi-omics

Learning objectives

Define a correlation edge; construct a symmetric adjacency matrix; calculate degree and inspect threshold sensitivity; separate association from regulation.

Biological question and data

Which hypothetical proteins show correlated abundance across the same animals?

The complete log2 protein table contains a supplied shared latent factor and independent noise. A correlation network is descriptive. Animal observations are paired across proteins, and edges are not directional or causal.

No biological download is required. Numerical inputs are hypothetical, not observed animal data. The runnable examples use base R. Count and abundance units are specified where used.

1. Align and calculate associations

Use animals as the observations and protein features as variables.

helper_dir <- if (dir.exists("R")) "R" else "../R"
stopifnot(file.exists(file.path(helper_dir,"teaching_data.R")))
source(file.path(helper_dir,"teaching_data.R"))
d <- teaching_data()
X <- t(d$protein)
stopifnot(identical(rownames(X),d$metadata$animal),!anyNA(X))
C <- cor(X)
round(C,2)

          protein01 protein02 protein03 protein04 protein05 protein06 protein07
protein01      1.00      0.09      0.23      0.24      0.07      0.08      0.25
protein02      0.09      1.00      0.49      0.42      0.42      0.48      0.61
protein03      0.23      0.49      1.00      0.66      0.60      0.73      0.71
protein04      0.24      0.42      0.66      1.00      0.57      0.61      0.73
protein05      0.07      0.42      0.60      0.57      1.00      0.52      0.59
protein06      0.08      0.48      0.73      0.61      0.52      1.00      0.77
protein07      0.25      0.61      0.71      0.73      0.59      0.77      1.00
protein08      0.24      0.59      0.74      0.68      0.70      0.78      0.85
          protein08
protein01      0.24
protein02      0.59
protein03      0.74
protein04      0.68
protein05      0.70
protein06      0.78
protein07      0.85
protein08      1.00

stopifnot(isSymmetric(C),all(abs(diag(C)-1)<1e-12))

Positive correlations are expected from the shared generating factor. The same pattern could arise from a shared technical factor in another dataset.

2. Define edges and calculate degree

Use a supplied absolute-correlation cutoff of .55. Remove self-edges and count each undirected edge once.

cutoff <- .55
A <- abs(C)>=cutoff; diag(A) <- FALSE
edges <- which(A & upper.tri(A),arr.ind=TRUE)
degree <- rowSums(A)
data.frame(protein=colnames(X),degree)

            protein degree
protein01 protein01      0
protein02 protein02      2
protein03 protein03      5
protein04 protein04      5
protein05 protein05      4
protein06 protein06      4
protein07 protein07      6
protein08 protein08      6

stopifnot(isSymmetric(A),sum(degree)==2*nrow(edges))
angle <- seq(0,2*pi,length.out=ncol(X)+1)[-(ncol(X)+1)]
pos <- cbind(cos(angle),sin(angle))
plot(pos,type="n",axes=FALSE,xlab="",ylab="",asp=1,
 xlim=c(-1.35,1.35),ylim=c(-1.35,1.35),main="Correlation network")
if(nrow(edges)>0) for(i in seq_len(nrow(edges))) {
 e <- edges[i,];segments(pos[e[1],1],pos[e[1],2],pos[e[2],1],pos[e[2],2],col="grey60")
}
points(pos,pch=19,col="#1b9e77",cex=1.5)
text(pos,labels=colnames(X),pos=3,cex=.7)

The picture is a layout of a defined graph, not a spatial map of interacting proteins.

3. Examine threshold dependence

Use the same variables and animals for each threshold.

threshold <- c(.3,.55,.75)
data.frame(threshold,edges=vapply(threshold,function(t)
 sum((abs(C)>=t)[upper.tri(C)]),integer(1)))

  threshold edges
1      0.30    21
2      0.55    16
3      0.75     3

Connectivity changes with the chosen threshold; a hub designation should state that dependence.

4. Generate a common-driver graph

Two features can correlate because both respond to a third variable. The following hypothetical example contains no direct connection between the two features.

set.seed(1103)
common <- rnorm(200)
feature_A <- common+rnorm(200,sd=.3)
feature_B <- common+rnorm(200,sd=.3)
c(marginal=cor(feature_A,feature_B),
 after_common_adjustment=cor(residuals(lm(feature_A~common)),
 residuals(lm(feature_B~common))))

               marginal after_common_adjustment 
             0.91972604              0.02861018 

The residual association need not be exactly zero. Adjustment here uses the known generator; in real data the relevant common drivers may be unobserved or measured with error.

5. Describe edge meaning and stability

Record the correlation measure, feature transformation, animal subset, threshold and missingness policy. A graph made from pairwise complete observations can have edges estimated from different animal subsets. That can change comparability.

Assess sampling sensitivity by resampling independent animals, then inspect threshold sensitivity separately. A stable edge remains an association under the stated procedure. It does not establish regulatory direction or rule out a stable batch effect.

Interpretation challenge

Why can a node linked to every other feature be a poor intervention target?

Suggested reasoning

Its degree may reflect a shared driver, technical structure or duplicated measurements. Connectivity under a correlation rule is not evidence that perturbing the node changes the others.

Interpretation and limits

Correlation can reflect diet, genotype, tissue mixture, batch or other common drivers. A network inferred from abundance is not equivalent to a curated interaction database. The threshold is supplied rather than selected to optimize a biological story; stability is assessed in the biomarker practical.

Exercises

1. Why is the diagonal removed?
2. Could a batch shift create an apparent module?
3. Why is the sum of degrees twice the edge count?

Suggested answers

1. A feature’s self-correlation is not an inter-feature relationship.
2. Yes; shared technical variation can correlate many measurements.
3. Each undirected edge contributes one degree to each endpoint.

Continue learning

B14 Pathways and networks.

Continue with network and biomarker validation. Return to the practical index or the unified map.