Gene Expression Data

The expression file must contain a column with gene identifiers (Gene Symbol, Entrez ID, Ensembl, etc.) and at least one column with expression values (e.g., normalized counts, log2 fold-change).

If some network genes have no expression value, levi assigns a neutral value (0.5), representing no change, and generates a log file listing the missing genes.

Expression data can be obtained from:

• Gene Expression Omnibus (GEO)
• Array Express
• The Cancer Genome Atlas (TCGA)
• Sequence Read Archive (SRA)

Biological Network

levi supports the following network file formats:

Network interaction data can be obtained from:

• STRING database
• miRBase

Networks with no coordinates

Every format above carries a position for each node, because levi draws its landscape over the plane the network occupies. Interaction databases and collaborators usually hand over something simpler: a list of pairs, with no coordinates at all.

leviFromEdges() computes them with igraph (Csárdi and Nepusz 2006) and returns the two tables levi() expects:

library(levi)

interactions <- data.frame(
    from = c("HUB", "HUB", "HUB", "HUB", "N1", "N2", "N3", "N4"),
    to   = c("N1", "N2", "N3", "N4", "N5", "N6", "N7", "N8"),
    stringsAsFactors = FALSE)

set.seed(42)                       # force-directed layouts are stochastic
net <- leviFromEdges(interactions, layout = "kk")

## Network laid out with 'kk': 9 nodes, 8 edges.

head(net$nodes)

##   name        x         y
## 1  HUB 50.50000  50.50000
## 2   N1 38.79074  76.13341
## 3   N2 24.86659  38.79074
## 4   N3 62.20926  24.86659
## 5   N4 76.13341  62.20926
## 6   N5 27.88855 100.00000

expression <- data.frame(
    ID      = c("HUB", paste0("N", 1:8)),
    Test    = c(200, 200, 200, 200, 200, 5, 5, 5, 5),
    Control = c(10, 10, 10, 10, 10, 200, 200, 200, 200))

res_edges <- levi(
    networkCoordinatesInput  = net$nodes,
    networkInteractionsInput = net$edges,
    fileTypeInput            = "stg",
    expressionInput          = expression,
    geneSymbolInput          = "ID",
    readExpColumn            = readExpColumn("Test-Control"),
    resolutionValueInput     = 30,
    smoothValueInput         = 50)

res_edges$scores[, c("Gene", "LandscapeScore", "Rank")]

##   Gene LandscapeScore Rank
## 1  HUB         0.9524    1
## 2   N1         0.9027    2
## 3   N2         0.9027    3
## 4   N3         0.9027    4
## 5   N4         0.9027    5
## 6   N5         0.0854    6
## 7   N6         0.0854    7
## 8   N7         0.0854    8
## 9   N8         0.0854    9

A caveat worth keeping in mind: the layout is an analytical choice, not only an aesthetic one. The same expression values arranged differently give a different landscape, because levi reads neighbourhoods. Where a curated layout already exists, prefer it over a computed one.

Graphical User Interface (GUI)

The GUI runs locally via Shiny and covers the same analysis as script mode: it calls levi() with the parameters chosen in the side panel, so a landscape built in the interface and one built from a script agree, p-values included.

library(levi)
LEVIui(browser = TRUE)   # Opens in the system browser
LEVIui(browser = FALSE)  # Opens inside RStudio

Results

The main panel shows the figure in two tabs and the tables in four.

• 2D landscape: the map, with a format selector (TIFF, BMP, JPEG, PNG) and Download Plot. Click and drag on the map to select an area: Expression area shows the summed score of the selected cells and the Genes tab lists the genes under the selection.
• 3D surface: the interactive surface, shown when 3D surface is ticked. Rotate it, then use the camera icon of the toolbar to save a PNG of the current view, or Download 3D (HTML) to keep the view in an interactive file. The camera of the current view is printed as R code, ready to paste into leviSave3D(camera = ...) to reproduce it from a script.
• Genes: the genes under the area selected on the 2D map. This tab is shown only while the 2D map is in front.
• Node scores: the landscape score of every node, ranked, with a CSV download. This is result$scores of script mode.
• Peaks and valleys: the detected extrema and the closest gene to each. This is result$peaks.
• Regions: the connected areas of the landscape beyond the neutral score, with their size, mass and peak score; after a permutation test, the regional p-value and whether the region passes the significance level. This is result$regions$summary.

Figure 2: The Settings tab of the side panel, the 3D surface and the regions table after a permutation test with 100 iterations.

A busy indicator is shown at the top of the page while the landscape or the permutation test is being computed.

Script Mode

library(levi)

template_network <- file.path(system.file(package = "levi"), "extdata",
                              "medusa.dat", fsep = .Platform$file.sep)

template_expression <- file.path(system.file(package = "levi"), "extdata",
                                 "expression.dat", fsep = .Platform$file.sep)

multicolor <- levi(
    networkCoordinatesInput = template_network,
    expressionInput        = template_expression,
    fileTypeInput          = "dat",
    geneSymbolInput        = "ID",
    readExpColumn          = readExpColumn("TumorCurrentSmoker-NormalNeverSmoker"),
    contrastValueInput     = 50,
    resolutionValueInput   = 20,
    zoomValueInput         = 50,
    smoothValueInput       = 5,
    contourLevi            = TRUE
)

## There are 1 nodes without expression value, see log in path: /tmp/Rtmp9k2SHA/TumorCurrentSmoker-NormalNeverSmoker/levi.log

base <- readExpColumn(
    "TumorCurrentSmoker-NormalNeverSmoker",
    "TumorFormerSmoker-NormalFormerSmoker"
)

levi(
    networkCoordinatesInput = template_network,
    expressionInput         = template_expression,
    fileTypeInput           = "dat",
    geneSymbolInput         = "ID",
    readExpColumn           = base,
    setcolor                = "pink_green",
    contourLevi             = FALSE
)

## There are 1 nodes without expression value, see log in path: /tmp/Rtmp9k2SHA/TumorCurrentSmoker-NormalNeverSmoker/levi.log

## There are 1 nodes without expression value, see log in path: /tmp/Rtmp9k2SHA/TumorFormerSmoker-NormalFormerSmoker/levi.log

# Microarray RMA — keep log2 scale, use logfc mode
res_rma <- levi(
    networkCoordinatesInput = template_network,
    expressionInput         = rma_expression_df,   # columns: ID, Tumor_log2, Normal_log2
    fileTypeInput           = "dat",
    geneSymbolInput         = "ID",
    readExpColumn           = readExpColumn("Tumor_log2-Normal_log2"),
    signal_mode             = "logfc",
    expressionLog           = FALSE   # keep in log2 scale; no 2^ transform
)

# scRNA-seq — single-column avg_log2FC with reduced steepness
res_sc <- levi(
    networkCoordinatesInput = template_network,
    expressionInput         = seurat_markers_df,  # column: GeneID, avg_log2FC
    fileTypeInput           = "dat",
    geneSymbolInput         = "GeneID",
    readExpColumn           = readExpColumn("avg_log2FC-avg_log2FC"),
    signal_mode             = "logfc",
    logfc_k                 = 0.5     # softer gradient for large FC values
)

# DESeq2 log2FoldChange as single column
res_deseq <- levi(
    networkCoordinatesInput = template_network,
    expressionInput         = leviFromDESeq2(dds_result),
    fileTypeInput           = "dat",
    geneSymbolInput         = "GeneID",
    readExpColumn           = readExpColumn("log2FoldChange-log2FoldChange"),
    signal_mode             = "logfc"
)

logfc_net  <- system.file("extdata", "logfc_network.dat", package = "levi")
logfc_expr <- system.file("extdata", "logfc_expression.dat", package = "levi")

# Control is 8 for every gene and Test runs from 3 to 10, so the logFC
# spans -5 to +2 on the log2 scale.
res_logfc <- levi(
    networkCoordinatesInput = logfc_net,
    expressionInput         = logfc_expr,
    fileTypeInput           = "dat",
    geneSymbolInput         = "ID",
    readExpColumn           = readExpColumn("Test-Control"),
    resolutionValueInput    = 40,
    signal_mode             = "logfc",
    expressionLog           = FALSE   # already log2; no back-transform
)

res_logfc$scores[, c("Gene", "LandscapeScore", "Rank")]

##   Gene LandscapeScore Rank
## 1   GF         0.8230    1
## 2   GE         0.3309    2
## 3   GD         0.1008    3
## 4   GC         0.0663    4
## 5   GB         0.0183    5
## 6   GA         0.0077    6

From DESeq2

library(DESeq2)
library(levi)

# After running DESeq2 analysis:
res <- results(dds, contrast = c("condition", "Tumor", "Normal"))

expr_df <- leviFromDESeq2(res, gene_col = "GeneSymbol")

levi(
    expressionInput         = expr_df,
    networkCoordinatesInput = template_network,
    fileTypeInput           = "dat",
    geneSymbolInput         = "GeneSymbol",
    readExpColumn           = readExpColumn("log2FoldChange-log2FoldChange"),
    signal_mode             = "logfc"
)

From edgeR

library(edgeR)
library(levi)

fit  <- glmQLFit(dge, design)
qlf  <- glmQLFTest(fit, coef = 2)

expr_df <- leviFromEdgeR(qlf, gene_col = "GeneSymbol")

levi(
    expressionInput         = expr_df,
    networkCoordinatesInput = template_network,
    fileTypeInput           = "dat",
    geneSymbolInput         = "GeneSymbol",
    readExpColumn           = readExpColumn("logFC-logFC"),
    signal_mode             = "logfc"
)

From limma

library(limma)
library(levi)

fit2    <- eBayes(fit)
expr_df <- leviFromLimma(fit2, coef = 1, gene_col = "GeneSymbol")

levi(
    expressionInput         = expr_df,
    networkCoordinatesInput = template_network,
    fileTypeInput           = "dat",
    geneSymbolInput         = "GeneSymbol",
    readExpColumn           = readExpColumn("logFC-logFC"),
    signal_mode             = "logfc"
)

Checking the adapters without the tools

Each adapter also accepts a plain data.frame carrying the columns the corresponding tool produces. That makes the conversion reproducible here, and is handy for testing a pipeline before installing the full stack:

# A DESeq2 results table has baseMean and log2FoldChange.
res_de <- data.frame(
    baseMean       = c(1200, 1100, 900),
    log2FoldChange = c(4.3, 4.1, -5.3),
    row.names      = c("HUB", "N1", "N5"))
leviFromDESeq2(res_de, gene_col = "GeneID")

##   GeneID baseMean log2FoldChange
## 1    HUB     1200            4.3
## 2     N1     1100            4.1
## 3     N5      900           -5.3

# An edgeR topTags table has logFC and logCPM.
tt_edger <- data.frame(
    logFC     = c(4.3, -5.1),
    logCPM    = c(10.2, 9.8),
    row.names = c("HUB", "N5"))
leviFromEdgeR(tt_edger, gene_col = "GeneID")

##   GeneID logCPM logFC
## 1    HUB   10.2   4.3
## 2     N5    9.8  -5.1

# A limma topTable has logFC and AveExpr.
tt_limma <- data.frame(
    logFC     = c(4.3, -5.1),
    AveExpr   = c(8.1, 7.7),
    row.names = c("HUB", "N5"))
leviFromLimma(tt_limma, gene_col = "GeneID")

##   GeneID AveExpr logFC
## 1    HUB     8.1   4.3
## 2     N5     7.7  -5.1

# A Seurat FindMarkers table has avg_log2FC and pct.2, which become
# Test and Control.
mk_seurat <- data.frame(
    avg_log2FC = c(2.5, -1.8),
    pct.2      = c(0.3, 0.8),
    row.names  = c("HUB", "N5"))
leviFromSeurat(mk_seurat, gene_col = "GeneID")

##   GeneID Control Test
## 1    HUB     0.3  2.5
## 2     N5     0.8 -1.8

Node landscape scores

Every levi() call returns a ranked table of landscape scores invisibly. A score close to 1 means the gene’s network neighborhood is strongly up-regulated in the test condition; close to 0 means strongly down-regulated.

The network weights the landscape by degree. Each edge adds a support point at its midpoint carrying the mean of its two endpoints, so a hub with twenty interactions surrounds itself with twenty extra points while a leaf adds one. The smoothed surface around a hub is therefore dominated by the hub and its neighbours, and a moderately changed hub can look more prominent than a strongly changed peripheral gene. The permutation tests keep the network fixed, so this weighting is part of their null and does not bias the p-values; it does affect what the figure emphasises, and result$scores should be read with the node’s degree in mind.

result <- levi(
    networkCoordinatesInput = template_network,
    expressionInput         = template_expression,
    fileTypeInput           = "dat",
    geneSymbolInput         = "ID",
    readExpColumn           = readExpColumn("TumorCurrentSmoker-NormalNeverSmoker"),
    contrastValueInput      = 50,
    resolutionValueInput    = 20,
    zoomValueInput          = 50,
    smoothValueInput        = 5
)

## There are 1 nodes without expression value, see log in path: /tmp/Rtmp9k2SHA/TumorCurrentSmoker-NormalNeverSmoker/levi.log

# Top 10 most altered genes in the network
head(result$scores, 10)

##     Gene         X         Y LandscapeScore Rank
## 1  OGDHL 0.3323616 0.1683213         0.7768    1
## 2   IDH1 0.2850717 0.1341222         0.6944    2
## 3   IDH2 0.2352682 0.1331376         0.6678    3
## 4   DLAT 0.1228036 0.3618685         0.6486    4
## 5     FH 0.3917084 0.2269471         0.6327    5
## 6  PDHA2 0.3384309 0.2020927         0.6325    6
## 7   SDHC 0.2797636 0.4590206         0.6221    7
## 8   MDH1 0.3203979 0.2105595         0.5974    8
## 9  PDHA1 0.2585272 0.1312074         0.5960    9
## 10 IDH3A 0.2843243 0.3034554         0.5852   10

The scores table contains columns Gene, X, Y, LandscapeScore (0–1), and Rank. Genes with scores > 0.7 are candidates for up-regulated hub regions; genes with scores < 0.3 are down-regulated hubs.

Automatic peak detection

Peaks (local maxima) and valleys (local minima) are detected automatically and returned in result$peaks. Each row reports the nearest gene, the matrix position, and the landscape score.

# Detected peaks (up-regulated regions) and valleys (down-regulated)
result$peaks

##   Type NearestGene MatrixRow MatrixCol  Score
## 1 peak       OGDHL        48        16 0.7910
## 2 peak        MDH1        42        25 0.6726
## 3 peak          CS        23        30 0.6705

Permutation significance test

To distinguish visually prominent features from statistical noise, use n_perm to run a permutation test. Measured expression values are randomly shuffled among network nodes, edge signals are recalculated and the landscape is rebuilt each time. The layout remains fixed.

By default (inference_unit = "region") the eight-connected regions beyond the neutral score are redetected in every permutation and each observed region is compared with the largest regional mass seen under the null (Westfall and Young 1993; Nichols and Holmes 2002), over both directions. result$regions$summary gains PSpatial and Significant, and significant regions are outlined in white on the plot. The legacy inference_unit = "cell" tests every grid cell and draws contours from BY-adjusted p-values; it is what the graphical interface uses. The vignette levi_inference explains which null answers which question.

result_sig <- levi(
    networkCoordinatesInput = template_network,
    expressionInput         = template_expression,
    fileTypeInput           = "dat",
    geneSymbolInput         = "ID",
    readExpColumn           = readExpColumn("TumorCurrentSmoker-NormalNeverSmoker"),
    n_perm    = 100,     # more permutations improve p-value resolution
    sig_level = 0.05,
    perm_side = "both"   # "over", "under", or "both" (default)
)

## There are 1 nodes without expression value, see log in path: /tmp/Rtmp9k2SHA/TumorCurrentSmoker-NormalNeverSmoker/levi.log

## Permutation test: 100 iterations for 'TumorCurrentSmoker-NormalNeverSmoker'...

# One row per region, with its permutation p-value:
result_sig$regions$summary[, c("Region", "Direction", "Cells", "Mass",
                               "PSpatial", "Significant")]

##    Region Direction Cells         Mass  PSpatial Significant
## 1 over_01      over  1137 0.0059606740 0.2574257       FALSE
## 2 over_03      over   288 0.0005096155 0.9405941       FALSE
## 3 over_02      over   107 0.0001561263 0.9702970       FALSE

With inference_unit = "cell", result$pvalues holds the $over and $under adjusted p-value matrices instead, and the landscape overlays two contours:

Use perm_side = "both" (default) to show both simultaneously.

Interpretation: genes inside the dashed contour are strong candidates for up-regulated pathway hubs; genes inside the dotted contour are down-regulated hubs. Both boundaries are calibrated against a null distribution from random label shuffling, controlling for network topology effects.

Beyond the landscape test: replication and graph TFCE

The node-label p-value above is conditional on the network and on its layout. It asks whether the observed arrangement of the same values is unusually clustered, so it says nothing about replication, and a different drawing of the same graph is a different test. Two consequences matter in practice:

• When most of a network responds (an a-priori pathway under a stimulus that hits all of it, for instance), any shuffling of the values produces regions of similar mass and the node-label p-value approaches 1, even though the effect is real and strong. The test is answering a spatial question, not a biological one.
• The landscape depends on the layout. Force-directed layouts are stochastic, and regions detected on one drawing overlap only partly with those detected on another.

When biological replicates are available, we recommend reporting the landscape together with two sample-label tests, which permute condition labels across replicates instead of values across nodes:

• leviReplicateInference() keeps the landscape and its regions but recomputes the fold-change from permuted labels, so PSpatial now answers whether the region replicates.
• leviGraphTFCEInference() needs no layout at all. It fits a limma moderated t per gene (Smyth 2004) and integrates it over every threshold with threshold-free cluster enhancement (Smith and Nichols 2009) on the network components, controlling the family-wise error rate with the maximum over genes and directions. It is the more powerful of the two and is unaffected by the layout, so it is the natural confirmatory test for a region seen on the landscape.

With four replicates per group the example below is small enough to enumerate all choose(8, 4) = 70 label arrangements, so the smallest attainable p-value is 1/70.

hub_net <- system.file("extdata", "hub_network.dat", package = "levi")
genes   <- c("HUB", paste0("N", 1:8))

# Four control and four treated replicates on a log2 scale; the hub and its
# first two neighbours respond.
set.seed(2026)
expression <- matrix(rnorm(9 * 8, mean = 6, sd = 0.3), 9, 8,
                     dimnames = list(genes, paste0("s", 1:8)))
groups <- rep(c("control", "treated"), each = 4)
expression[c("HUB", "N1", "N2"), groups == "treated"] <-
    expression[c("HUB", "N1", "N2"), groups == "treated"] + 1.5

# Landscape regions tested against biological replication
rep_res <- leviReplicateInference(
    expression, groups, test = "treated", control = "control",
    networkCoordinatesInput = hub_net, fileTypeInput = "dat",
    resolutionValueInput = 20, smoothValueInput = 30, seed = 1)
rep_res$regions$summary[, c("Region", "Direction", "Mass", "PSpatial",
                            "Significant")]

##    Region Direction        Mass   PSpatial Significant
## 1 over_03      over 0.033651593 0.02857143        TRUE
## 2 over_05      over 0.003513050 0.48571429       FALSE
## 3 over_04      over 0.002689098 0.48571429       FALSE
## 4 over_02      over 0.002476661 0.48571429       FALSE
## 5 over_01      over 0.002083346 0.48571429       FALSE

# Layout-free confirmation: graph TFCE with the same null
tfce <- leviGraphTFCEInference(
    expression, groups, hub_net, test = "treated", control = "control",
    seed = 1)
tfce$statistic[order(tfce$statistic$PGlobal), ][1:4, ]

##     Gene        T        TFCE      POver PUnder    PGlobal GlobalSignificant
## HUB  HUB 8.219189 264.0536523 0.01428571      1 0.01428571              TRUE
## N1    N1 6.796803 180.1085001 0.01428571      1 0.02857143              TRUE
## N2    N2 6.741273 180.1085001 0.01428571      1 0.02857143              TRUE
## N4    N4 1.009317   0.7337502 0.58571429      1 0.80000000             FALSE

Both tests find the hub and its two responding neighbours; the node-label test alone could not have told us that this replicates. For paired designs pass blocks; for single-cell data use the leviSingleCell* functions, which permute donors rather than cells. The vignette levi_inference explains the three nulls in detail and levi_validation reports their calibration and power.

Available datasets

Bimodal dataset

Two distinct hubs (over-expressed cluster A and under-expressed cluster B) connected by a neutral bridge gene. The expected landscape shows a peak on the left and a valley on the right — ideal for validating leviDiff.

bimodal_net  <- file.path(system.file(package="levi"), "extdata",
                           "bimodal_network.dat")
bimodal_expr <- file.path(system.file(package="levi"), "extdata",
                           "bimodal_expression.dat")

res_bimodal <- levi(
    networkCoordinatesInput = bimodal_net,
    expressionInput         = bimodal_expr,
    fileTypeInput           = "dat",
    geneSymbolInput         = "ID",
    readExpColumn           = readExpColumn("Test-Control"),
    contrastValueInput      = 50,
    resolutionValueInput    = 20,
    zoomValueInput          = 50,
    smoothValueInput        = 5,
    contourLevi             = TRUE
)

## Warning: `stat_contour()`: Zero contours were generated

## Warning in min(x): no non-missing arguments to min; returning Inf

## Warning in max(x): no non-missing arguments to max; returning -Inf

signal_mode = “logfc” with log2 microarray intensities

The logfc_expression dataset simulates RMA-normalized microarray data (log2 intensities ranging 5–8), where signal_mode = "logfc" produces a wider score range than signal_mode = "ratio" because it captures the direction and magnitude of fold-change explicitly.

logfc_net  <- file.path(system.file(package="levi"), "extdata",
                         "logfc_network.dat")
logfc_expr <- file.path(system.file(package="levi"), "extdata",
                         "logfc_expression.dat")

res_logfc <- levi(
    networkCoordinatesInput = logfc_net,
    expressionInput         = logfc_expr,
    fileTypeInput           = "dat",
    geneSymbolInput         = "ID",
    readExpColumn           = readExpColumn("Test-Control"),
    signal_mode             = "logfc",
    contrastValueInput      = 50,
    resolutionValueInput    = 20,
    zoomValueInput          = 50,
    smoothValueInput        = 5,
    contourLevi             = TRUE
)

cat("Score range with logfc mode:",
    round(diff(range(res_logfc$scores$LandscapeScore)), 3), "\n")

## Score range with logfc mode: 0.874

Multi-comparison (batch mode) and leviDiff

The hub_multicomp_expression dataset has three conditions for the hub network. Using batch mode produces one landscape per comparison; leviDiff() subtracts them cell-by-cell to reveal regional changes between conditions.

hub_net <- file.path(system.file(package="levi"), "extdata", "hub_network.dat")
mc_expr <- file.path(system.file(package="levi"), "extdata",
                     "hub_multicomp_expression.dat")

# Batch mode: two comparisons in one call
res_multi <- levi(
    networkCoordinatesInput = hub_net,
    expressionInput         = mc_expr,
    fileTypeInput           = "dat",
    geneSymbolInput         = "ID",
    readExpColumn           = readExpColumn("Cond_A-Cond_B", "Cond_A-Cond_C"),
    resolutionValueInput    = 20,
    smoothValueInput        = 5
)

# Side-by-side
leviGrid(res_multi)

# Regional difference between conditions B and C
leviDiff(res_multi[[1]], res_multi[[2]],
         label_a = "Cond_A vs Cond_B",
         label_b = "Cond_A vs Cond_C")

RNA-seq: airway (DESeq2 + STRING)

The airway dataset (Himes et al. 2014) contains RNA-seq counts from airway smooth muscle cells treated with dexamethasone (DEX) vs untreated controls across 4 cell lines. DESeq2 identifies differentially expressed genes; leviFromSTRING() builds the interaction network automatically.

BiocManager::install(c("airway", "DESeq2", "STRINGdb"))
library(airway); library(DESeq2); library(levi)

data(airway)
dds <- DESeqDataSet(airway, design = ~cell + dex)
dds <- DESeq(dds)
res <- results(dds, contrast = c("dex", "trt", "untrt"))

# Convert to levi input
expr_df <- leviFromDESeq2(res)

# Top 80 DE genes — build STRING network
top80 <- head(expr_df$GeneID[order(abs(expr_df$log2FoldChange),
                                    decreasing = TRUE)], 80)
set.seed(42)
net <- leviFromSTRING(top80, species = 9606, score_threshold = 400)

# Landscape: which network hubs respond to DEX treatment?
levi(
    expressionInput          = expr_df,
    networkCoordinatesInput  = net$nodes,
    networkInteractionsInput = net$edges,
    fileTypeInput            = "stg",
    geneSymbolInput          = "GeneID",
    readExpColumn            = readExpColumn("log2FoldChange-log2FoldChange"),
    signal_mode              = "logfc",
    n_perm                   = 500,
    perm_side                = "both"
)

Microarray: ALL (limma + STRING)

The ALL dataset (Chiaretti et al. 2004) contains Affymetrix HG-U95Av2 microarray data from 128 Acute Lymphoblastic Leukemia patients. RMA-normalized values are in log2 scale, making signal_mode = "logfc" the natural choice.

BiocManager::install(c("ALL", "limma", "STRINGdb"))
library(ALL); library(limma); library(levi)

data(ALL)

# B-cell vs T-cell subtypes
design  <- model.matrix(~0 + ALL$BT)
colnames(design) <- c("B", "T")
contrast <- makeContrasts(B - T, levels = design)
fit  <- lmFit(ALL, design)
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)

expr_df <- leviFromLimma(fit2, coef = 1)

set.seed(7)
net <- leviFromSTRING(expr_df$GeneID, species = 9606,
                      score_threshold = 700,   # high confidence
                      layout = "kk")           # Kamada-Kawai for <100 nodes

levi(
    expressionInput          = expr_df,
    networkCoordinatesInput  = net$nodes,
    networkInteractionsInput = net$edges,
    fileTypeInput            = "stg",
    geneSymbolInput          = "GeneID",
    readExpColumn            = readExpColumn("logFC-logFC"),
    signal_mode              = "logfc",
    logfc_k                  = 2        # tighter FC range in microarray
)

scRNA-seq: PBMC 10x (Seurat + STRING)

Single-cell RNA-seq data from 3k PBMCs of a healthy donor (10x Genomics 2016). FindMarkers() computes avg_log2FC values per gene for a given cell type contrast, which are fed directly into levi using single-column logFC mode.

BiocManager::install(c("TENxPBMCData", "STRINGdb"))
install.packages("Seurat")
library(TENxPBMCData); library(Seurat); library(levi)

pbmc_sce <- TENxPBMCData("pbmc3k")
pbmc <- as.Seurat(pbmc_sce)

# Standard Seurat pipeline
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc)
pbmc <- FindNeighbors(pbmc)
pbmc <- FindClusters(pbmc, resolution = 0.5)

# Markers between two clusters (e.g., T cells vs B cells after annotation)
markers <- FindMarkers(pbmc,
    ident.1 = "CD4 T cells",
    ident.2 = "B cells",
    min.pct = 0.25)

expr_df <- leviFromSeurat(markers, gene_col = "row.names")

set.seed(21)
net <- leviFromSTRING(
    rownames(markers), species = 9606,
    score_threshold = 400,
    layout = "fr"
)

levi(
    expressionInput          = expr_df,
    networkCoordinatesInput  = net$nodes,
    networkInteractionsInput = net$edges,
    fileTypeInput            = "stg",
    geneSymbolInput          = "GeneID",
    readExpColumn            = readExpColumn("avg_log2FC-avg_log2FC"),
    signal_mode              = "logfc",
    logfc_k                  = 0.5      # softer gradient for large FC values
)

Saving and reusing a downloaded network

leviFromSTRING() returns a plain R list, so the network can be saved to disk and reloaded in future sessions — avoiding repeated downloads.

# --- Download once ---
set.seed(42)
net <- leviFromSTRING(mapk_genes, species = 9606, score_threshold = 400)

# Option 1: save as TSV files — readable by levi() directly via file path
write.table(net$nodes, "string_nodes.tsv",
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(net$edges, "string_edges.tsv",
            sep = "\t", row.names = FALSE, quote = FALSE)

# Option 2: save the full R object (preserves igraph + coordinates)
saveRDS(net, "string_network.rds")

# --- Reload in a later session ---

# From TSV files (fileTypeInput = "stg" accepts file paths or data.frames)
levi(
    networkCoordinatesInput  = "string_nodes.tsv",
    networkInteractionsInput = "string_edges.tsv",
    fileTypeInput            = "stg",
    ...
)

# From the RDS object
net <- readRDS("string_network.rds")
levi(
    networkCoordinatesInput  = net$nodes,
    networkInteractionsInput = net$edges,
    fileTypeInput            = "stg",
    ...
)

Tip — persistent STRINGdb cache: pass a fixed input_directory to leviFromSTRING() so the raw STRING files are stored locally and reused across sessions without re-downloading:

net <- leviFromSTRING(
    genes           = my_genes,
    species         = 9606,
    input_directory = "~/.stringdb_cache"   # persists between R sessions
)

Layout guide

Species codes

Full pipeline example (RNA-seq + STRING)

library(airway)
library(DESeq2)
library(levi)

# 1. Differential expression
data(airway)
dds <- DESeqDataSet(airway, design = ~cell + dex)
dds <- DESeq(dds)
res <- results(dds, contrast = c("dex", "trt", "untrt"))

# 2. Build expression data.frame
expr_df <- leviFromDESeq2(res)

# 3. Retrieve STRING network for top DE genes
top_genes <- head(expr_df$GeneID[order(abs(expr_df$log2FoldChange),
                                        decreasing = TRUE)], 100)
set.seed(42)
net <- leviFromSTRING(top_genes, species = 9606, score_threshold = 400)

# 4. Visualise landscape
levi(
    expressionInput          = expr_df,
    networkCoordinatesInput  = net$nodes,
    networkInteractionsInput = net$edges,
    fileTypeInput            = "stg",
    geneSymbolInput          = "GeneID",
    readExpColumn            = readExpColumn("log2FoldChange-log2FoldChange"),
    signal_mode              = "logfc",
    n_perm                   = 500,
    perm_side                = "both"
)

Signal contracts and reproducibility

Signal interpretation: ratio preserves Test/(Test + Control), with no min-max rescaling; equal nonzero inputs give 0.5. expressionLog = TRUE back-transforms log2 inputs only in this mode. A single ratio column means abundance/(abundance + 1), not a comparison with a control. logfc accepts two log-scale columns or one already computed logFC, mapping zero to 0.5. zscore centres on the mean logFC of measured network support points, not on biological absence of change. Missing measurements are assigned 0.5 and listed in result$metadata. Gaussian smoothing mixes neighbouring signals, so these baseline statements apply to the input signals and to uniformly neutral networks.

Permutation inference is conditional on the fixed network and layout. Measured gene values are shuffled as pairs and edge midpoints are recalculated; missing positions stay fixed. In the default regional mode the maximum regional mass over both directions is the reference statistic, which controls the search across regions without a further adjustment. In cell mode both tails over occupied cells form one multiple-testing family per comparison, adjusted with p_adjust_method (“BY” by default); result$raw_pvalues retains the unadjusted values. perm_side selects the displayed side without changing either family. This is not a test of differential expression between biological replicates; increasing n_perm alone does not validate inferential use. Calibration across networks, layouts and missingness patterns still requires simulation studies.

Call set.seed() before permutation runs. result$metadata records the RNG state, network, coordinates, signal mode, grid settings and software versions. leviDiff() rejects incompatible metadata or grid coordinates. For DESeq2, edgeR and limma adapters, select the logFC column against itself with signal_mode = "logfc"; abundance annotations are not control measurements. For edgeR, select the contrast in glmLRT() or glmQLFTest() before calling leviFromEdgeR() on the resulting test object. KEGG uses the supplied universe, converting both selected genes and background to ENTREZID with the same OrgDb.