## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.align = "center", fig.width = 6, fig.height = 5, message = FALSE, warning = FALSE ) library(levi) ## ----install, eval=FALSE------------------------------------------------------ # if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # BiocManager::install("levi") # # # Optional — needed for specific features # BiocManager::install(c( # "DESeq2", "edgeR", "limma", # DE adapters # "SummarizedExperiment", # leviFromSE # "STRINGdb", # leviFromSTRING # "clusterProfiler", "org.Hs.eg.db", # leviEnrich # "airway" # example dataset # )) # install.packages(c("Seurat", "ggrepel", "patchwork", "plotly")) ## ----readexp_demo------------------------------------------------------------- # Single comparison readExpColumn("TumorCurrentSmoker-NormalNeverSmoker") # Two comparisons — levi() will return a list of two results readExpColumn( "TumorCurrentSmoker-NormalNeverSmoker", "TumorFormerSmoker-NormalFormerSmoker" ) ## ----signal_mode_demo, fig.height=4, fig.width=10----------------------------- logfc_net <- file.path(system.file(package="levi"), "extdata", "logfc_network.dat") logfc_expr <- file.path(system.file(package="levi"), "extdata", "logfc_expression.dat") base_call <- list(networkCoordinatesInput = logfc_net, expressionInput = logfc_expr, fileTypeInput = "dat", geneSymbolInput = "ID", readExpColumn = readExpColumn("Test-Control"), contrastValueInput = 50, resolutionValueInput = 20, zoomValueInput = 50, smoothValueInput = 5) res_ratio <- do.call(levi, c(base_call, list(signal_mode = "ratio"))) res_logfc <- do.call(levi, c(base_call, list(signal_mode = "logfc"))) res_zscore <- do.call(levi, c(base_call, list(signal_mode = "zscore"))) cat(sprintf( "Score range ratio: %.3f | logfc: %.3f | zscore: %.3f\n", diff(range(res_ratio$scores$LandscapeScore)), diff(range(res_logfc$scores$LandscapeScore)), diff(range(res_zscore$scores$LandscapeScore)) )) ## ----hub---------------------------------------------------------------------- hub_net <- file.path(system.file(package="levi"), "extdata", "hub_network.dat") hub_expr <- file.path(system.file(package="levi"), "extdata", "hub_expression.dat") res_hub <- levi( networkCoordinatesInput = hub_net, expressionInput = hub_expr, fileTypeInput = "dat", geneSymbolInput = "ID", readExpColumn = readExpColumn("Test-Control"), contrastValueInput = 50, resolutionValueInput = 20, zoomValueInput = 50, smoothValueInput = 5, contourLevi = TRUE ) ## ----hub_scores--------------------------------------------------------------- res_hub$scores[, c("Gene", "LandscapeScore", "Rank")] ## ----gradient----------------------------------------------------------------- grad_net <- file.path(system.file(package="levi"), "extdata", "gradient_network.dat") grad_expr <- file.path(system.file(package="levi"), "extdata", "gradient_expression.dat") res_grad <- levi( networkCoordinatesInput = grad_net, expressionInput = grad_expr, fileTypeInput = "dat", geneSymbolInput = "ID", readExpColumn = readExpColumn("Test-Control"), contrastValueInput = 50, resolutionValueInput = 20, zoomValueInput = 50, smoothValueInput = 5 ) scores_ord <- res_grad$scores[order(res_grad$scores$Rank), c("Gene", "LandscapeScore")] scores_ord ## ----bimodal------------------------------------------------------------------ bim_net <- file.path(system.file(package="levi"), "extdata", "bimodal_network.dat") bim_expr <- file.path(system.file(package="levi"), "extdata", "bimodal_expression.dat") res_bim <- levi( networkCoordinatesInput = bim_net, expressionInput = bim_expr, fileTypeInput = "dat", geneSymbolInput = "ID", readExpColumn = readExpColumn("Test-Control"), contrastValueInput = 50, resolutionValueInput = 20, zoomValueInput = 50, smoothValueInput = 5, contourLevi = TRUE ) ## ----bimodal_check------------------------------------------------------------ cat("A_HUB score:", round(res_bim$scores$LandscapeScore[ res_bim$scores$Gene == "A_HUB"], 3), "\n") cat("B_HUB score:", round(res_bim$scores$LandscapeScore[ res_bim$scores$Gene == "B_HUB"], 3), "\n") ## ----flat--------------------------------------------------------------------- flat_net <- file.path(system.file(package="levi"), "extdata", "flat_network.dat") flat_expr <- file.path(system.file(package="levi"), "extdata", "flat_expression.dat") res_flat <- levi( networkCoordinatesInput = flat_net, expressionInput = flat_expr, fileTypeInput = "dat", geneSymbolInput = "ID", readExpColumn = readExpColumn("Test-Control"), contrastValueInput = 50, resolutionValueInput = 20, zoomValueInput = 50, smoothValueInput = 5 ) cat(sprintf("Score SD = %.4f (expected < 0.15)\n", sd(res_flat$scores$LandscapeScore))) ## ----sparse------------------------------------------------------------------- sparse_net <- file.path(system.file(package="levi"), "extdata", "sparse_network.dat") sparse_expr <- file.path(system.file(package="levi"), "extdata", "sparse_expression.dat") res_sparse <- levi( networkCoordinatesInput = sparse_net, expressionInput = sparse_expr, fileTypeInput = "dat", geneSymbolInput = "ID", readExpColumn = readExpColumn("Test-Control"), contrastValueInput = 50, resolutionValueInput = 20, zoomValueInput = 50, smoothValueInput = 5 ) cat("Nodes in scores table:", nrow(res_sparse$scores), "(expected 15)\n") ## ----result_structure--------------------------------------------------------- # res_hub was computed above; the same six fields come back from every call. names(res_hub) str(res_hub$scores) # Gene, X, Y, LandscapeScore, Rank str(res_hub$peaks) # Type (peak/valley), NearestGene, ..., Score str(res_hub$landscape) # the plotted surface: Var1, Var2, z res_hub$pvalues # NULL here, because n_perm = 0 class(res_hub$plot) # the ggplot object ## ----scores_demo-------------------------------------------------------------- # Top 5 genes by landscape score head(res_hub$scores[, c("Gene", "LandscapeScore", "Rank")], 5) # Detected peaks and valleys if (!is.null(res_bim$peaks) && nrow(res_bim$peaks) > 0) res_bim$peaks[, c("Type", "NearestGene", "Score")] ## ----perm--------------------------------------------------------------------- set.seed(42) res_perm <- levi( networkCoordinatesInput = grad_net, expressionInput = grad_expr, fileTypeInput = "dat", geneSymbolInput = "ID", readExpColumn = readExpColumn("Test-Control"), contrastValueInput = 50, resolutionValueInput = 20, zoomValueInput = 50, smoothValueInput = 5, contourLevi = TRUE, n_perm = 50, perm_side = "both", sig_level = 0.05 ) ## ----perm-regions------------------------------------------------------------- res_perm$regions$summary[, c("Region", "Direction", "Cells", "Mass", "PSpatial", "Significant")] ## ----batch, fig.show="hold", out.width="50%"---------------------------------- 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") res_list <- levi( networkCoordinatesInput = hub_net, expressionInput = mc_expr, fileTypeInput = "dat", geneSymbolInput = "ID", readExpColumn = readExpColumn("Cond_A-Cond_B", "Cond_A-Cond_C"), contrastValueInput = 50, resolutionValueInput = 20, zoomValueInput = 50, smoothValueInput = 5 ) cat("Number of comparisons:", length(res_list), "\n") cat("Cond_A-Cond_B score range:", round(diff(range(res_list[[1]]$scores$LandscapeScore)), 3), "\n") cat("Cond_A-Cond_C score range:", round(diff(range(res_list[[2]]$scores$LandscapeScore)), 3), "\n") ## ----levigrid, fig.show="hold", out.width="50%"------------------------------- leviGrid(res_list, titles = c("Cond A vs B", "Cond A vs C"), ncol = 2) ## ----levidiff----------------------------------------------------------------- leviDiff( res_list[[1]], res_list[[2]], label_a = "Cond_A vs Cond_B", label_b = "Cond_A vs Cond_C" ) ## ----deseq2_adapter, eval=FALSE----------------------------------------------- # library(DESeq2) # # dds <- DESeqDataSetFromMatrix(counts, colData, design = ~condition) # dds <- DESeq(dds) # res_de <- results(dds, contrast = c("condition", "treated", "untreated")) # # expr_df <- leviFromDESeq2(res_de, gene_col = "GeneID") # # Columns: GeneID, baseMean (abundance annotation), log2FoldChange (logFC signal) # # levi( # expressionInput = expr_df, # # ... # readExpColumn = readExpColumn("log2FoldChange-log2FoldChange"), # signal_mode = "logfc" # ) ## ----edger_adapter, eval=FALSE------------------------------------------------ # library(edgeR) # # dge <- DGEList(counts = counts, group = group) # dge <- calcNormFactors(dge) # fit <- glmQLFit(dge, design) # res <- glmQLFTest(fit, coef = 2) # # expr_df <- leviFromEdgeR(res, gene_col = "GeneID") # # Columns: GeneID, logCPM (abundance annotation), logFC (logFC signal) ## ----limma_adapter, eval=FALSE------------------------------------------------ # library(limma) # # fit <- lmFit(eset, design) # fit2 <- contrasts.fit(fit, makeContrasts(TvsC = Tumor - Control, # levels = design)) # fit2 <- eBayes(fit2) # # expr_df <- leviFromLimma(fit2, coef = "TvsC", gene_col = "GeneID") # # Columns: GeneID, AveExpr (abundance annotation), logFC (logFC signal) ## ----seurat_adapter, eval=FALSE----------------------------------------------- # library(Seurat) # # markers <- FindMarkers(seurat_obj, # ident.1 = "CD4_T", # ident.2 = "B_cell", # min.pct = 0.25) # # expr_df <- leviFromSeurat(markers, gene_col = "GeneID") # # Columns: GeneID, Control (pct.2 annotation), Test (avg_log2FC signal) ## ----se_adapter, eval=FALSE--------------------------------------------------- # library(SummarizedExperiment) # # # Option A: condition labels from colData # expr_df <- leviFromSE( # se = my_se, # assay_name = "counts", # condition_col = "treatment", # test_level = "treated", # ctrl_level = "control", # gene_col = "GeneID" # ) # # # Option B: explicit sample names # expr_df <- leviFromSE( # se = my_se, # test_col = c("Sample1", "Sample3"), # ctrl_col = c("Sample2", "Sample4") # ) ## ----adapters_live------------------------------------------------------------ genes <- c("HUB", paste0("N", 1:8)) # DESeq2-shaped results table res_de <- data.frame( baseMean = rep(1000, 9), log2FoldChange = c(4.3, 4.1, 4.4, 4.2, 4.3, -5.3, -5.1, -5.4, -5.2), row.names = genes) expr_de <- leviFromDESeq2(res_de, gene_col = "GeneID") head(expr_de, 3) # edgeR-shaped and limma-shaped tables leviFromEdgeR(data.frame(logFC = 4.3, logCPM = 10.2, row.names = "HUB"), gene_col = "GeneID") leviFromLimma(data.frame(logFC = 4.3, AveExpr = 8.1, row.names = "HUB"), gene_col = "GeneID") # Seurat-shaped markers: pct.2 becomes Control, avg_log2FC becomes Test leviFromSeurat(data.frame(avg_log2FC = 2.5, pct.2 = 0.3, row.names = "HUB"), gene_col = "GeneID") # A SummarizedExperiment, aggregated by condition label counts <- matrix(c(200, 200, 200, 200, 200, 5, 5, 5, 5, 190, 210, 195, 205, 200, 6, 4, 5, 5, 10, 10, 10, 10, 10, 200, 200, 200, 200), nrow = 9, dimnames = list(genes, c("t1", "t2", "n1"))) se <- SummarizedExperiment::SummarizedExperiment( assays = list(counts = counts), colData = data.frame(condition = c("Tumor", "Tumor", "Normal"), row.names = colnames(counts))) expr_se <- leviFromSE(se, assay_name = "counts", condition_col = "condition", test_level = "Tumor", ctrl_level = "Normal", gene_col = "GeneID") head(expr_se, 3) ## ----adapters_to_levi, fig.height=5, fig.width=5------------------------------ hub_net <- system.file("extdata", "hub_network.dat", package = "levi") res_from_de <- levi( expressionInput = expr_de, networkCoordinatesInput = hub_net, fileTypeInput = "dat", geneSymbolInput = "GeneID", readExpColumn = readExpColumn("log2FoldChange-log2FoldChange"), resolutionValueInput = 40, signal_mode = "logfc") head(res_from_de$scores[, c("Gene", "LandscapeScore", "Rank")], 4) ## ----string_basic, eval=FALSE------------------------------------------------- # library(levi) # # BiocManager::install("STRINGdb") # # # MAPK pathway genes (human) # mapk_genes <- c("EGFR", "KRAS", "BRAF", "MAP2K1", "MAPK1", # "MAPK3", "RPS6KA1", "MYC", "JUN", "FOS") # # set.seed(42) # layout is stochastic; fix seed for reproducibility # net <- leviFromSTRING( # genes = mapk_genes, # species = 9606, # human # score_threshold = 400, # medium confidence # layout = "fr" # Fruchterman-Reingold # ) # # # net$nodes : data.frame — name, x, y # # net$edges : data.frame — V1, V2 # # net$graph : igraph object # # levi( # expressionInput = my_de_results, # networkCoordinatesInput = net$nodes, # networkInteractionsInput = net$edges, # fileTypeInput = "stg", # geneSymbolInput = "GeneID", # readExpColumn = readExpColumn("log2FoldChange-log2FoldChange"), # signal_mode = "logfc" # ) ## ----string_save, eval=FALSE-------------------------------------------------- # # Save to TSV — readable directly by levi() as file paths # 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) # # # Save full R object (preserves igraph + layout) # saveRDS(net, "string_network.rds") # # # Reload in future sessions # net2 <- readRDS("string_network.rds") # levi(networkCoordinatesInput = net2$nodes, # networkInteractionsInput = net2$edges, # fileTypeInput = "stg", ...) # # # Or use STRINGdb local cache to avoid re-downloading raw files # net3 <- leviFromSTRING(genes, input_directory = "~/.stringdb_cache") ## ----airway, eval=FALSE------------------------------------------------------- # 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")) # # 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) # # 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" # ) ## ----all_leukemia, eval=FALSE------------------------------------------------- # BiocManager::install(c("ALL", "limma", "STRINGdb")) # library(ALL); library(limma); library(levi) # # data(ALL) # 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") # # levi( # expressionInput = expr_df, # networkCoordinatesInput = net$nodes, # networkInteractionsInput = net$edges, # fileTypeInput = "stg", # geneSymbolInput = "GeneID", # readExpColumn = readExpColumn("logFC-logFC"), # signal_mode = "logfc", # logfc_k = 2 # ) ## ----pbmc, eval=FALSE--------------------------------------------------------- # BiocManager::install("TENxPBMCData") # install.packages("Seurat") # library(TENxPBMCData); library(Seurat); library(levi) # # pbmc_sce <- TENxPBMCData("pbmc3k") # pbmc <- as.Seurat(pbmc_sce) # pbmc <- NormalizeData(pbmc) # pbmc <- FindVariableFeatures(pbmc) # pbmc <- ScaleData(pbmc) # pbmc <- RunPCA(pbmc) # pbmc <- FindNeighbors(pbmc) # pbmc <- FindClusters(pbmc, resolution = 0.5) # # markers <- FindMarkers(pbmc, # ident.1 = "CD4 T cells", # ident.2 = "B cells", # min.pct = 0.25) # expr_df <- leviFromSeurat(markers, gene_col = "GeneID") # # 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 # ) ## ----enrich, eval=FALSE------------------------------------------------------- # BiocManager::install(c("clusterProfiler", "org.Hs.eg.db")) # # enrich_res <- leviEnrich( # result = res_hub, # top_n = 20, # top-scoring genes (peaks) # bottom_n = 20, # bottom-scoring genes (valleys) # organism = "hsa", # KEGG organism code # orgdb = "org.Hs.eg.db", # keytype = "SYMBOL", # pval_cutoff = 0.05, # types = c("GO_BP", "KEGG") # ) # # # enrich_res$over — enrichment on peak genes # # enrich_res$under — enrichment on valley genes # # # Dotplots are printed automatically # # Access results programmatically: # head(as.data.frame(enrich_res$over$GO_BP)) ## ----plot3d, eval=FALSE------------------------------------------------------- # install.packages("plotly") # # levi( # networkCoordinatesInput = hub_net, # expressionInput = hub_expr, # fileTypeInput = "dat", # geneSymbolInput = "ID", # readExpColumn = readExpColumn("Test-Control"), # contrastValueInput = 50, # resolutionValueInput = 20, # zoomValueInput = 50, # smoothValueInput = 5, # plot3d = TRUE # ) ## ----gui, eval=FALSE---------------------------------------------------------- # LEVIui(browser = FALSE) # open in RStudio Viewer # LEVIui(browser = TRUE) # open in system browser ## ----palettes, eval=FALSE----------------------------------------------------- # # Multicolor (default) # levi(..., setcolor = "default") # # # Two-color options # levi(..., setcolor = "purple_pink") # levi(..., setcolor = "green_blue") # levi(..., setcolor = "blue_yellow") # levi(..., setcolor = "pink_green") # levi(..., setcolor = "orange_purple") # levi(..., setcolor = "green_marine") ## ----session------------------------------------------------------------------ sessionInfo()