This file demonstrates the examples showcase static tasks spanning 10 omics domains and corresponding code responses generated by large language model.
Below are 10 representative bioinformatics tasks from our 1,002-task benchmark, covering genomics, transcriptomics, epigenomics, proteomics, metabolomics, single-cell omics, spatial omics, microbiomics, pharmacogenomics, and multi-omics integration.
The following R code demonstrates how to generate LLM responses for benchmark tasks using the llmhelper package:
# Install llmhelper package
if (!requireNamespace("devtools", quietly = TRUE)) {
install.packages("devtools")
}
devtools::install_github('SolvingLab/OmixBench', subdir = 'llmhelper')
# Load required libraries
library(llmhelper)
library(glue)
# Configure LLM client
response_client <- llm_provider(
model = 'o4-mini',
temperature = 0.2,
max_tokens = 20000,
timeout = 500,
stream = FALSE,
verbose = FALSE
)
# Define tasks
task_examples <- c(
"Identify tumor mutational burden (TMB) in a melanoma patient cohort by analyzing VCF files from whole exome sequencing, and generate a summary table of TMB values per patient.",
"Perform WGCNA on RNA-seq data from a breast cancer patient cohort to identify co-expression modules associated with tumor subtypes and generate network plots displaying module eigengene relationships.",
"Identify enriched transcription factor motifs in ATAC-seq peak regions derived from a breast cancer patient cohort and output a summary table comparing motif frequencies against a whole-genome background.",
"Identify hub proteins using a protein-protein interaction network constructed from mass spectrometry-based proteomics data in a breast cancer cohort, and generate a summary table of centrality metrics to prioritize potential therapeutic targets.",
"Construct a metabolic network based on GC-MS metabolomics profiles from a type 2 diabetes cohort to reveal key metabolic hubs, producing a corresponding network diagram and summary table of pathway connectivity.",
"Apply iClusterPlus to integrate RNA-seq expression data and DNA methylation profiles from a breast cancer patient cohort to identify tumor subtypes and produce a heatmap illustrating cluster assignments.",
"Identify cell differentiation trajectories and branch events in single-cell RNA-seq data from acute myeloid leukemia patients, and generate a detailed trajectory plot for visualizing pseudo-time progression.",
"Identify spatial ligand-receptor pairs using RNA-seq based spatial omics data from a tumor microenvironment, and generate a heat map summarizing intercellular interaction strengths for a subset of tumor subtypes.",
"Visualize co-occurrence networks of microbial taxa using 16S rRNA gene sequencing data from an inflammatory bowel disease (IBD) patient cohort, generating an interactive network plot that highlights both positive and negative correlations.",
"Generate dose–response curves for multiple targeted therapies using viability data from prostate cancer cell lines treated at varying concentrations; output should include individual plots with fitted curves and EC50 annotations."
)
# Generate responses for each task
for (task in task_examples) {
cat(paste0('\n### Task: ', task, '\n\n'))
response <- get_llm_response(
set_prompt(
system = "You are an AI assistant specialized in bioinformatics. Provide complete R code solutions.",
user = glue::glue("Please provide R code to solve the following bioinformatics task:\n\n{task}")
),
response_client,
max_retries = 1,
verbose = TRUE
)
cat('\n')
}
## Task: Identify tumor mutational burden (TMB) in a melanoma patient cohort by analyzing VCF files from whole exome sequencing, and generate a summary table of TMB values per patient.
## --- Sending request to LLM provider (o4-mini): ---
## Please provide R code to solve the following bioinformatics task:
##
## Identify tumor mutational burden (TMB) in a melanoma patient cohort by analyzing VCF files from whole exome sequencing, and generate a summary table of TMB values per patient.
##
## You must provide a text response.
## --- Receiving response from LLM provider: ---
## Here is a self-contained R script that
##
## 1. Reads all VCFs in a directory
## 2. Filters for high‐quality SNVs in coding (exonic) regions
## 3. Counts the mutations per sample
## 4. Calculates TMB as mutations per megabase of target exome
## 5. Writes out a summary table of TMB per patient
##
## ```r
## # Install required packages if not already installed
## if (!requireNamespace("BiocManager", quietly=TRUE))
## install.packages("BiocManager")
##
## pkgList <- c("VariantAnnotation", "GenomicFeatures", "GenomicRanges", "dplyr", "data.table")
## for (pkg in pkgList) {
## if (!requireNamespace(pkg, quietly=TRUE))
## BiocManager::install(pkg, ask=FALSE)
## }
##
## library(VariantAnnotation)
## library(GenomicFeatures)
## library(GenomicRanges)
## library(dplyr)
## library(data.table)
##
## # 1. Load a TxDb for human coding regions (adjust genome/build as needed)
## # Here we assume hg19/UCSC; for hg38 use TxDb.Hsapiens.UCSC.hg38.knownGene
## txdb <- makeTxDbFromUCSC(genome="hg19", tablename="knownGene")
##
## # Extract all exonic ranges
## exons_gr <- exons(txdb)
##
## # Calculate total exome size in megabases
## total_exome_bp <- sum(width(reduce(exons_gr)))
## total_exome_mb <- total_exome_bp / 1e6
## message("Total exome size (MB): ", round(total_exome_mb, 2))
##
## # 2. Point to your directory of VCFs
## vcf_dir <- "path/to/vcf_directory" # <--- set this to your folder
## vcf_files <- list.files(vcf_dir, pattern="\\.vcf(\\.gz)?$", full.names=TRUE)
##
## # Prepare a data.table to hold results
## results <- data.table(sample=character(), n_mutations=integer(), TMB=numeric())
##
## # 3. Loop through VCFs
## for (vcf_file in vcf_files) {
## # Derive sample name from file
## sample_name <- tools::file_path_sans_ext(basename(vcf_file))
##
## # Read VCF (assume single-sample VCFs)
## vcf <- readVcf(vcf_file, genome="hg19")
##
## # Keep only PASS variants
## vcf <- vcf[fixed(vcf)$FILTER == "PASS", ]
##
## # Keep only SNVs
## is_snv <- width(ref(vcf)) == 1 & width(unlist(alt(vcf))) == 1
## vcf <- vcf[is_snv, ]
##
## # Map VCF rowRanges to GRanges
## vr <- rowRanges(vcf)
##
## # Find overlaps with exons
## hits <- findOverlaps(vr, exons_gr)
## vr_exonic <- vr[unique(queryHits(hits))]
##
## # Count exonic SNVs
## n_mut <- length(vr_exonic)
##
## # Calculate TMB = mutations per MB
## tmb_val <- n_mut / total_exome_mb
##
## # Add to results
## results <- rbind(results, data.table(
## sample = sample_name,
## n_mutations = n_mut,
## TMB = round(tmb_val, 2)
## ))
##
## message(sprintf("Processed %s: %d mutations, TMB=%.2f",
## sample_name, n_mut, tmb_val))
## }
##
## # 4. Save summary table
## out_file <- "TMB_summary_per_sample.csv"
## fwrite(results, file=out_file)
##
## message("Summary written to: ", out_file)
## ```
##
## How it works:
##
## • Loads a UCSC TxDb to fetch all exonic (coding) ranges
## • Calculates the total exome footprint in megabases
## • Iterates over each VCF file, keeping only PASS SNVs
## • Intersects variants with exonic regions and counts them
## • Computes TMB = (number of exonic SNVs) / (exome size in MB)
## • Writes a CSV with columns: sample, n_mutations, TMB
## ......
## ......Note: The llmhelper package facilitates LLM integration in R environments. For detailed documentation, see:
-
Package source: https://github.com/SolvingLab/OmixBench/llmhelper
-
User manual: https://github.com/SolvingLab/OmixBench/llmhelper_1.0.0.pdf
Domain: Genomics
Task: Identify tumor mutational burden (TMB) in a melanoma patient cohort by analyzing VCF files from whole exome sequencing, and generate a summary table of TMB values per patient.
View Generated Code
# Install required packages if not already installed
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
pkgList <- c("VariantAnnotation", "GenomicFeatures", "GenomicRanges", "dplyr", "data.table")
for (pkg in pkgList) {
if (!requireNamespace(pkg, quietly=TRUE))
BiocManager::install(pkg, ask=FALSE)
}
library(VariantAnnotation)
library(GenomicFeatures)
library(GenomicRanges)
library(dplyr)
library(data.table)
# Load a TxDb for human coding regions
txdb <- makeTxDbFromUCSC(genome="hg19", tablename="knownGene")
exons_gr <- exons(txdb)
# Calculate total exome size in megabases
total_exome_bp <- sum(width(reduce(exons_gr)))
total_exome_mb <- total_exome_bp / 1e6
message("Total exome size (MB): ", round(total_exome_mb, 2))
# Point to your directory of VCFs
vcf_dir <- "path/to/vcf_directory"
vcf_files <- list.files(vcf_dir, pattern="\\.vcf(\\.gz)?$", full.names=TRUE)
results <- data.table(sample=character(), n_mutations=integer(), TMB=numeric())
# Loop through VCFs
for (vcf_file in vcf_files) {
sample_name <- tools::file_path_sans_ext(basename(vcf_file))
vcf <- readVcf(vcf_file, genome="hg19")
vcf <- vcf[fixed(vcf)$FILTER == "PASS", ]
is_snv <- width(ref(vcf)) == 1 & width(unlist(alt(vcf))) == 1
vcf <- vcf[is_snv, ]
vr <- rowRanges(vcf)
hits <- findOverlaps(vr, exons_gr)
vr_exonic <- vr[unique(queryHits(hits))]
n_mut <- length(vr_exonic)
tmb_val <- n_mut / total_exome_mb
results <- rbind(results, data.table(
sample = sample_name,
n_mutations = n_mut,
TMB = round(tmb_val, 2)
))
message(sprintf("Processed %s: %d mutations, TMB=%.2f",
sample_name, n_mut, tmb_val))
}
# Save summary table
out_file <- "TMB_summary_per_sample.csv"
fwrite(results, file=out_file)
message("Summary written to: ", out_file)Domain: Transcriptomics
Task: Perform WGCNA on RNA-seq data from a breast cancer patient cohort to identify co-expression modules associated with tumor subtypes and generate network plots displaying module eigengene relationships.
View Generated Code
# Load libraries
if (!require("WGCNA")) {
install.packages("WGCNA")
library(WGCNA)
}
options(stringsAsFactors = FALSE)
allowWGCNAThreads()
# Read in expression data and clinical traits
exprFile <- "breast_cancer_expr.csv"
traitFile <- "breast_cancer_traits.csv"
datExpr0 <- read.csv(exprFile, row.names = 1, check.names = FALSE)
datExpr <- as.data.frame(t(datExpr0))
datTraits <- read.csv(traitFile, row.names = 1)
commonSamples <- intersect(rownames(datExpr), rownames(datTraits))
datExpr <- datExpr[commonSamples, ]
datTraits <- datTraits[commonSamples, ]
# Data cleaning
gsg <- goodSamplesGenes(datExpr, verbose = 3)
if (!gsg$allOK) {
if (any(!gsg$goodGenes)) {
datExpr <- datExpr[, gsg$goodGenes]
}
if (any(!gsg$goodSamples)) {
datExpr <- datExpr[gsg$goodSamples, ]
datTraits <- datTraits[gsg$goodSamples, ]
}
}
# Sample clustering
sampleTree <- hclust(dist(datExpr), method = "average")
plot(sampleTree, main = "Sample clustering to detect outliers",
sub = "", xlab = "", cex.lab = 1.2, cex.axis = 1.2, cex.main = 1.4)
# Choose soft-thresholding power
powers <- c(1:20)
sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
par(mfrow = c(1,2))
cex1 <- 0.9
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3]) * sft$fitIndices[,2],
xlab="Soft Threshold (power)", ylab="Scale Free Topology Model Fit, signed R^2",
type="n", main = "Scale independence")
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3]) * sft$fitIndices[,2],
labels=powers, cex=cex1, col="red")
abline(h=0.80, col="red")
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)", ylab="Mean connectivity", type="n",
main = "Mean connectivity")
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1, col="red")
softPower <- 6
# Construct network and detect modules
adjacency <- adjacency(datExpr, power = softPower)
TOM <- TOMsimilarity(adjacency)
dissTOM <- 1 - TOM
geneTree <- hclust(as.dist(dissTOM), method = "average")
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04)
minModuleSize <- 30
dynamicMods <- cutreeDynamic(dendro = geneTree, distM = dissTOM,
deepSplit = 2, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize)
moduleColors <- labels2colors(dynamicMods)
plotDendroAndColors(geneTree, moduleColors,
"Dynamic Tree Cut", dendroLabels = FALSE,
hang = 0.03, addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors")
# Calculate module eigengenes
MEList <- moduleEigengenes(datExpr, colors = moduleColors)
MEs <- MEList$eigengenes
MEDiss <- 1 - cor(MEs)
METree <- hclust(as.dist(MEDiss), method = "average")
plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")
mergeCutHeight <- 0.25
merge <- mergeCloseModules(datExpr, moduleColors, cutHeight = mergeCutHeight, verbose = 3)
mergedColors <- merge$colors
mergedMEs <- merge$newMEs
moduleColors <- mergedColors
MEs <- mergedMEs
plotDendroAndColors(geneTree, cbind(dynamic = labels2colors(dynamicMods),
merged = moduleColors),
c("Dynamic Tree Cut", "Merged modules"),
hang = 0.03, addGuide = TRUE, guideHang = 0.05)
# Relate modules to external traits
MEs <- moduleEigengenes(datExpr, moduleColors)$eigengenes
moduleTraitCor <- cor(MEs, datTraits, use = "p")
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(datExpr))
textMatrix <- paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "")
par(mar = c(6, 8.5, 3, 3))
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = colnames(datTraits),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = blueWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
zlim = c(-1,1),
main = "Module–Trait Relationships")
# Plot eigengene network
plotEigengeneNetworks(MEs, "", marDendro = c(2,4,2,2),
marHeatmap = c(3,4,2,2), setLayout = FALSE)
# Export module to Cytoscape
moduleOfInterest <- "blue"
probes <- colnames(datExpr)
inModule <- (moduleColors == moduleOfInterest)
modProbes <- probes[inModule]
modTOM <- TOM[inModule, inModule]
edgeFile <- paste0("CytoscapeEdges_", moduleOfInterest, ".txt")
nodeFile <- paste0("CytoscapeNodes_", moduleOfInterest, ".txt")
exportNetworkToCytoscape(modTOM,
edgeFile = edgeFile,
nodeFile = nodeFile,
weighted = TRUE,
threshold = 0.02,
nodeNames = modProbes,
altNodeNames = modProbes,
nodeAttr = moduleColors[inModule])Domain: Epigenomics
Task: Identify enriched transcription factor motifs in ATAC-seq peak regions derived from a breast cancer patient cohort and output a summary table comparing motif frequencies against a whole-genome background.
View Generated Code
# Install and load required packages
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
pkgs <- c(
"BSgenome.Hsapiens.UCSC.hg38",
"GenomicRanges",
"rtracklayer",
"TFBSTools",
"JASPAR2020",
"motifmatchr",
"dplyr"
)
for (p in pkgs) {
if (!requireNamespace(p, quietly=TRUE))
BiocManager::install(p, ask=FALSE)
library(p, character.only=TRUE)
}
# Import ATAC-seq peaks
peaks_file <- "atac_peaks.bed"
peaks <- rtracklayer::import(peaks_file, format="BED")
# Define reference genome
genome <- BSgenome.Hsapiens.UCSC.hg38
# Fetch human TF PWMs from JASPAR2020
opts <- list(species="Homo sapiens", all_versions=FALSE)
motifs <- TFBSTools::getMatrixSet(JASPAR2020, opts)
message("Number of motifs loaded: ", length(motifs))
# Build length-matched random genomic background
set.seed(42)
peak_lengths <- width(peaks)
n_peaks <- length(peaks)
seqlens_all <- seqlengths(genome)[names(seqlengths(genome))]
chr_bg <- sample(
names(seqlens_all),
n_peaks,
replace = TRUE,
prob = seqlens_all / sum(seqlens_all)
)
bg_starts <- mapply(function(chr, w) {
max_start <- seqlens_all[chr] - w + 1
sample.int(max_start, 1)
}, chr_bg, peak_lengths)
bg_regions <- GenomicRanges::GRanges(
seqnames = chr_bg,
ranges = IRanges::IRanges(start = bg_starts, width = peak_lengths)
)
message("Background regions generated: ", length(bg_regions))
# Scan both peaks and background for motif presence
peak_matches <- motifmatchr::matchMotifs(
pwms = motifs,
subject = peaks,
genome = genome
)
bg_matches <- motifmatchr::matchMotifs(
pwms = motifs,
subject = bg_regions,
genome = genome
)
peak_mat <- assay(peak_matches, "motifMatches")
bg_mat <- assay(bg_matches, "motifMatches")
# Count how many regions contain each motif
peak_counts <- colSums(peak_mat)
bg_counts <- colSums(bg_mat)
N_peaks <- nrow(peak_mat)
N_bg <- nrow(bg_mat)
# Run Fisher's exact test
library(dplyr)
results <- data.frame(
motif_id = names(peak_counts),
motif_name = sapply(motifs[names(peak_counts)], function(x) x@name),
peak_with = as.integer(peak_counts),
peak_without = as.integer(N_peaks - peak_counts),
bg_with = as.integer(bg_counts),
bg_without = as.integer(N_bg - bg_counts),
stringsAsFactors = FALSE
) %>%
rowwise() %>%
mutate(
a = peak_with; b = peak_without;
c = bg_with; d = bg_without;
fish = list(fisher.test(matrix(c(a, b, c, d), nrow=2))),
odds_ratio = fish[[1]]$estimate,
p_value = fish[[1]]$p.value
) %>%
ungroup() %>%
mutate(
adj_p_value = p.adjust(p_value, method="BH"),
fold_enrichment = (peak_with / N_peaks) / (bg_with / N_bg)
) %>%
arrange(adj_p_value) %>%
select(
motif_id, motif_name,
peak_with, peak_without,
bg_with, bg_without,
odds_ratio, fold_enrichment,
p_value, adj_p_value
)
# Write summary table
out_file <- "motif_enrichment_summary.csv"
write.csv(results, file = out_file, row.names = FALSE)
message("Results saved to: ", out_file)Domain: Proteomics
Task: Identify hub proteins using a protein-protein interaction network constructed from mass spectrometry-based proteomics data in a breast cancer cohort, and generate a summary table of centrality metrics to prioritize potential therapeutic targets.
View Generated Code
library(readr)
library(dplyr)
library(igraph)
library(STRINGdb)
# Read in proteomics data
proteo_df <- read_csv("breast_cancer_proteomics.csv")
# Map Uniprot IDs to STRING internal IDs
string_db <- STRINGdb$new(
version="11.5",
species=9606,
score_threshold=700,
input_directory=""
)
mapped <- string_db$map(
proteo_df,
"Protein_ID",
removeUnmappedRows = TRUE
)
# Retrieve PPI edges
ppi_edges <- string_db$get_interactions(mapped$STRING_id)
# Build undirected igraph object
g <- graph_from_data_frame(
d = ppi_edges[, c("from","to","combined_score")],
directed = FALSE
)
E(g)$weight <- ppi_edges$combined_score
# Compute centrality metrics
centrality_df <- tibble(
STRING_id = V(g)$name,
degree = degree(g),
betweenness = betweenness(g, weights=1/E(g)$weight, normalized=TRUE),
closeness = closeness(g, weights=1/E(g)$weight, normalized=TRUE),
eigenvector = eigen_centrality(g, weights=E(g)$weight)$vector
)
# Merge back original Protein IDs
summary_df <- centrality_df %>%
left_join(
mapped %>% select(STRING_id, Protein_ID, gene_name),
by = "STRING_id"
) %>%
select(Protein_ID, gene_name, everything(), -STRING_id) %>%
arrange(desc(degree))
# Write output
write_csv(summary_df, "hub_protein_centrality_summary.csv")
message("Top 10 hub proteins by degree:")
print(head(summary_df, 10))Domain: Metabolomics
Task: Construct a metabolic network based on GC-MS metabolomics profiles from a type 2 diabetes cohort to reveal key metabolic hubs, producing a corresponding network diagram and summary table of pathway connectivity.
View Generated Code
required.pkgs <- c("tidyverse","Hmisc","igraph","ggraph","KEGGREST")
for (pkg in required.pkgs) {
if (!requireNamespace(pkg, quietly=TRUE)) install.packages(pkg)
}
library(tidyverse)
library(Hmisc)
library(igraph)
library(ggraph)
library(KEGGREST)
# Load data
met_file <- "metabolomics_data.csv"
meta_file <- "sample_metadata.csv"
map_file <- "metabolite_map.csv"
met.df <- read_csv(met_file)
meta.df <- read_csv(meta_file)
map.df <- read_csv(map_file)
met.mat <- met.df %>%
column_to_rownames("SampleID") %>%
as.matrix()
# Filter T2D samples
td2_samples <- meta.df %>%
filter(Group == "T2D") %>%
pull(SampleID)
met.t2d <- met.mat[td2_samples, , drop=FALSE]
# Compute Spearman correlations
rc <- Hmisc::rcorr(met.t2d, type="spearman")
cor.mat <- rc$r
p.mat <- rc$P
# Threshold to define edges
rho.thresh <- 0.6
p.thresh <- 0.05
edge_list <- which(abs(cor.mat) > rho.thresh & p.mat < p.thresh, arr.ind=TRUE)
edge_df <- tibble(
from = rownames(cor.mat)[edge_list[,1]],
to = colnames(cor.mat)[edge_list[,2]],
weight = cor.mat[edge_list]
) %>%
filter(as.character(from) < as.character(to))
# Build igraph object
g <- graph_from_data_frame(edge_df, directed=FALSE, vertices=map.df$Metabolite)
V(g)$degree <- igraph::degree(g)
V(g)$betweenness <- igraph::betweenness(g)
# Annotate with KEGG pathways
get_paths <- function(kegg_id) {
pw <- try(keggLink("pathway", paste0("cpd:", kegg_id)), silent=TRUE)
if (inherits(pw, "try-error") || length(pw)==0) return(NA_character_)
pw.ids <- sub("path:", "", pw)
pw.names <- sapply(pw.ids, function(x) {
info <- keggGet(paste0("path:", x))
if (!is.null(info[[1]]$NAME)) info[[1]]$NAME else NA_character_
})
paste0(pw.ids, ":", pw.names, collapse="|")
}
vertex_df <- as_tibble(as_data_frame(g, what="vertices")) %>%
rename(Metabolite = name) %>%
left_join(map.df, by="Metabolite")
path.lookup <- vertex_df %>%
distinct(KEGG_ID) %>%
mutate(Pathways = map_chr(KEGG_ID, ~ ifelse(is.na(.x), NA, get_paths(.x))))
vertex_df <- vertex_df %>%
left_join(path.lookup, by="KEGG_ID")
for (i in seq_len(nrow(vertex_df))) {
vname <- vertex_df$Metabolite[i]
V(g)[vname]$KEGG_ID <- vertex_df$KEGG_ID[i]
V(g)[vname]$Pathways <- vertex_df$Pathways[i]
V(g)[vname]$degree <- vertex_df$degree[i]
V(g)[vname]$betweenness <- vertex_df$betweenness[i]
}
# Pathway connectivity summary
path_edge_df <- edge_df %>%
left_join(vertex_df %>% select(Metabolite, Pathways), by=c("from"="Metabolite")) %>%
rename(pw_from = Pathways) %>%
left_join(vertex_df %>% select(Metabolite, Pathways), by=c("to"="Metabolite")) %>%
rename(pw_to = Pathways) %>%
mutate(pw_from = str_split(pw_from, "\\|"),
pw_to = str_split(pw_to, "\\|")) %>%
unnest(pw_from, pw_to) %>%
filter(!is.na(pw_from), !is.na(pw_to))
path_stats <- vertex_df %>%
filter(!is.na(Pathways)) %>%
separate_rows(Pathways, sep="\\|") %>%
group_by(Pathways) %>%
summarise(
nMetabolites = n(),
avg_degree = mean(degree),
avg_betweenness = mean(betweenness)
) %>%
arrange(desc(nMetabolites))
inter_pw <- path_edge_df %>%
distinct(from,to,pw_from,pw_to) %>%
filter(pw_from != pw_to) %>%
group_by(pw_from, pw_to) %>%
tally(name="edge_count") %>%
ungroup() %>%
rename(pathway1 = pw_from, pathway2 = pw_to)
print(path_stats)
print(inter_pw)
# Plot network
vertex_plot <- vertex_df %>%
mutate(primary_pw = str_split(Pathways, "\\|") %>% map_chr(~ .x[1])) %>%
replace_na(list(primary_pw="unmapped"))
V(g)$primary_pw <- vertex_plot$primary_pw[match(V(g)$name, vertex_plot$Metabolite)]
set.seed(123)
ggraph(g, layout="fr") +
geom_edge_link(aes(edge_alpha=abs(weight), edge_width=abs(weight)),
edge_colour="grey50") +
geom_node_point(aes(size=degree, color=primary_pw)) +
geom_node_text(aes(label=name), repel=TRUE, size=2) +
scale_size_continuous(range=c(2,8)) +
theme_void() +
labs(title="Metabolite Correlation Network in T2D",
subtitle="Edges: |Spearman ρ| > 0.6, p < 0.05",
color="Primary KEGG Pathway",
size="Degree")
write_csv(path_stats, "pathway_stats.csv")
write_csv(inter_pw, "inter_pathway_edges.csv")Domain: Multi-omics Integration
Task: Apply iClusterPlus to integrate RNA-seq expression data and DNA methylation profiles from a breast cancer patient cohort to identify tumor subtypes and produce a heatmap illustrating cluster assignments.
View Generated Code
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!requireNamespace("iClusterPlus", quietly = TRUE))
BiocManager::install("iClusterPlus")
install.packages(c("matrixStats","pheatmap"))
library(iClusterPlus)
library(matrixStats)
library(pheatmap)
# Load data
expr_counts <- read.csv("rna_counts.csv", row.names = 1, check.names = FALSE)
meth_beta <- read.csv("methylation_beta.csv", row.names = 1, check.names = FALSE)
common_samples <- intersect(colnames(expr_counts), colnames(meth_beta))
expr_counts <- expr_counts[, common_samples]
meth_beta <- meth_beta[, common_samples]
# Pre-processing
library(edgeR)
dge <- DGEList(counts = expr_counts)
dge <- calcNormFactors(dge)
log2cpm <- cpm(dge, log = TRUE, prior.count=1)
Mvals <- log2( meth_beta / (1 - meth_beta + 1e-8) + 1e-8 )
# Select top variable features
topN <- 500
gene_vars <- rowVars(log2cpm)
meth_vars <- rowVars(Mvals)
top_genes <- names(sort(gene_vars, decreasing=TRUE))[1:topN]
top_cpgs <- names(sort(meth_vars, decreasing=TRUE))[1:topN]
expr_sub <- log2cpm[top_genes, ]
meth_sub <- Mvals[top_cpgs, ]
expr_scaled <- t(scale(t(expr_sub)))
meth_scaled <- t(scale(t(meth_sub)))
data_list <- list(
RNAseq = expr_scaled,
Meth = meth_scaled
)
# Tune iClusterPlus
set.seed(123)
tune.res <- tune.iClusterPlus(
dt1 = data_list[[1]],
dt2 = data_list[[2]],
type = c("gaussian","gaussian"),
K = 3,
n.lambda = 5,
scale = TRUE,
maxiter = 20,
eps = 1e-3,
trace = FALSE
)
cv.errors <- tune.res$cv.error
mean.err <- rowMeans(cv.errors)
best.idx <- which.min(mean.err)
best.lambda.list <- tune.res$lambda[[best.idx]]
# Run final iClusterPlus
iclust.res <- iClusterPlus(
dt1 = data_list[[1]],
dt2 = data_list[[2]],
type = c("gaussian","gaussian"),
K = 3,
lambda = best.lambda.list,
scale = TRUE,
maxiter = 50,
eps = 1e-3,
trace = TRUE
)
clusters <- iclust.res$clusters
table(clusters)
# Build combined matrix for heatmap
combined_mat <- rbind(expr_scaled, meth_scaled)
rownames(combined_mat)[1:topN] <- paste0("G_", top_genes)
rownames(combined_mat)[(topN+1):(2*topN)] <- paste0("C_", top_cpgs)
annotation_col <- data.frame(
iCluster = factor(clusters)
)
rownames(annotation_col) <- colnames(combined_mat)
ann_colors <- list(
iCluster = setNames(
RColorBrewer::brewer.pal(3, "Set1"),
sort(unique(clusters))
)
)
pheatmap(
combined_mat,
show_rownames = FALSE,
show_colnames = FALSE,
annotation_col = annotation_col,
annotation_colors = ann_colors,
scale = "row",
clustering_method = "ward.D2",
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
main = "iClusterPlus integrated\nRNAseq + Methylation"
)Domain: Single-cell Omics
Task: Identify cell differentiation trajectories and branch events in single-cell RNA-seq data from acute myeloid leukemia patients, and generate a detailed trajectory plot for visualizing pseudo-time progression.
View Generated Code
if (!requireNamespace("Seurat", quietly=TRUE)) install.packages("Seurat")
if (!requireNamespace("monocle3", quietly=TRUE)) {
remotes::install_github("cole-trapnell-lab/monocle3")
}
if (!requireNamespace("patchwork", quietly=TRUE)) install.packages("patchwork")
if (!requireNamespace("ggplot2", quietly=TRUE)) install.packages("ggplot2")
library(Seurat)
library(monocle3)
library(ggplot2)
library(patchwork)
# Load 10× data
data_dir <- "path/to/your/10x_data"
counts <- Read10X(data.dir = data_dir)
# Build Seurat object
seurat_obj <- CreateSeuratObject(counts = counts,
project = "AML_scRNA",
min.cells = 3,
min.features = 200)
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")
seurat_obj <- subset(seurat_obj,
subset = nFeature_RNA > 200 &
nFeature_RNA < 2500 &
percent.mt < 5)
# Standard workflow
seurat_obj <- NormalizeData(seurat_obj, normalization.method = "LogNormalize",
scale.factor = 10000)
seurat_obj <- FindVariableFeatures(seurat_obj, selection.method = "vst",
nfeatures = 2000)
seurat_obj <- ScaleData(seurat_obj, features = rownames(seurat_obj))
seurat_obj <- RunPCA(seurat_obj, features = VariableFeatures(seurat_obj))
ElbowPlot(seurat_obj)
seurat_obj <- FindNeighbors(seurat_obj, dims = 1:20)
seurat_obj <- FindClusters(seurat_obj, resolution = 0.5)
seurat_obj <- RunUMAP(seurat_obj, dims = 1:20)
DimPlot(seurat_obj, reduction = "umap", label = TRUE) + NoLegend()
# Convert to Monocle3
cds <- as.cell_data_set(seurat_obj)
cds <- preprocess_cds(cds, num_dim = 50)
cds <- reduce_dimension(cds, reduction_method = "UMAP")
cds <- cluster_cells(cds, reduction_method = "UMAP")
# Learn trajectory graph
cds <- learn_graph(cds)
# Order cells in pseudotime
root_cells <- colnames(cds)[colData(cds)$seurat_clusters == "0"]
cds <- order_cells(cds, root_cells = root_cells)
# Plot trajectory
p1 <- plot_cells(cds,
color_cells_by = "pseudotime",
label_cell_groups = FALSE,
label_leaves = TRUE,
label_branch_points = TRUE,
graph_label_size = 4) +
ggtitle("Monocle3 Trajectory: Pseudotime Progression")
p2 <- plot_cells(cds,
color_cells_by = "seurat_clusters",
label_cell_groups = TRUE,
label_leaves = FALSE,
label_branch_points = FALSE,
graph_label_size = 4) +
ggtitle("Monocle3 Trajectory: Cell Clusters")
(p1 | p2) + plot_annotation(tag_levels = "A")
# Identify branch-associated genes
pr_graph_test_res <- graph_test(cds, neighbor_graph = "principal_graph", cores = 4)
sig_genes <- rownames(subset(pr_graph_test_res, q_value < 0.05))
length(sig_genes)
top_genes <- head(sig_genes[order(pr_graph_test_res[sig_genes, "morans_I"],
decreasing = TRUE)], 6)
plot_genes_in_pseudotime(cds[top_genes, ],
color_cells_by = "seurat_clusters") +
ggtitle("Top Branch-associated Genes Over Pseudotime")Domain: Spatial Omics
Task: Identify spatial ligand-receptor pairs using RNA-seq based spatial omics data from a tumor microenvironment, and generate a heat map summarizing intercellular interaction strengths for a subset of tumor subtypes.
View Generated Code
library(Seurat)
library(CellChat)
library(patchwork)
library(pheatmap)
library(dplyr)
# Load spatial data
data_dir <- "/path/to/your/spatial/data/"
se.sp <- Load10X_Spatial(
data.dir = data_dir,
filename = "filtered_feature_bc_matrix.h5"
)
# Preprocessing
se.sp <- SCTransform(se.sp, assay = "Spatial", verbose = FALSE)
se.sp <- RunPCA(se.sp, assay = "SCT", verbose = FALSE)
se.sp <- FindNeighbors(se.sp, dims = 1:30)
se.sp <- FindClusters(se.sp, resolution = 0.5)
se.sp <- RunUMAP(se.sp, dims = 1:30)
DimPlot(se.sp, group.by = "seurat_clusters", label = TRUE) + ggtitle("UMAP clusters")
# Define subtypes
subtypes <- c("SubtypeA", "SubtypeB", "SubtypeC")
table(se.sp$TumorSubtype)
se.list <- lapply(subtypes, function(sbt) {
subset(se.sp, subset = TumorSubtype == sbt)
})
names(se.list) <- subtypes
# Run CellChat
run_cellchat <- function(se_obj, group.by = "seurat_clusters") {
cc <- createCellChat(object = se_obj, group.by = group.by)
cc@DB <- CellChatDB.human
cc <- subsetData(cc)
cc <- identifyOverExpressedGenes(cc)
cc <- identifyOverExpressedInteractions(cc)
cc <- computeCommunProb(cc)
cc <- filterCommunication(cc, min.cells = 10)
cc <- computeCommunProbPathway(cc)
cc <- aggregateNet(cc)
return(cc)
}
cellchat.list <- lapply(se.list, run_cellchat)
# Extract interaction matrices
interaction.mats <- lapply(cellchat.list, function(cc) {
m <- cc@net$weight
m[is.na(m)] <- 0
return(as.matrix(m))
})
common.ct <- Reduce(intersect, lapply(interaction.mats, rownames))
interaction.mats <- lapply(interaction.mats, function(m) {
m2 <- m[common.ct, common.ct]
return(m2)
})
# Build summary table
pair.df <- do.call(rbind, lapply(names(interaction.mats), function(sbt) {
m <- interaction.mats[[sbt]]
v <- as.vector(m)
names(v) <- as.vector(outer(rownames(m), colnames(m), paste, sep = "_to_"))
data.frame(Subtype = sbt, Pair = names(v), Strength = v, row.names = NULL)
}))
wide.df <- pivot_wider(pair.df, names_from = Subtype, values_from = Strength,
values_fill = 0)
mat <- as.matrix(wide.df[,-1])
rownames(mat) <- wide.df$Pair
# Select top variable interactions
vars <- apply(mat, 1, var)
topN <- 50
sel <- names(sort(vars, decreasing = TRUE))[1:topN]
mat.top <- mat[sel, ]
# Plot heatmap
pheatmap(
mat.top,
scale = "row",
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "complete",
color = colorRampPalette(c("navy","white","firebrick3"))(50),
show_rownames = TRUE,
show_colnames = TRUE,
main = "Intercellular Interaction Strengths\n(top 50 most variable pairs)"
)Domain: Microbiomics
Task: Visualize co-occurrence networks of microbial taxa using 16S rRNA gene sequencing data from an inflammatory bowel disease (IBD) patient cohort, generating an interactive network plot that highlights both positive and negative correlations.
View Generated Code
library(phyloseq)
library(SpiecEasi)
library(igraph)
library(visNetwork)
library(dplyr)
# Assume phyloseq object 'ps' is loaded
# ps <- readRDS("your_phyloseq_object.rds")
# Pre-processing
ps_rel <- transform_sample_counts(ps, function(x) x / sum(x))
ps_filt <- filter_taxa(ps_rel,
function(x) sum(x > 1e-3) > 0.2 * length(x),
prune = TRUE)
# Network inference with SPIEC-EASI
se.mb <- spiec.easi(ps_filt,
method = "mb",
lambda.min.ratio = 1e-2,
nlambda = 20,
sel.crit = "stars",
pulsar.params = list(rep.num = 20))
beta_mat <- symBeta(getOptBeta(se.mb), mode = "maxabs")
# Build edge list
ij <- which(beta_mat != 0, arr.ind = TRUE)
edges <- data.frame(
from = rownames(beta_mat)[ij[,1]],
to = colnames(beta_mat)[ij[,2]],
weight = beta_mat[ij],
stringsAsFactors = FALSE
)
edges <- edges %>%
filter(from < to) %>%
mutate(
correlation = ifelse(weight > 0, "positive", "negative"),
color = ifelse(weight > 0, "red", "blue"),
width = 2 + 8 * abs(weight)
)
# Build node list
tax_table_df <- as.data.frame(tax_table(ps_filt))
otu_mat <- as.matrix(otu_table(ps_filt))
mean_abund <- rowMeans(otu_mat)
nodes <- data.frame(
id = names(mean_abund),
label = ifelse(is.na(tax_table_df$Genus),
rownames(tax_table_df),
tax_table_df$Genus),
title = paste0(
"<p><b>Taxon:</b> ", rownames(tax_table_df), "<br>",
"<b>Mean rel.abund:</b> ",
sprintf("%.3f", mean_abund), "</p>"
),
value = mean_abund,
stringsAsFactors = FALSE
)
# Interactive visualization
visNetwork(nodes, edges, width = "100%", height = "800px") %>%
visNodes(
shape = "dot",
color = list(
background = "#97C2FC",
border = "#2B7CE9",
highlight = "#FFD700"
)
) %>%
visEdges(
color = list(
highlight = "orange"
),
smooth = FALSE
) %>%
visOptions(
highlightNearest = list(enabled = TRUE, degree = 1, hover = TRUE),
nodesIdSelection = TRUE
) %>%
visLayout(randomSeed = 123) %>%
visPhysics(
solver = "forceAtlas2Based",
forceAtlas2Based = list(
gravitationalConstant = -50
),
stabilization = TRUE
)Domain: Pharmacogenomics
Task: Generate dose–response curves for multiple targeted therapies using viability data from prostate cancer cell lines treated at varying concentrations; output should include individual plots with fitted curves and EC50 annotations.
View Generated Code
if (!requireNamespace("drc", quietly=TRUE)) install.packages("drc")
if (!requireNamespace("tidyverse", quietly=TRUE)) install.packages("tidyverse")
library(drc)
library(tidyverse)
library(purrr)
# Read viability data
viability_data <- read_csv("viability_data.csv")
glimpse(viability_data)
# Nest, fit models, compute EC50 and generate plots
results_tbl <-
viability_data %>%
group_by(CellLine, Drug) %>%
nest() %>%
mutate(
drm_model = map(data, ~ drm(
Viability ~ Concentration,
data = .x,
fct = LL.4(names = c("Slope", "Lower", "Upper", "EC50"))
)),
EC50 = map_dbl(drm_model, ~ ED(.x, 50, interval="delta")[1, "Estimate"]),
newdata = map(data, ~ {
rng <- range(.x$Concentration)
tibble(Concentration = 10^seq(log10(rng[1]), log10(rng[2]), length.out = 100))
}),
preds = map2(drm_model, newdata, ~ {
tibble(
Concentration = .y$Concentration,
Fit = predict(.x, newdata = .y)
)
}),
plot = pmap(
list(data, preds, EC50, CellLine, Drug),
function(df_raw, df_pred, ec50_val, cell, drug) {
ggplot() +
geom_point(data = df_raw,
aes(x = Concentration, y = Viability),
color = "blue", size = 2) +
geom_line(data = df_pred,
aes(x = Concentration, y = Fit),
color = "red", size = 1) +
scale_x_log10() +
labs(
title = paste(cell, "–", drug),
x = "Drug Concentration (µM, log scale)",
y = "Cell Viability (%)"
) +
annotate("text",
x = min(df_raw$Concentration)*1.2,
y = max(df_raw$Viability),
label= paste0("EC50 = ", signif(ec50_val, 3)),
hjust= 0, size = 4, color = "darkred") +
theme_minimal()
}
)
) %>%
ungroup()
# Save plots
output_dir <- "dose_response_plots"
if (!dir.exists(output_dir)) dir.create(output_dir)
results_tbl %>%
mutate(
filename = paste0(output_dir, "/",
CellLine, "_", Drug, "_dose_response.png")
) %>%
pull(plot, filename) %>%
imap(~ ggsave(filename = .y, plot = .x,
width = 6, height = 4, dpi = 300))If you use this benchmark in your research, please cite:
@article{liu2025systematic,
title={Systematic Evaluation and Strategic Optimization of Large Language Models for Multi-omics Analysis},
author={Liu, Zaoqu and Wu, Yushuai and Yang, Jingkuan and others},
journal={In preparation},
year={2025}
}For questions or collaboration inquiries, please contact:
- Zaoqu Liu: liuzaoqu@163.com