diff --git a/README.md b/README.md index 22c89ed8..56e3ba5e 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ RamplotR is an open-source R Shiny application that brings backbone geometry, a | --- | --- | | **Explore a structure** | Interactive φ/ψ plots with several reference-density datasets, native RamplotR density regions and a parallel six-class Rama8000 standard validation that matches the current cctbx/Phenix categories. | | **Inspect residues in context** | Synchronized Ramachandran plot, searchable residue table, all-chain sequence navigator and NGL 3D viewer. The issue queue distinguishes Rama8000 outliers from native RamplotR `Not allowed` regions and missing angles. | -| **Work with predicted models** | Import AlphaFold DB models by UniProt accession or upload AlphaFold 2/3, ColabFold and ESMFold structures. Examine pLDDT, and view a linked PAE heatmap when compatible confidence data are available. | +| **Work with predicted models** | Import AlphaFold DB models or upload AlphaFold 2/3, ColabFold and ESMFold structures. Examine pLDDT/PAE and analyse multiple AF2/ColabFold or ESMFold seeds as a prediction ensemble with residue-level φ/ψ, Rama8000 and confidence agreement. | | **Examine structural geometry** | Explore peptide ω, side-chain χ1 and descriptive Cβ measurements. Optionally attach the matching deposited structure's official wwPDB validation report for independent rotamer, clash and geometry annotations. | | **Inspect experimental evidence** | Overlay a local CCP4/MRC cryo-EM map in the 3D viewer as a qualitative aid, without uploading the map to a separate service. | | **Compare models** | Sequence-align chains from two structures; use the Conformational Change Explorer to navigate residue-level wrapped φ/ψ displacement, then inspect paired residues in linked 2D/3D views alongside native and Rama8000 category changes. | @@ -105,6 +105,7 @@ The default RamplotR teal contour palette provides consistent, recognisable publ - [AlphaFold, ColabFold and ESMFold confidence analysis](docs/prediction-confidence.md) - [Geometry, official wwPDB evidence, cryo-EM overlays, ensembles and batch mode](docs/structural-verification.md) - [Rama8000 standard validation and direct wwPDB comparison](docs/rama8000-validation.md) +- [Prediction ensemble analysis](docs/prediction-ensembles.md) - [Independent wwPDB angle-validation protocol](docs/wwpdb-validation.md) · [Results](docs/validation-results.md) - [Performance and large-structure benchmarks](docs/benchmark-results.md) · [Scaling results](docs/scaling-results.md) - [Shinylive browser deployment and one.com caching](docs/shinylive-deployment.md) diff --git a/docs/inspection-user-guide.md b/docs/inspection-user-guide.md index 546ff22c..9e381a0f 100644 --- a/docs/inspection-user-guide.md +++ b/docs/inspection-user-guide.md @@ -83,6 +83,29 @@ comparison can exceed the alignment size limit; select shorter chains instead. Identical Ramachandran coordinates do not imply identical Cartesian structure, and an angular difference alone is not evidence of a clinically meaningful change. Distinct models from the same structure are related observations, not independent experiments. +## Prediction ensembles + +When the loaded structure is a prediction, the **Summary** tab exposes a +prediction-ensemble workflow even if the current coordinate file contains only +one model. Upload additional AF2/ColabFold, ESMFold or other models whose +B-factor field contains pLDDT. The currently loaded compatible model can be +included as one ensemble member. + +RamplotR reports circular phi/psi spread, residue coverage, Rama8000 agreement +and pLDDT spread across the uploaded models. A compact **Prediction variability +map** ranks residues by the larger circular SD of phi or psi. Clicking a map +cell or ensemble-table row selects the corresponding residue in the existing +inspector and linked 2D/3D views. + +The variability map is a navigation tool. Model-to-model disagreement is +described as prediction uncertainty or heterogeneity, not molecular dynamics. +Duplicate coordinate files are rejected, and the model export records labels, +declared source and coordinate MD5 hashes. + +AlphaFold 3 ensemble confidence is intentionally not inferred from B-factors in +this first version because each AF3 model needs its matching atom/token +confidence sidecar. + ## Publication exports and reproducibility The Summary tab includes an export disclosure. **Vector SVG** produces an editable figure; **High-resolution PNG** is suitable for manuscripts; **HTML report** includes the plot, summary counts, selected residue data and analysis settings. The residue CSV exports the current table filters; the comparison CSV exports the aligned comparison. diff --git a/docs/prediction-ensembles.md b/docs/prediction-ensembles.md new file mode 100644 index 00000000..8ab5e684 --- /dev/null +++ b/docs/prediction-ensembles.md @@ -0,0 +1,118 @@ +# Prediction ensemble analysis + +RamplotR can analyse multiple independently generated prediction models as a +single **prediction ensemble**. The purpose is to expose where prediction +seeds/models agree or disagree in local backbone geometry and confidence. + +This is deliberately not described as molecular dynamics. Variation between +predicted models can reflect model uncertainty, alternative plausible +solutions, stochastic sampling, different seeds or pipeline settings. It is +not experimental evidence that a protein moves between those conformations. + +## Supported inputs + +The first implementation accepts one coordinate file per model for: + +- AlphaFold 2 / ColabFold; +- ESMFold; +- other prediction models that store pLDDT in the B-factor field. + +The currently loaded compatible prediction can optionally be included as one +ensemble member. + +AlphaFold 3 is intentionally excluded from automatic ensemble confidence +analysis for now. AF3 confidence is atom/token based and needs its matching +confidence JSON sidecar; RamplotR does not silently reinterpret an AF3 model as +AF2. + +Each uploaded file must contain exactly one structural model. Up to 30 models +are analysed per run. + +## Residue matching + +Models are matched by: + +- chain; +- residue number; +- insertion code; +- residue identity. + +Missing residues remain missing and never receive fabricated coordinates, +angles or confidence scores. The Summary view reports how many residues are +present in every analysed model. + +Duplicate coordinate files are rejected using MD5 hashes so repeated copies of +one model cannot artificially inflate ensemble agreement. + +## Per-residue statistics + +For every matched residue RamplotR reports: + +- circular mean phi and psi; +- circular standard deviation of phi and psi; +- number of models contributing each angle; +- native RamplotR density-region agreement; +- Rama8000 Favored/Allowed/Outlier mode and agreement; +- mean, standard deviation, minimum and maximum pLDDT. + +Circular statistics are essential because +179 degrees and -179 degrees are +neighbours rather than opposite conformations. + +The prediction variability map uses the larger of phi-SD and psi-SD only as a +navigation aid: + +| Maximum circular SD | Display | +| --- | --- | +| <5 degrees | stable | +| 5-15 degrees | moderate | +| 15-30 degrees | variable | +| >=30 degrees | high | + +These bands are interface thresholds, not statistical significance cutoffs. + +A separate outline marks residues whose Rama8000 category differs between +models. + +## Model-level provenance + +The model summary records: + +- model/file label; +- residue count; +- number of residues with finite phi/psi; +- Rama8000 outlier count; +- mean and minimum pLDDT; +- declared prediction source; +- coordinate-file MD5 hash when available. + +The two downloadable CSV files therefore preserve both residue-level ensemble +results and model-level provenance. + +## Interpretation + +A useful pattern is a residue with: + +- large phi/psi spread; +- changing Rama8000 category; +- and/or large pLDDT variability across seeds. + +That combination identifies a position worth inspecting in the 2D plot and 3D +structure. It does **not** by itself identify a biologically relevant +conformational switch. + +Conversely, high pLDDT in every model does not imply that all models agree on +the same backbone geometry. RamplotR keeps confidence and geometric agreement +as separate measurements. + +## Current limitations + +- models are matched by residue identity rather than sequence realignment; +- AF3 multi-model confidence sidecars are not yet supported; +- ensemble members are not superposed or clustered in 3D yet; +- the current map summarizes local backbone variability, not global domain + motion; +- prediction ensembles are not a substitute for experimental ensemble or + dynamics data. + +Pairwise structural differences can be explored separately in the Compare tab, +including the Conformational Change Explorer. diff --git a/shinyRam/R/ensemble.R b/shinyRam/R/ensemble.R index 1d35b975..6b6cf08e 100644 --- a/shinyRam/R/ensemble.R +++ b/shinyRam/R/ensemble.R @@ -39,13 +39,34 @@ ram_ensemble_summary <- function(models) { membership <- matrix(FALSE,nrow=n,ncol=nmodels) for(i in seq_along(models)) membership[,i] <- !is.na(match(all_keys,keys[[i]])) - if(!n) return(data.frame(chain=character(),resi=integer(), + + optional_all <- function(field) + all(vapply(models,function(tbl) field %in% names(tbl),logical(1))) + + if(!n) { + empty <- data.frame(chain=character(),resi=integer(), insertion_code=character(),resn=character(),models_present=integer(), phi_models=integer(),psi_models=integer(), phi_mean=numeric(),phi_sd=numeric(),psi_mean=numeric(),psi_sd=numeric(), classified_models=integer(),class_consistency=numeric(), region_mode=character(),changes_class=logical(), - stringsAsFactors=FALSE)) + stringsAsFactors=FALSE) + if(optional_all("rama8000_region")) { + empty$rama8000_models <- integer() + empty$rama8000_consistency <- numeric() + empty$rama8000_mode <- character() + empty$rama8000_changes <- logical() + } + if(optional_all("plddt")) { + empty$plddt_models <- integer() + empty$plddt_mean <- numeric() + empty$plddt_sd <- numeric() + empty$plddt_min <- numeric() + empty$plddt_max <- numeric() + } + return(empty) + } + get <- function(col, type="numeric") { result <- matrix(if(type=="character") NA_character_ else NA_real_, nrow=n,ncol=nmodels) @@ -56,45 +77,147 @@ ram_ensemble_summary <- function(models) { } result } - phi <- get("phi");psi <- get("psi");region <- get("region","character") - # Use the first occurrence of each ID across all models, not row number. + phi <- get("phi"); psi <- get("psi"); region <- get("region","character") + ids <- do.call(rbind,lapply(models,function(x) x[,required[1:4],drop=FALSE])) - first <- !duplicated(unlist(keys,use.names=FALSE)) + flat_keys <- unlist(keys,use.names=FALSE) + first <- !duplicated(flat_keys) info <- ids[first,,drop=FALSE] - info <- info[match(all_keys,unlist(keys,use.names=FALSE)[first]),,drop=FALSE] + info <- info[match(all_keys,flat_keys[first]),,drop=FALSE] + calc <- function(matrix_values) { result <- t(vapply(seq_len(n),function(i) ram_ensemble_circular(matrix_values[i,]),numeric(2))) colnames(result) <- c("mean","sd") result } - ph <- calc(phi);ps <- calc(psi) - region_mode <- vapply(seq_len(n),function(i) { - values <- region[i,]; values <- values[!is.na(values)] - if(!length(values)) NA_character_ else - names(sort(table(values),decreasing=TRUE))[[1L]] - },character(1)) - agreement <- vapply(seq_len(n),function(i) { - values <- region[i,];values <- values[!is.na(values)] - if(!length(values)) NA_real_ else max(table(values))/length(values) - },numeric(1)) + mode_and_consistency <- function(matrix_values) { + mode <- vapply(seq_len(n),function(i) { + values <- matrix_values[i,]; values <- values[!is.na(values)] + if(!length(values)) NA_character_ else + names(sort(table(values),decreasing=TRUE))[[1L]] + },character(1)) + consistency <- vapply(seq_len(n),function(i) { + values <- matrix_values[i,]; values <- values[!is.na(values)] + if(!length(values)) NA_real_ else max(table(values))/length(values) + },numeric(1)) + list(mode=mode,consistency=consistency, + count=as.integer(rowSums(!is.na(matrix_values)))) + } + + ph <- calc(phi); ps <- calc(psi) + native <- mode_and_consistency(region) out <- data.frame(info, models_present=as.integer(rowSums(membership)), phi_models=as.integer(rowSums(is.finite(phi))), psi_models=as.integer(rowSums(is.finite(psi))), phi_mean=ph[,"mean"],phi_sd=ph[,"sd"], psi_mean=ps[,"mean"],psi_sd=ps[,"sd"], - classified_models=as.integer(rowSums(!is.na(region))), - class_consistency=as.numeric(agreement),region_mode=region_mode, - changes_class=!is.na(agreement) & agreement<1, + classified_models=native$count, + class_consistency=as.numeric(native$consistency), + region_mode=native$mode, + changes_class=!is.na(native$consistency) & native$consistency<1, stringsAsFactors=FALSE,check.names=FALSE) - out[order(-as.integer(out$changes_class), - -pmax(replace(out$phi_sd,is.na(out$phi_sd),0), - replace(out$psi_sd,is.na(out$psi_sd),0)), + + if(optional_all("rama8000_region")) { + standard <- get("rama8000_region","character") + stat <- mode_and_consistency(standard) + out$rama8000_models <- stat$count + out$rama8000_consistency <- stat$consistency + out$rama8000_mode <- stat$mode + out$rama8000_changes <- !is.na(stat$consistency) & stat$consistency<1 + } + + if(optional_all("plddt")) { + confidence <- get("plddt") + finite_count <- rowSums(is.finite(confidence)) + safe_stat <- function(fun) vapply(seq_len(n),function(i) { + values <- confidence[i,is.finite(confidence[i,])] + if(!length(values)) NA_real_ else fun(values) + },numeric(1)) + out$plddt_models <- as.integer(finite_count) + out$plddt_mean <- safe_stat(base::mean) + out$plddt_sd <- vapply(seq_len(n),function(i) { + values <- confidence[i,is.finite(confidence[i,])] + if(length(values)<2L) NA_real_ else stats::sd(values) + },numeric(1)) + out$plddt_min <- safe_stat(base::min) + out$plddt_max <- safe_stat(base::max) + } + + spread <- pmax(replace(out$phi_sd,is.na(out$phi_sd),0), + replace(out$psi_sd,is.na(out$psi_sd),0)) + standard_change <- if("rama8000_changes" %in% names(out)) + as.integer(out$rama8000_changes) else rep(0L,nrow(out)) + out[order(-standard_change,-as.integer(out$changes_class),-spread, out$chain,out$resi,out$insertion_code),,drop=FALSE] } +ram_prediction_ensemble_analyze <- function(pdbs, classifier, source, + labels = NULL, + max_models = 30L) { + if(!is.list(pdbs) || length(pdbs) < 2L) + stop("A prediction ensemble requires at least two predicted structures.") + permitted <- c("alphafold2","esmfold","other_prediction") + if(length(source)!=1L || !source %in% permitted) + stop("Prediction ensembles currently support AF2/ColabFold, ESMFold or other pLDDT-in-B-factor models.") + max_models <- suppressWarnings(as.integer(max_models)) + if(length(max_models)!=1L || is.na(max_models) || + max_models < 2L || max_models > 30L) + stop("Prediction ensemble max_models must be between 2 and 30.") + count <- min(length(pdbs),max_models) + if(is.null(labels)) labels <- paste("Model",seq_along(pdbs)) + labels <- as.character(labels) + if(length(labels)!=length(pdbs) || any(!nzchar(labels))) + stop("Every prediction model needs a label.") + + models <- lapply(seq_len(count),function(i) { + pdb <- ram_model_at(pdbs[[i]],1L) + torsions <- ram_extract_torsions(pdb) + classified <- classifier(torsions) + confidence <- ram_prediction_from_atoms(pdb,torsions,source) + keys <- ram_prediction_key(classified$chain,classified$resi, + classified$insertion_code) + confidence_keys <- ram_prediction_key(confidence$chain,confidence$resi, + confidence$insertion_code) + if(anyDuplicated(confidence_keys)) + stop("Prediction residue identifiers must be unique.") + classified$plddt <- confidence$plddt[match(keys,confidence_keys)] + classified$confidence_category <- ram_plddt_category(classified$plddt) + classified$model_label <- labels[[i]] + classified + }) + summary <- ram_ensemble_summary(models) + model_summary <- do.call(rbind,lapply(seq_len(count),function(i) { + table <- models[[i]] + finite_angles <- is.finite(table$phi) & is.finite(table$psi) + data.frame( + model=labels[[i]], + residues=nrow(table), + finite_phi_psi=sum(finite_angles), + rama8000_outliers=if ("rama8000_region" %in% names(table)) + sum(table$rama8000_region=="Outlier",na.rm=TRUE) else NA_integer_, + plddt_mean=if ("plddt" %in% names(table) && any(is.finite(table$plddt))) + base::mean(table$plddt[is.finite(table$plddt)]) else NA_real_, + plddt_min=if ("plddt" %in% names(table) && any(is.finite(table$plddt))) + base::min(table$plddt[is.finite(table$plddt)]) else NA_real_, + stringsAsFactors=FALSE + ) + })) + list( + summary=summary, + model_summary=model_summary, + models=models, + labels=labels[seq_len(count)], + analyzed_models=count, + available_models=length(pdbs), + common_residues=sum(summary$models_present==count), + source=source, + limited=count1L) tagList( + tags$h4("Models stored in this structure"), + tags$p(class="ram-confidence-explainer", + "Model variation is matched by chain, residue and insertion code. Circular statistics correctly handle the -180°/180° boundary; models with missing coordinates contribute only observed angles."), + tags$div(class="ram-ensemble-actions", + actionButton("calculateEnsemble","Analyse structural models", + class="btn-primary btn-sm"), + downloadButton("downloadEnsemble","Export structural ensemble CSV") + ), + uiOutput("ensembleResultSummary"), + tags$div(class="ram-residue-table",DT::DTOutput("ensembleRows")) ), - uiOutput("ensembleResultSummary"), - tags$div(class="ram-residue-table",DT::DTOutput("ensembleRows")) + if(is_prediction) tagList( + if(structure$nmodels>1L) tags$hr(), + tags$div(class="ram-prediction-ensemble-head", + tags$h4("Prediction ensemble"), + tags$p(class="ram-confidence-explainer", + "Upload independently generated AF2/ColabFold, ESMFold or other pLDDT-in-B-factor models. RamplotR compares model-to-model geometry and confidence; this variation is prediction uncertainty/heterogeneity, not experimental dynamics.") + ), + tags$div(class="ram-prediction-ensemble-controls", + selectInput("predictionEnsembleSource","Prediction model type", + choices=c("AlphaFold 2 / ColabFold"="alphafold2", + "ESMFold"="esmfold", + "Other model with pLDDT in B-factor"="other_prediction"), + selected=if(structure$declared_source %in% + c("esmfold","other_prediction")) structure$declared_source + else "alphafold2", + selectize=FALSE), + fileInput("predictionEnsembleFiles", + "Additional prediction models", + multiple=TRUE, + accept=c(".pdb",".ent",".cif",".mmcif",".mcif")), + if(!identical(structure$declared_source,"alphafold3")) + checkboxInput("includeLoadedPrediction", + paste("Include currently loaded model:",structure$name),value=TRUE) + else + tags$p(class="ram-confidence-warning", + "The loaded AlphaFold 3 model is not auto-added: matching atom-confidence JSON is required for ensemble confidence analysis."), + actionButton("calculatePredictionEnsemble", + "Analyse prediction ensemble",class="btn-primary btn-sm") + ), + tags$p(class="ram-field-hint", + "AlphaFold 3 ensembles are not accepted in this first version because per-model atom confidence needs its matching JSON sidecar; they are not silently treated as AF2."), + uiOutput("predictionEnsembleSummary"), + uiOutput("predictionEnsembleTrack"), + tags$div(class="ram-residue-table", + DT::DTOutput("predictionEnsembleRows")), + tags$div(class="ram-ensemble-actions", + downloadButton("downloadPredictionEnsemble", + "Export prediction ensemble CSV"), + downloadButton("downloadPredictionEnsembleModels", + "Export model summary CSV") + ) + ) ) ) }) @@ -1872,11 +1924,13 @@ server <- function(input, output, session) { withProgress(message="Analysing compatible ensemble models",value=0.2,{ result <- tryCatch( ram_ensemble_analyze(structure$pdb,max_models=min(30L,structure$nmodels), - classifier=function(torsions) - ram_classify_torsions(torsions, + classifier=function(torsions) { + classified <- ram_classify_torsions(torsions, reference_dir=file.path("static",input$bgtype), selected_reference=plot_reference(),mode=input$validationMode, - threshold_fn=ram_density_thresholds)), + threshold_fn=ram_density_thresholds) + ram_rama8000_classify(classified,file.path("static","rama8000")) + }), error=function(e) { showNotification(conditionMessage(e),type="error",duration=12) NULL @@ -1943,6 +1997,317 @@ server <- function(input, output, session) { file,row.names=FALSE,na="") ) + prediction_ensemble_input_key <- reactive({ + structure <- req(loaded()) + uploaded <- input$predictionEnsembleFiles + file_signature <- if(is.null(uploaded) || !nrow(uploaded)) "" else + paste(uploaded$name,uploaded$size,uploaded$type,uploaded$datapath, + sep=":",collapse="|") + include_loaded <- isTRUE(input$includeLoadedPrediction) && + !identical(structure$declared_source,"alphafold3") + paste( + if(is.null(input$predictionEnsembleSource)) "" else input$predictionEnsembleSource, + include_loaded, + if(include_loaded) current_model() else "", + file_signature, + sep="::" + ) + }) + + prediction_ensemble_matches <- reactive({ + value <- prediction_ensemble_results() + if(is.null(value)) return(NULL) + structure <- req(loaded()) + if(!identical(value$key,structure$key) || + !identical(value$mode,input$validationMode) || + !identical(value$reference,input$bgtype) || + !identical(value$background,input$background) || + !identical(value$input_key,prediction_ensemble_input_key())) + return(NULL) + value$result + }) + + observeEvent(loaded(), { + prediction_ensemble_results(NULL) + }, ignoreInit=TRUE) + + observeEvent(input$calculatePredictionEnsemble, { + structure <- req(loaded()) + source <- req(input$predictionEnsembleSource) + permitted <- c("alphafold2","esmfold","other_prediction") + if(!source %in% permitted) return() + + uploaded <- input$predictionEnsembleFiles + include_loaded <- isTRUE(input$includeLoadedPrediction) && + !identical(structure$declared_source,"alphafold3") + source_loaded <- if(identical(structure$declared_source,"alphafold_db")) + "alphafold2" else structure$declared_source + + if(include_loaded && !identical(source_loaded,source)) { + showNotification( + paste0("The loaded model is declared as ",source_loaded, + " but the ensemble is configured as ",source, + ". Choose the matching model type or exclude the loaded model."), + type="error",duration=14) + return() + } + + file_count <- if(is.null(uploaded)) 0L else nrow(uploaded) + if(file_count + as.integer(include_loaded) < 2L) { + showNotification( + "A prediction ensemble needs at least two models. Upload another model or include the loaded prediction.", + type="warning",duration=12) + return() + } + + withProgress(message="Analysing prediction ensemble",value=0.05,{ + pdbs <- list() + labels <- character() + hashes <- character() + structure_models <- integer() + input_roles <- character() + + if(include_loaded) { + selected_model <- current_model() + pdbs[[length(pdbs)+1L]] <- ram_model_at(structure$pdb,selected_model) + labels <- c(labels, + if(structure$nmodels>1L) + sprintf("%s [model %s]",structure$name,selected_model) + else structure$name) + hashes <- c(hashes, + if(is.character(structure$source_id) && + length(structure$source_id)==1L && + file.exists(structure$source_id)) + unname(tools::md5sum(structure$source_id)) else NA_character_) + structure_models <- c(structure_models,selected_model) + input_roles <- c(input_roles,"loaded") + } + + if(file_count) { + for(i in seq_len(file_count)) { + incProgress(0.35/max(1L,file_count), + detail=paste("Loading",uploaded$name[[i]])) + model <- tryCatch( + ram_load_structure( + path=uploaded$datapath[[i]], + original_name=uploaded$name[[i]] + ), + error=function(e) e + ) + if(inherits(model,"error")) { + showNotification( + paste(uploaded$name[[i]],conditionMessage(model),sep=": "), + type="error",duration=14) + return() + } + if(ram_model_count(model)!=1L) { + showNotification( + paste(uploaded$name[[i]], + "contains multiple structural models. Prediction-ensemble uploads must contain one model per file."), + type="error",duration=14) + return() + } + pdbs[[length(pdbs)+1L]] <- model + labels <- c(labels, + tools::file_path_sans_ext(basename(uploaded$name[[i]]))) + hashes <- c(hashes,unname(tools::md5sum(uploaded$datapath[[i]]))) + structure_models <- c(structure_models,1L) + input_roles <- c(input_roles,"uploaded") + } + } + + known_hashes <- hashes[!is.na(hashes) & nzchar(hashes)] + if(anyDuplicated(known_hashes)) { + showNotification( + "The ensemble contains duplicate coordinate files. Remove duplicate seeds/models before analysing agreement.", + type="error",duration=14) + return() + } + + result <- tryCatch( + ram_prediction_ensemble_analyze( + pdbs,source=source,labels=labels,max_models=30L, + classifier=function(torsions) { + classified <- ram_classify_torsions( + torsions, + reference_dir=file.path("static",input$bgtype), + selected_reference=plot_reference(), + mode=input$validationMode, + threshold_fn=ram_density_thresholds + ) + ram_rama8000_classify( + classified,file.path("static","rama8000")) + } + ), + error=function(e) { + showNotification(conditionMessage(e),type="error",duration=14) + NULL + } + ) + if(is.null(result)) return() + result$provenance <- data.frame( + model=result$labels, + source=result$source, + input_role=input_roles[seq_len(result$analyzed_models)], + structure_model=structure_models[seq_len(result$analyzed_models)], + coordinate_md5=hashes[seq_len(result$analyzed_models)], + stringsAsFactors=FALSE + ) + prediction_ensemble_results(list( + key=structure$key,mode=input$validationMode, + reference=input$bgtype,background=input$background, + input_key=prediction_ensemble_input_key(), + result=result + )) + incProgress(0.6,detail="Summarising model agreement") + }) + },ignoreInit=TRUE) + + output$predictionEnsembleSummary <- renderUI({ + result <- prediction_ensemble_matches() + if(is.null(result)) return(tags$p(class="ram-field-hint", + "Upload at least two compatible prediction models and run the ensemble analysis.")) + data <- result$summary + standard_changes <- if("rama8000_changes" %in% names(data)) + sum(data$rama8000_changes,na.rm=TRUE) else 0L + angular_variable <- sum( + pmax(data$phi_sd,data$psi_sd,na.rm=TRUE)>=20,na.rm=TRUE) + confidence_variable <- if("plddt_sd" %in% names(data)) + sum(is.finite(data$plddt_sd) & data$plddt_sd>=10) else 0L + tags$div( + tags$div(class="ram-confidence-metrics", + tags$span(class="ram-confidence-metric", + sprintf("%s models analysed",result$analyzed_models)), + tags$span(class="ram-confidence-metric", + sprintf("%s residues present in every model",result$common_residues)), + tags$span(class="ram-confidence-metric", + sprintf("%s residues with Rama8000 disagreement",standard_changes)), + tags$span(class="ram-confidence-metric", + sprintf("%s residues with ≥20° angular SD",angular_variable)), + tags$span(class="ram-confidence-metric", + sprintf("%s residues with pLDDT SD ≥10",confidence_variable)) + ), + tags$p(class="ram-confidence-explainer", + "These values quantify disagreement among prediction models/seeds. They do not demonstrate molecular motion or experimental conformational heterogeneity."), + if(result$limited) + tags$p(class="ram-confidence-warning", + "Only the first 30 models were analysed.") + ) + }) + + output$predictionEnsembleTrack <- renderUI({ + result <- prediction_ensemble_matches() + if(is.null(result) || !nrow(result$summary)) return(NULL) + data <- result$summary + spread <- pmax(data$phi_sd,data$psi_sd,na.rm=TRUE) + spread[!is.finite(data$phi_sd) & !is.finite(data$psi_sd)] <- NA_real_ + band <- ifelse(!is.finite(spread),"unavailable", + ifelse(spread<5,"stable", + ifelse(spread<15,"moderate", + ifelse(spread<30,"variable","high")))) + cells <- lapply(seq_len(nrow(data)),function(i) { + label <- paste0(data$resn[[i]]," ",data$chain[[i]],":", + data$resi[[i]],data$insertion_code[[i]]) + standard <- if("rama8000_changes" %in% names(data) && + isTRUE(data$rama8000_changes[[i]])) + " · Rama8000 category differs across models" else "" + tags$button(type="button", + class=paste("ram-ensemble-cell", + paste0("ram-ensemble-",band[[i]]), + if(nzchar(standard)) "has-standard-change" else ""), + "data-chain"=data$chain[[i]], + "data-resi"=data$resi[[i]], + "data-insertion"=data$insertion_code[[i]], + title=paste0(label," · angular SD ", + if(is.finite(spread[[i]])) sprintf("%.1f°",spread[[i]]) else "N/A", + if("plddt_mean" %in% names(data) && is.finite(data$plddt_mean[[i]])) + sprintf(" · mean pLDDT %.1f",data$plddt_mean[[i]]) else "", + standard), + "aria-label"=paste("Inspect",label,"from prediction ensemble") + ) + }) + tags$section(class="ram-ensemble-track-panel", + tags$div(class="ram-ensemble-track-head", + tags$strong("Prediction variability map"), + tags$span("max circular SD of φ or ψ per residue") + ), + tags$div(class="ram-ensemble-track",role="group", + "aria-label"="Prediction ensemble residue variability",cells), + tags$div(class="ram-ensemble-track-legend", + tags$span(class="ram-ensemble-stable","<5°"), + tags$span(class="ram-ensemble-moderate","5–15°"), + tags$span(class="ram-ensemble-variable","15–30°"), + tags$span(class="ram-ensemble-high","≥30°"), + tags$span(class="ram-ensemble-standard-mark", + "outline = Rama8000 disagreement")) + ) + }) + + output$predictionEnsembleRows <- DT::renderDT({ + result <- prediction_ensemble_matches() + req(result) + data <- result$summary + if(!nrow(data)) return(DT::datatable(data,rownames=FALSE)) + fields <- c("chain","resi","insertion_code","resn","models_present", + "phi_sd","psi_sd","rama8000_mode","rama8000_consistency", + "plddt_mean","plddt_sd","plddt_min","plddt_max") + fields <- fields[fields %in% names(data)] + shown <- data[,fields,drop=FALSE] + for(field in intersect(c("phi_sd","psi_sd","plddt_mean","plddt_sd", + "plddt_min","plddt_max"),names(shown))) + shown[[field]] <- round(shown[[field]],1L) + if("rama8000_consistency" %in% names(shown)) + shown$rama8000_consistency <- round(100*shown$rama8000_consistency,1L) + names(shown) <- c( + chain="Chain",resi="Residue",insertion_code="Ins.",resn="AA", + models_present="Models",phi_sd="φ SD (°)",psi_sd="ψ SD (°)", + rama8000_mode="Rama8000 mode", + rama8000_consistency="Rama8000 agreement (%)", + plddt_mean="pLDDT mean",plddt_sd="pLDDT SD", + plddt_min="pLDDT min",plddt_max="pLDDT max" + )[names(shown)] + DT::datatable(shown,rownames=FALSE,selection="single", + options=list(pageLength=12,scrollX=TRUE,autoWidth=FALSE,dom="ftip"), + class="compact stripe hover") + },server=FALSE) + + observeEvent(input$predictionEnsembleRows_rows_selected, { + data <- req(prediction_ensemble_matches())$summary + ix <- input$predictionEnsembleRows_rows_selected[[1L]] + if(!length(ix) || !is.finite(ix) || ix<1L || ix>nrow(data)) return() + row <- data[ix,,drop=FALSE] + select_from(list(chain=as.character(row$chain[[1L]]), + resi=as.integer(row$resi[[1L]]), + insertion_code=as.character(row$insertion_code[[1L]]))) + }) + + observeEvent(input$ramPredictionEnsemblePick, { + select_from(input$ramPredictionEnsemblePick) + },ignoreInit=TRUE) + + output$downloadPredictionEnsemble <- downloadHandler( + filename=function() safe_filename("prediction-ensemble-residues.csv"), + content=function(file) utils::write.csv( + req(prediction_ensemble_matches())$summary,file,row.names=FALSE,na="") + ) + output$downloadPredictionEnsembleModels <- downloadHandler( + filename=function() safe_filename("prediction-ensemble-models.csv"), + content=function(file) { + result <- req(prediction_ensemble_matches()) + models <- result$model_summary + if(!is.null(result$provenance)) { + provenance <- result$provenance + if(nrow(provenance)!=nrow(models)) + stop("Prediction ensemble provenance no longer matches model order.") + models$source <- provenance$source + models$input_role <- provenance$input_role + models$structure_model <- provenance$structure_model + models$coordinate_md5 <- provenance$coordinate_md5 + } + utils::write.csv(models,file,row.names=FALSE,na="") + } + ) + output$summary <- renderUI({ data <- displayed() eligible <- data[!is.na(data$region) & !data$resn %in% c("GLY", "PRO"), diff --git a/shinyRam/www/custom.js b/shinyRam/www/custom.js index 8a7000a2..88bc7bfb 100644 --- a/shinyRam/www/custom.js +++ b/shinyRam/www/custom.js @@ -167,6 +167,20 @@ // A delegated sequence click survives Shiny's HTML re-rendering, including // when the user switches chains or changes scientific reference datasets. + document.addEventListener("click", function (event) { + const button = event.target && event.target.closest && + event.target.closest(".ram-ensemble-cell"); + if (!button) return; + const resi = Number(button.dataset.resi); + if (!Number.isInteger(resi) || !window.Shiny || !window.Shiny.setInputValue) + return; + window.Shiny.setInputValue("ramPredictionEnsemblePick", { + chain: String(button.dataset.chain || ""), + resi, + insertion_code: String(button.dataset.insertion || "") + }, { priority: "event" }); + }); + document.addEventListener("click", function (event) { const button = event.target && event.target.closest && event.target.closest(".ram-seq-res"); diff --git a/shinyRam/www/styles.css b/shinyRam/www/styles.css index a6e11148..96b9fcf3 100644 --- a/shinyRam/www/styles.css +++ b/shinyRam/www/styles.css @@ -1163,6 +1163,49 @@ body > .container-fluid { max-width: none; padding: 0; } .ram-sequence-scroll-hint { display:none; } } +/* Prediction-ensemble controls and variability map. These encode model-to-model + disagreement only; the interface explicitly avoids presenting it as dynamics. */ +.ram-prediction-ensemble-head h4 { margin:0 0 4px; color:#184f57; } +.ram-prediction-ensemble-controls { + display:grid; grid-template-columns:minmax(190px,1fr) minmax(220px,1.4fr); + gap:9px 14px; align-items:end; margin:10px 0 8px; +} +.ram-prediction-ensemble-controls .form-group { margin-bottom:0; } +.ram-prediction-ensemble-controls .checkbox { margin:4px 0 0; } +.ram-prediction-ensemble-controls .btn { justify-self:start; } +.ram-ensemble-track-panel { margin:12px 0 15px; padding:11px 12px; + border:1px solid #dce9e6; border-radius:9px; background:#fbfdfc; } +.ram-ensemble-track-head { display:flex; align-items:baseline; gap:8px; + flex-wrap:wrap; color:#5e777d; font-size:10px; } +.ram-ensemble-track-head strong { color:#245c63; font-size:12px; } +.ram-ensemble-track { display:flex; gap:1px; overflow-x:auto; + padding:7px 2px 8px; margin-top:5px; scrollbar-width:thin; } +.ram-ensemble-cell { position:relative; width:7px; min-width:7px; height:22px; + border:0; border-radius:2px; padding:0; cursor:pointer; } +.ram-ensemble-stable { background:#7fb9af; } +.ram-ensemble-moderate { background:#d5bd61; } +.ram-ensemble-variable { background:#dc8d4c; } +.ram-ensemble-high { background:#cb5848; } +.ram-ensemble-unavailable { background:#c9d4d5; } +.ram-ensemble-cell.has-standard-change { + outline:1px solid #71333a; outline-offset:-1px; +} +.ram-ensemble-cell:hover,.ram-ensemble-cell:focus-visible { + outline:2px solid #135f66; outline-offset:1px; z-index:2; +} +.ram-ensemble-track-legend { display:flex; flex-wrap:wrap; gap:5px; + align-items:center; color:#60777d; font-size:9px; } +.ram-ensemble-track-legend span { padding:3px 7px; border-radius:999px; + border:1px solid rgba(40,70,72,.12); color:#294f55; font-weight:700; } +.ram-ensemble-track-legend .ram-ensemble-stable { background:#dcefeb; } +.ram-ensemble-track-legend .ram-ensemble-moderate { background:#f1e8bd; } +.ram-ensemble-track-legend .ram-ensemble-variable { background:#f2c995; } +.ram-ensemble-track-legend .ram-ensemble-high { background:#e99d8d; } +.ram-ensemble-standard-mark { background:#fff; box-shadow:inset 0 0 0 1px #71333a; } +@media(max-width:760px) { + .ram-prediction-ensemble-controls { grid-template-columns:1fr; } +} + /* Compact residue-level comparison track. Colour is ordinal only: the exact wrapped delta-phi/delta-psi values remain the scientific measurements. */ .ram-change-explorer { margin:12px 0 16px; padding:12px 13px; diff --git a/tests/ensemble.R b/tests/ensemble.R index 6bfa1a41..64ddf3b9 100644 --- a/tests/ensemble.R +++ b/tests/ensemble.R @@ -22,6 +22,22 @@ assert(one$phi_models==2L && one$psi_models==2L && assert(two$changes_class && two$class_consistency==0.5 && two$classified_models==2L, "Classification changes across models must be explicit.") + +m$rama8000_region <- c("Favored","Allowed") +m$plddt <- c(95,88) +other$rama8000_region <- ifelse(other$resi==1L,"Outlier","Allowed") +other$plddt <- ifelse(other$resi==1L,91,92) +extended <- ram_ensemble_summary(list(m,other)) +e1 <- extended[extended$resi==1L,,drop=FALSE] +e2 <- extended[extended$resi==2L,,drop=FALSE] +assert(e1$rama8000_changes && e1$rama8000_consistency==0.5 && + e1$rama8000_models==2L, + "Prediction ensembles must expose Rama8000 category disagreement.") +assert(!e2$rama8000_changes && e2$rama8000_consistency==1, + "Stable Rama8000 categories must remain explicit.") +assert(e1$plddt_models==2L && abs(e1$plddt_mean-93)<1e-12 && + e1$plddt_min==91 && e1$plddt_max==95, + "Prediction ensemble pLDDT summary changed.") only_one <- other[other$resi!=1L,,drop=FALSE] partial <- ram_ensemble_summary(list(m,only_one)) p <- partial[partial$resi==1L,,drop=FALSE] @@ -30,4 +46,64 @@ assert(p$models_present==1L && p$phi_models==1L && bad <- rbind(m,m[1,,drop=FALSE]) assert(inherits(try(ram_ensemble_summary(list(bad)),silent=TRUE), "try-error"),"Ambiguous IDs must be rejected.") +# Prediction-ensemble orchestration is tested with small stubs so this unit test +# remains independent of Bio3D and file I/O. +ram_model_at <- function(pdb,index=1L) pdb +ram_extract_torsions <- function(pdb) pdb$torsions +ram_prediction_key <- function(chain,resi,insertion_code="") { + paste(chain,resi,insertion_code,sep="\r") +} +ram_plddt_category <- function(value) { + ifelse(value>=90,"Very high",ifelse(value>=70,"Confident", + ifelse(value>=50,"Low","Very low"))) +} +ram_prediction_from_atoms <- function(pdb,torsions,source) { + data.frame(chain=torsions$chain,resi=torsions$resi, + insertion_code=torsions$insertion_code, + plddt=pdb$plddt,stringsAsFactors=FALSE) +} +make_prediction <- function(phi_shift=0,plddt=c(92,81)) list( + torsions=data.frame(chain="A",resi=1:2,insertion_code="", + resn=c("ALA","SER"),phi=c(-60,-80)+phi_shift, + psi=c(-45,150),stringsAsFactors=FALSE), + plddt=plddt +) +prediction_classifier <- function(torsions) { + torsions$region <- c("Favoured","Allowed") + torsions$rama8000_region <- if(torsions$phi[[1L]] < -55) + c("Favored","Allowed") else c("Allowed","Allowed") + torsions +} +prediction <- ram_prediction_ensemble_analyze( + list(make_prediction(0,c(94,80)),make_prediction(12,c(88,84))), + classifier=prediction_classifier,source="alphafold2", + labels=c("seed-1","seed-2") +) +assert(prediction$analyzed_models==2L && prediction$available_models==2L && + prediction$common_residues==2L, + "Prediction ensemble model/residue coverage changed.") +assert(identical(prediction$labels,c("seed-1","seed-2")) && + identical(as.character(prediction$model_summary$model), + c("seed-1","seed-2")), + "Prediction ensemble model labels must be preserved.") +p1 <- prediction$summary[prediction$summary$resi==1L,,drop=FALSE] +assert(p1$plddt_models==2L && abs(p1$plddt_mean-91)<1e-12 && + p1$rama8000_changes && p1$rama8000_consistency==0.5, + "Prediction ensemble must combine pLDDT and Rama8000 disagreement.") +assert(abs(prediction$model_summary$plddt_mean[[1L]]-87)<1e-12 && + prediction$model_summary$residues[[1L]]==2L, + "Per-model prediction provenance summary changed.") +assert(inherits(try(ram_prediction_ensemble_analyze( + list(make_prediction(),make_prediction()),prediction_classifier, + source="alphafold3"),silent=TRUE),"try-error"), + "AF3 must not be silently interpreted as a B-factor prediction ensemble.") +assert(inherits(try(ram_prediction_ensemble_analyze( + list(make_prediction()),prediction_classifier, + source="esmfold"),silent=TRUE),"try-error"), + "Prediction ensemble helper must require at least two models.") +assert(inherits(try(ram_prediction_ensemble_analyze( + list(make_prediction(),make_prediction()),prediction_classifier, + source="esmfold",max_models=31L),silent=TRUE),"try-error"), + "Prediction ensemble helper must enforce the documented 30-model maximum.") + message("Circular ensemble geometry and residue-alignment tests passed.") diff --git a/tests/prediction-browser.cjs b/tests/prediction-browser.cjs index a8a94a09..56d98669 100644 --- a/tests/prediction-browser.cjs +++ b/tests/prediction-browser.cjs @@ -12,6 +12,10 @@ const puppeteer = require("puppeteer-core"); ? line.slice(0,60) + " 87.00" + line.slice(66) : line).join("\n"); const fixture = path.resolve(output,"synthetic-esmfold.pdb"); fs.writeFileSync(fixture,pdb); + const pdbSeed2 = raw.split(/\r?\n/).map(line => line.startsWith("ATOM ") + ? line.slice(0,60) + " 67.00" + line.slice(66) : line).join("\n"); + const fixtureSeed2 = path.resolve(output,"synthetic-esmfold-seed2.pdb"); + fs.writeFileSync(fixtureSeed2,pdbSeed2); const atomLines = pdb.split("\n").filter(line => line.startsWith("ATOM ")); const unique = [...new Set(atomLines.map(line => line.slice(21,22) + "|" + line.slice(22,26).trim() + "|" + @@ -167,6 +171,57 @@ const puppeteer = require("puppeteer-core"); ".ram-confidence-metrics")?.textContent.includes("0.82")); await page.screenshot({path:path.join(output,"prediction-af3.png"), fullPage:true}); - console.log("Live ESMFold, AlphaFold2 and AlphaFold3 prediction views passed."); + + // Prediction ensemble: two separate model files, deliberately carrying + // different synthetic pLDDT values. This tests multi-file upload, + // aggregation, variability rendering and map -> residue selection. + await page.click('.nav-tabs a[data-value="summary"]'); + await page.waitForSelector("#ram-ensemble-panel",{timeout:20000}); + const ensemblePanel = await page.$("#ram-ensemble-panel"); + if (!(await page.evaluate(el=>el.open,ensemblePanel))) + await page.click("#ram-ensemble-panel > summary"); + await page.waitForSelector("#predictionEnsembleSource",{timeout:12000}); + await page.select("#predictionEnsembleSource","esmfold"); + const ensembleUpload = await page.$("#predictionEnsembleFiles"); + await ensembleUpload.uploadFile(fixture,fixtureSeed2); + await page.waitForFunction(() => { + const input=document.getElementById("predictionEnsembleFiles"); + return input && input.files && input.files.length===2; + },{timeout:10000}); + await new Promise(done=>setTimeout(done,1200)); + await page.click("#calculatePredictionEnsemble"); + await page.waitForFunction(() => { + const panel=document.querySelector("#predictionEnsembleSummary"); + return panel && panel.textContent.includes("2 models analysed") && + document.querySelectorAll(".ram-ensemble-cell").length>20 && + document.querySelectorAll("#predictionEnsembleRows tbody tr").length>0; + },{timeout:35000}); + const ensembleState=await page.evaluate(() => { + const summary=document.querySelector("#predictionEnsembleSummary").textContent; + const rows=Array.from(document.querySelectorAll( + "#predictionEnsembleRows tbody tr")); + const first=rows[0] ? Array.from(rows[0].querySelectorAll("td")) + .map(td=>td.textContent.trim()) : []; + return { + summary, + cells:document.querySelectorAll(".ram-ensemble-cell").length, + first + }; + }); + const disagreement = ensembleState.summary.match( + /(\d+) residues with pLDDT SD ≥10/); + assert.ok(disagreement && Number(disagreement[1]) > 0, + "Synthetic 20-point pLDDT differences must produce a nonzero disagreement count."); + assert.ok(ensembleState.cells>20, + "Prediction ensemble variability map should cover the protein."); + await page.click(".ram-ensemble-cell"); + await page.waitForFunction(() => + document.querySelector("#selectedResidueInfo strong"), + {timeout:15000}); + await page.screenshot({ + path:path.join(output,"prediction-ensemble.png"),fullPage:true + }); + + console.log("Live ESMFold, AlphaFold2, AlphaFold3 and prediction-ensemble views passed."); } finally {await browser.close();} })().catch(e=>{console.error(e);process.exitCode=1;});