Add structure-group conformational comparison

.github/workflows/scientific-tests.yml

@@ -57,6 +57,7 @@ jobs:
           Rscript tests/geometry.R
           Rscript tests/experimental.R
           Rscript tests/ensemble.R
+          Rscript tests/group-comparison.R
           Rscript tests/batch.R
         shell: bash
       - name: Verify AlphaFold and ESMFold confidence formats

README.md

@@ -19,7 +19,7 @@ RamplotR is an open-source R Shiny application that brings backbone geometry, a
 | **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. |
+| **Compare models** | Sequence-align two structures with the Conformational Change Explorer, or compare biologically defined structure sets (for example apo/holo or WT/mutant) using residue-level circular φ/ψ means, dispersion and between-group backbone shifts. |
 | **Publish or automate** | Export SVG and high-resolution PNG figures, filtered CSV tables and standalone HTML reports. Run the offline R command-line tool on individual files or a directory of structures. |
 
 Advanced analysis stays in collapsible panels or dedicated comparison/summary views, keeping the everyday 2D/3D inspection screen uncluttered.
@@ -106,6 +106,7 @@ The default RamplotR teal contour palette provides consistent, recognisable publ
 - [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)
+- [Structure-group conformational comparison](docs/group-conformation-comparison.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)

docs/group-conformation-comparison.md

@@ -0,0 +1,115 @@
+# Structure-group conformational comparison
+
+RamplotR can compare two **sets** of related protein structures at residue level.
+Typical examples include apo versus ligand-bound structures, wild type versus
+mutants, or experimental versus predicted model sets.
+
+The analysis is deliberately conformation-centred. It compares circular
+backbone phi/psi distributions after sequence alignment instead of reducing
+each structure pair to a single Cartesian RMSD.
+
+## Workflow
+
+1. Load a representative structure in RamplotR.
+2. Open **Compare** and expand **Compare groups of structures**.
+3. Select the protein chain that defines the reference residue coordinate
+   system.
+4. Define Group A and Group B labels.
+5. Optionally include the currently loaded structure in Group A.
+6. Upload one or more PDB/mmCIF files for each group.
+7. Run **Analyse groups**.
+
+Each uploaded file contributes model 1. RamplotR chooses the best matching
+protein chain independently for every file using sequence identity and
+reference-chain coverage. The default acceptance thresholds are 70% identity
+and 70% reference coverage; both can be changed before analysis.
+
+## Residue-level calculations
+
+Every accepted candidate chain is globally sequence-aligned to the selected
+reference chain. Candidate phi/psi values are then represented on the reference
+residue coordinate system. Insertions without a reference residue are not
+invented as comparable positions.
+
+For each residue and each group RamplotR calculates:
+
+- number of models contributing phi and psi;
+- circular mean phi and psi;
+- circular standard deviation of phi and psi;
+- modal Rama8000 category and its within-group consistency.
+
+The between-group effect is reported as wrapped differences between the two
+circular means:
+
+```
+delta_phi = wrap(phi_mean_B - phi_mean_A)
+delta_psi = wrap(psi_mean_B - psi_mean_A)
+backbone_shift = sqrt(delta_phi^2 + delta_psi^2)
+```
+
+Wrapping is performed independently across the -180/180-degree boundary.
+
+The combined backbone shift is a **navigation effect size**, not a statistical
+significance score and not a Cartesian distance.
+
+## Consistent-shift marker
+
+A residue is highlighted as a low-dispersion consistent shift when:
+
+- both groups contribute at least two finite phi and psi observations;
+- the between-group combined shift is at least 30 degrees; and
+- the largest within-group circular SD across phi and psi is at most 15
+  degrees.
+
+These thresholds intentionally identify residues worth inspecting. They do not
+establish a biological effect or a p-value. The exact group means, angular
+differences and within-group dispersion remain visible in the table.
+
+A second marker indicates when the modal Rama8000 category differs between the
+two groups.
+
+## Interpretation
+
+Group-level conformational differences can reflect genuine structural states,
+but can also arise from:
+
+- ligands or cofactors;
+- construct boundaries and engineered mutations;
+- crystal packing;
+- cryo-EM classification;
+- differences in experimental resolution;
+- prediction uncertainty;
+- different domain arrangements or oligomeric states.
+
+Sequence alignment alone therefore does not prove that two groups are
+biologically exchangeable. The exported member table records the automatically
+selected chain, sequence identity and coverage for every input structure so
+these assumptions can be audited.
+
+For prediction ensembles, use RamplotR's dedicated prediction-ensemble
+workflow when the main question is seed/model uncertainty. Group comparison is
+more appropriate when the researcher has already defined biologically
+meaningful sets such as apo/holo or WT/mutant.
+
+## Exports
+
+The group comparison can export:
+
+- one residue-level CSV containing circular means, SDs, wrapped differences,
+  combined displacement and Rama8000 mode changes;
+- one member CSV recording each structure label, selected chain, identity,
+  reference coverage, candidate coverage and aligned residue count.
+
+## Testing
+
+Unit tests cover:
+
+- automatic best-chain selection;
+- rejection of unrelated chains at explicit thresholds;
+- circular means across +179/-179 degrees;
+- wrapped between-group differences;
+- low-dispersion consistent-shift detection;
+- Rama8000 modal-category changes.
+
+The browser test additionally compares 1CRN against an identical uploaded
+1CRN group and requires every residue to remain in the smallest shift band.

shinyRam/R/group-comparison.R

@@ -0,0 +1,192 @@
+# Group-level conformational comparison.
+#
+# Structures are sequence-aligned to a single reference chain. Per-group
+# circular phi/psi summaries reuse the same statistics as ensemble analysis.
+# Between-group shifts are descriptive effect sizes for navigation; no
+# significance test is implied.
+
+ram_alignment_quality <- function(reference, candidate) {
+  pairing <- ram_align_residues(reference, candidate)
+  both <- !is.na(pairing$index_a) & !is.na(pairing$index_b)
+  aligned <- sum(both)
+  if (!aligned) return(list(
+    identity=0, reference_coverage=0, candidate_coverage=0,
+    aligned=0L, score=0
+  ))
+  ra <- toupper(as.character(reference$resn[pairing$index_a[both]]))
+  rb <- toupper(as.character(candidate$resn[pairing$index_b[both]]))
+  identity <- mean(ra == rb)
+  ref_cov <- aligned / max(1L,nrow(reference))
+  cand_cov <- aligned / max(1L,nrow(candidate))
+  list(
+    identity=identity,
+    reference_coverage=ref_cov,
+    candidate_coverage=cand_cov,
+    aligned=as.integer(aligned),
+    score=identity * sqrt(ref_cov*cand_cov)
+  )
+}
+
+ram_best_chain_match <- function(reference, candidate,
+                                 min_identity=0.30,
+                                 min_reference_coverage=0.50) {
+  if (!nrow(reference) || !nrow(candidate))
+    stop("Reference and candidate structures need protein residues.")
+  chain_values <- unique(as.character(candidate$chain))
+  if (!length(chain_values)) stop("Candidate structure has no protein chain.")
+  scored <- lapply(chain_values,function(chain) {
+    table <- candidate[as.character(candidate$chain)==chain,,drop=FALSE]
+    quality <- ram_alignment_quality(reference,table)
+    data.frame(
+      chain=chain,
+      identity=quality$identity,
+      reference_coverage=quality$reference_coverage,
+      candidate_coverage=quality$candidate_coverage,
+      aligned=quality$aligned,
+      score=quality$score,
+      stringsAsFactors=FALSE
+    )
+  })
+  scores <- do.call(rbind,scored)
+  order_index <- order(-scores$score,-scores$identity,
+                       -scores$reference_coverage,scores$chain)
+  scores <- scores[order_index,,drop=FALSE]
+  best <- scores[1L,,drop=FALSE]
+  if (!is.finite(best$identity) || best$identity < min_identity ||
+      !is.finite(best$reference_coverage) ||
+      best$reference_coverage < min_reference_coverage) {
+    stop(sprintf(
+      "No candidate chain met the minimum alignment criteria (identity %.0f%%, reference coverage %.0f%%). Best chain %s: %.1f%% identity, %.1f%% reference coverage.",
+      100*min_identity,100*min_reference_coverage,best$chain,
+      100*best$identity,100*best$reference_coverage
+    ))
+  }
+  list(chain=best$chain, metrics=best, all=scores)
+}
+
+ram_map_chain_to_reference <- function(reference, candidate,
+                                       source_label,
+                                       source_chain=unique(candidate$chain)[[1L]]) {
+  pairing <- ram_align_residues(reference,candidate)
+  both <- !is.na(pairing$index_a) & !is.na(pairing$index_b)
+  pairing <- pairing[both,,drop=FALSE]
+  if (!nrow(pairing)) return(reference[0,,drop=FALSE])
+  ref_idx <- pairing$index_a
+  cand_idx <- pairing$index_b
+  required <- c("chain","resi","insertion_code","resn","phi","psi","region")
+  if (!all(required %in% names(reference)) ||
+      !all(required %in% names(candidate)))
+    stop("Mapped structures require residue identifiers, phi/psi and regions.")
+  out <- reference[ref_idx,c("chain","resi","insertion_code","resn"),
+                   drop=FALSE]
+  out$phi <- candidate$phi[cand_idx]
+  out$psi <- candidate$psi[cand_idx]
+  out$region <- candidate$region[cand_idx]
+  optional <- c("rama8000_region","rama8000_group","rama8000_score","plddt")
+  for(field in optional)
+    if(field %in% names(candidate)) out[[field]] <- candidate[[field]][cand_idx]
+  out$source_label <- as.character(source_label)
+  out$source_chain <- as.character(source_chain)
+  out
+}
+
+ram_prepare_structure_group <- function(reference, structures, labels=NULL,
+                                        min_identity=0.30,
+                                        min_reference_coverage=0.50) {
+  if (!is.list(structures) || !length(structures))
+    stop("Each comparison group needs at least one structure.")
+  if (is.null(labels)) labels <- paste("Structure",seq_along(structures))
+  labels <- as.character(labels)
+  if (length(labels)!=length(structures) || any(!nzchar(labels)))
+    stop("Every group structure needs a non-empty label.")
+  mapped <- vector("list",length(structures))
+  model_info <- vector("list",length(structures))
+  for(i in seq_along(structures)) {
+    table <- structures[[i]]
+    best <- ram_best_chain_match(
+      reference,table,min_identity,min_reference_coverage)
+    chain <- table[as.character(table$chain)==best$chain,,drop=FALSE]
+    mapped[[i]] <- ram_map_chain_to_reference(
+      reference,chain,labels[[i]],best$chain)
+    model_info[[i]] <- data.frame(
+      model=labels[[i]],
+      chain=best$chain,
+      identity=best$metrics$identity,
+      reference_coverage=best$metrics$reference_coverage,
+      candidate_coverage=best$metrics$candidate_coverage,
+      aligned=best$metrics$aligned,
+      stringsAsFactors=FALSE
+    )
+  }
+  list(models=mapped, model_summary=do.call(rbind,model_info))
+}
+
+ram_group_conformation_compare <- function(reference, group_a, group_b,
+                                           label_a="Group A",
+                                           label_b="Group B") {
+  if (!is.list(group_a) || !length(group_a) ||
+      !is.list(group_b) || !length(group_b))
+    stop("Both groups need at least one mapped structure.")
+  a <- ram_ensemble_summary(group_a)
+  b <- ram_ensemble_summary(group_b)
+  key <- function(data)
+    paste(data$chain,data$resi,data$insertion_code,toupper(data$resn),sep="\r")
+  ka <- key(a); kb <- key(b)
+  all_keys <- unique(c(ka,kb))
+  ia <- match(all_keys,ka); ib <- match(all_keys,kb)
+  template <- rbind(
+    a[!duplicated(ka),c("chain","resi","insertion_code","resn"),drop=FALSE],
+    b[!duplicated(kb),c("chain","resi","insertion_code","resn"),drop=FALSE]
+  )
+  kt <- key(template)
+  ids <- template[match(all_keys,kt),,drop=FALSE]
+  pick <- function(data,index,field,missing=NA_real_) {
+    out <- rep(missing,length(index))
+    good <- !is.na(index)
+    if(any(good) && field %in% names(data)) out[good] <- data[[field]][index[good]]
+    out
+  }
+  out <- ids
+  numeric_fields <- c("phi_models","psi_models","phi_mean","phi_sd",
+                      "psi_mean","psi_sd","rama8000_models",
+                      "rama8000_consistency")
+  char_fields <- c("rama8000_mode")
+  for(field in numeric_fields) {
+    out[[paste0("a_",field)]] <- pick(a,ia,field)
+    out[[paste0("b_",field)]] <- pick(b,ib,field)
+  }
+  for(field in char_fields) {
+    out[[paste0("a_",field)]] <- pick(a,ia,field,NA_character_)
+    out[[paste0("b_",field)]] <- pick(b,ib,field,NA_character_)
+  }
+  out$delta_phi <- ram_angular_difference(out$a_phi_mean,out$b_phi_mean)
+  out$delta_psi <- ram_angular_difference(out$a_psi_mean,out$b_psi_mean)
+  out$angular_displacement <- ram_backbone_angular_displacement(
+    out$delta_phi,out$delta_psi)
+  out$shift_band <- ram_backbone_shift_band(out$angular_displacement)
+  max_finite <- function(...) {
+    values <- cbind(...)
+    apply(values,1L,function(row) {
+      row <- row[is.finite(row)]
+      if(!length(row)) NA_real_ else max(row)
+    })
+  }
+  out$max_within_group_sd <- max_finite(
+    out$a_phi_sd,out$a_psi_sd,out$b_phi_sd,out$b_psi_sd)
+  enough <- out$a_phi_models>=2 & out$a_psi_models>=2 &
+            out$b_phi_models>=2 & out$b_psi_models>=2
+  out$consistent_shift <- enough &
+    is.finite(out$angular_displacement) &
+    out$angular_displacement>=30 &
+    is.finite(out$max_within_group_sd) &
+    out$max_within_group_sd<=15
+  out$rama8000_mode_changed <- !is.na(out$a_rama8000_mode) &
+    !is.na(out$b_rama8000_mode) &
+    out$a_rama8000_mode != out$b_rama8000_mode
+  out$group_a <- label_a
+  out$group_b <- label_b
+  out[order(-as.integer(out$consistent_shift),
+            -replace(out$angular_displacement,
+                     !is.finite(out$angular_displacement),-Inf),
+            out$chain,out$resi,out$insertion_code),,drop=FALSE]
+}