Why Signed Networks Are Superior for Plant Transcriptomics

WGCNA supports three distinct network adjacency topologies:

• Unsigned: $a_{ij} = |\text{cor}(x_i, x_j)|^\beta$. Treats strong positive co-expression ($+0.9$) and strong negative inverse expression ($-0.9$) identically. In plant stress cascades, this conflates stress activators with stress repressors within the same module.
• Signed (Recommended): $a_{ij} = |0.5 \times (1 + \text{cor}(x_i, x_j))|^\beta$. Maps correlations from $[-1, +1]$ to $[0, 1]$. Positively correlated genes ($+0.9 \to a_{ij} \approx 0.95^\beta$) cluster tightly, while negatively correlated or inverse regulatory relationships ($ -0.9 \to a_{ij} \approx 0.05^\beta$) are mapped near zero.
• Signed Hybrid: $a_{ij} = \text{cor}(x_i, x_j)^\beta$ if $\text{cor} > 0$ else $0$.

# Enable parallel processing across CPU cores
enableWGCNAThreads()

# Evaluate powers 1 through 20
powers <- c(1:10, seq(12, 20, by = 2))

# Execute scale-free topology analysis
sft <- pickSoftThreshold(
    expr_mat,
    powerVector = powers,
    networkType = "signed",
    RsquaredCut = 0.85,
    verbose = 5
)

fit <- sft$fitIndices

# Dual Diagnostic Plotting
par(mfrow = c(1, 2))

# Panel 1: Scale Independence (Signed R^2 fit)
plot(
    fit[, 1],
    -sign(fit[, 3]) * fit[, 2],
    xlab = "Soft-Thresholding Power (β)",
    ylab = "Scale-Free Topology Model Fit (Signed R²)",
    type = "n",
    main = "Scale Independence"
)
text(
    fit[, 1],
    -sign(fit[, 3]) * fit[, 2],
    labels = fit[, 1],
    col = "red",
    font = 2
)
abline(h = 0.85, col = "red", lty = 2)

# Panel 2: Mean Connectivity
plot(
    fit[, 1],
    fit[, 5],
    xlab = "Soft-Thresholding Power (β)",
    ylab = "Mean Node Connectivity",
    type = "b",
    pch = 19,
    col = "blue",
    main = "Mean Connectivity"
)

# Extract automated power estimate
soft_power <- sft$powerEstimate
if (is.na(soft_power)) {
    message("⚠ pickSoftThreshold returned NA. Inspect diagnostic curves before manual selection.")
} else {
    cat("✓ Suggested power from pickSoftThreshold: β =", soft_power, "\\n")
}

Diagnostic Audit When Scale-Free Fit Fails ($R^2 < 0.85$)

1. Sample Size Check: Does your experiment include fewer than 15 samples? Correlation estimates from $N < 12$ libraries are statistically noisy, which impairs scale-free curve fitting.
2. Global Batch Confounding: Run PCA on expr_mat. If PC1 accounts for $> 60\%$ of variance and separates sequencing dates rather than biological treatments, batch effects are flattening scale-free curvature.
3. Asymptotic Curve Inspection: Check whether $R^2$ plateaus at $0.80\text{--}0.84$ with a negative slope. If the curve flattens comfortably near $0.80$ while mean connectivity remains healthy ($30 \le k \le 150$), that power is scientifically justifiable.
4. Empirical Fallback Guidelines: As established by Zhang and Horvath (2005), for signed networks with $N < 20$ samples, empirical benchmark $\beta = 12$ is standard; for $N \ge 20$, $\beta = 9$ or $\beta = 6$ (unsigned) are widely validated.

# If automated powerEstimate is NA, assign justified manual power:
# soft_power <- 12 # Justified for signed network with N = 16 plant drought samples

# Programmatic validation of power parameter
stopifnot(
    length(soft_power) == 1,
    is.finite(soft_power),
    soft_power > 0
)
cat("Active Soft-Threshold Power validated: β =", soft_power, "\\n")