{"id":13926,"date":"2026-01-28T15:31:23","date_gmt":"2026-01-28T15:31:23","guid":{"rendered":"https:\/\/www.blopig.com\/blog\/?p=13926"},"modified":"2026-02-03T16:20:35","modified_gmt":"2026-02-03T16:20:35","slug":"advanced-pymol-visualization-for-weighted-structural-ensembles-part-1-ensemble-comparison","status":"publish","type":"post","link":"https:\/\/www.blopig.com\/blog\/2026\/01\/advanced-pymol-visualization-for-weighted-structural-ensembles-part-1-ensemble-comparison\/","title":{"rendered":"Advanced PyMOL Visualization for Weighted Structural Ensembles (Part 1): Ensemble Comparison"},"content":{"rendered":"\n

When working with structural ensembles from molecular dynamics, AlphaFold2 subsampling, or ensemble reweighting against experimental data, you quickly run into visualization problems. Many of these problems standard PyMOL tutorials don’t address: what do you do when there’s no single reference structure?<\/strong><\/p>\n\n\n\n

In this two-part series, I’ll share the PyMOL techniques I’ve developed for visualizing weighted ensembles where multiple conformational states coexist. Part 1 covers reference state handling, RMSD-based coloring, and cluster visualization. Part 2 will tackle efficient SASA surface generation for large ensembles. To the best of my knowledge, this is the most advanced PyMOL guide EVER.<\/strong><\/p>\n\n\n\n

The code snippets here are extracted from full scripts attached at the end of this post. All examples use two systems: TeaA<\/strong> (a membrane transporter with distinct open\/closed states) and MoPrP<\/strong> (mouse Prion Protein with partially unfolded forms).<\/p>\n\n\n\n\n\n\n\n

The Multi-Reference Problem<\/h2>\n\n\n\n

Traditional structure visualization assumes you have one<\/em> reference to align everything against. But ensemble methods often produce populations spanning multiple metastable states. For TeaA, we have open and closed conformations. For MoPrP, we have a folded state plus two partially unfolded forms (PUF1 and PUF2).<\/p>\n\n\n\n

Trying to show “how well our ensemble matches the data” requires thinking carefully about:<\/p>\n\n\n\n

    \n
  1. How to display multiple references without visual clutter<\/li>\n\n\n\n
  2. How to align structures when the reference itself is an ensemble<\/li>\n\n\n\n
  3. How to color structures by their conformational similarity to different states<\/li>\n<\/ol>\n\n\n\n

    Distinguishing Reference States with Cartoon Styles<\/h2>\n\n\n\n

    When overlaying multiple reference structures, color alone isn’t enough\u2014especially in grayscale-friendly publications. For multiple references, cartoon dash<\/strong> is a great choice: provides solid cartoons for the ensemble structures of interest. In my experience this is often clearer than transparency.<\/p>\n\n\n\n

    # Load reference structures\ncmd.load(\"TeaA_ref_closed_state.pdb\", \"closed_ref\")\ncmd.load(\"TeaA_ref_open_state.pdb\", \"open_ref\")\n\n# Style references as dashed cartoons\ncmd.show_as(\"cartoon\", \"closed_ref\")\ncmd.cartoon(\"dash\", \"closed_ref\")\ncmd.set(\"dash_width\", 0.1, \"closed_ref\")\ncmd.color(\"white\", \"closed_ref\")\ncmd.set(\"cartoon_transparency\", 0.4, \"closed_ref\")\n\ncmd.show_as(\"cartoon\", \"open_ref\")\ncmd.cartoon(\"dash\", \"open_ref\")\ncmd.set(\"dash_width\", 0.1, \"open_ref\")\ncmd.color(\"black\", \"open_ref\")\ncmd.set(\"cartoon_transparency\", 0.4, \"open_ref\")\n\n# Show all states if reference is an NMR ensemble\ncmd.set(\"all_states\", 1, \"closed_ref\")<\/pre>\n\n\n\n
    The dashed style immediately signals “this is a reference” while the transparency prevents it from dominating the visualization. For two-state systems like TeaA, I use a white\/black or blue\/orange color scheme that remains distinguishable even without color.
    <\/p>\n\n\n\n
    \"\"<\/a><\/figure>\n\n\n\n
    Left: Grayscale-friendly visualization using dash width variation. Right: Color-coded open (orange) and closed (blue) references.<\/em><\/p>\n\n\n\n
    Kabsch Alignment to Frame-Averaged Coordinates<\/h2>\n\n\n\n
    When your reference is itself an ensemble (like an NMR structure or MD trajectory), aligning to a single frame introduces bias. Instead, compute the mean coordinates<\/strong> across all reference frames and align to that.<\/p>\n\n\n\n
    import numpy as np\nfrom pymol import cmd\n\ndef get_reference_target_coords(reference_pdb, align_selection):\n    \"\"\"\n    Get target alignment coordinates from reference PDB.\n    If multi-state, returns frame-averaged coordinates.\n    \"\"\"\n    ref_obj = \"__temp_ref__\"\n    cmd.load(reference_pdb, ref_obj)\n    num_states = cmd.count_states(ref_obj)\n    \n    sel = f\"{ref_obj} and {align_selection}\"\n    \n    if num_states > 1:\n        # Compute frame average for multi-state references\n        sum_coords = None\n        count = 0\n        \n        for i in range(1, num_states + 1):\n            coords = cmd.get_coords(sel, state=i)\n            if coords is not None and len(coords) > 0:\n                if sum_coords is None:\n                    sum_coords = np.zeros_like(coords)\n                sum_coords += coords\n                count += 1\n        \n        target_coords = sum_coords \/ count if count > 0 else None\n    else:\n        target_coords = cmd.get_coords(sel, state=1)\n    \n    cmd.delete(ref_obj)\n    return target_coords<\/pre>\n\n\n\n
    The alignment selection matters enormously. For TeaA, I align on stable secondary structure elements that don’t change between open\/closed:
    <\/p>\n\n\n\n
    # TeaA: align on stable helices, not the mobile domains\nALIGN_SELECTION = \"(resi 226-256 or resi 16-27 or resi 66-69) and name CA\"<\/pre>\n\n\n\n
    For MoPrP, I use the stable helix region (residues 79-93) that remains folded even in partially unfolded forms:

    <\/p>\n\n\n\n
    # MoPrP: align on the stable helix core\nALIGN_START = 79\nALIGN_END = 93<\/pre>\n\n\n\n
    Now we can align every frame of our trajectory to these averaged reference coordinates using the Kabsch algorithm:

    <\/p>\n\n\n\n
    def kabsch_alignment(mobile_coords, reference_coords):\n    \"\"\"Kabsch algorithm to find optimal rotation matrix.\"\"\"\n    mobile_centered = mobile_coords - mobile_coords.mean(axis=0)\n    reference_centered = reference_coords - reference_coords.mean(axis=0)\n    \n    H = np.dot(mobile_centered.T, reference_centered)\n    U, S, Vt = np.linalg.svd(H)\n    R = np.dot(Vt.T, U.T)\n    \n    # Handle reflection case\n    if np.linalg.det(R) < 0:\n        Vt[-1, :] *= -1\n        R = np.dot(Vt.T, U.T)\n    \n    t = reference_coords.mean(axis=0) - np.dot(mobile_coords.mean(axis=0), R.T)\n    return R, t\n\n\ndef align_trajectory_to_coords(obj_name, target_coords, align_selection):\n    \"\"\"Align every state of a PyMOL object to target coordinates.\"\"\"\n    num_states = cmd.count_states(obj_name)\n    sel_align = f\"{obj_name} and {align_selection}\"\n    \n    for i in range(1, num_states + 1):\n        mobile_align_coords = cmd.get_coords(sel_align, state=i)\n        R, t = kabsch_alignment(mobile_align_coords, target_coords)\n        \n        # Transform ALL atoms, not just alignment selection\n        all_coords = cmd.get_coords(obj_name, state=i)\n        new_coords = np.dot(all_coords, R.T) + t\n        cmd.load_coords(new_coords, obj_name, state=i)<\/pre>\n\n\n\n
    RMSD-Based Coloring of Top Weighted Structures<\/h2>\n\n\n\n
    For weighted ensembles, you often want to show the highest-weight structures colored by their RMSD to a reference state. This reveals whether your top structures cluster near one conformational state or span multiple states.<\/p>\n\n\n\n
    def show_top_structures_by_rmsd(obj_name, weights, n=20, reference_obj=None):\n    \"\"\"\n    Extract top N structures by weight, color by RMSD to reference.\n    \"\"\"\n    top_indices = np.argsort(weights)[::-1][:n]\n    \n    # Calculate RMSD to reference for each top structure\n    rmsds = []\n    for idx in top_indices:\n        state = int(idx) + 1\n        temp_name = f\"temp_rmsd_{state}\"\n        cmd.create(temp_name, obj_name, state, 1)\n        \n        # Align and calculate RMSD\n        cmd.align(temp_name, reference_obj, cycles=0)\n        rmsd = cmd.rms_cur(f\"{temp_name} and name CA\", \n                          f\"{reference_obj} and name CA\")\n        rmsds.append(rmsd)\n        cmd.delete(temp_name)\n    \n    # Create individual objects colored by RMSD\n    rmsd_min, rmsd_max = min(rmsds), max(rmsds)\n    \n    for i, idx in enumerate(top_indices):\n        state = int(idx) + 1\n        name = f\"{obj_name}_top_{i+1:02d}\"\n        cmd.create(name, obj_name, state, 1)\n        \n        # Store RMSD in B-factor for spectrum coloring\n        cmd.alter(name, f\"b={rmsds[i]}\")\n        \n        # Set transparency based on weight\n        w = weights[idx]\n        transparency = max(0.1, min(0.9, 0.5 - w * 2))\n        cmd.set(\"cartoon_transparency\", transparency, name)\n        \n        # Visualization settings\n        cmd.show(\"cartoon\", name)\n        cmd.cartoon(\"tube\", name)\n        cmd.set(\"cartoon_tube_radius\", 0.2, name)\n        \n        # Color by RMSD: cyan (low\/similar) -> yellow (high\/different)\n        cmd.spectrum(\"b\", \"cyan white grey yellow\", name, \n                    minimum=rmsd_min, maximum=rmsd_max)\n    \n    # Group all top structures\n    members = \" \".join([f\"{obj_name}_top_{i+1:02d}\" for i in range(n)])\n    cmd.group(f\"{obj_name}_top_{n}_group\", members)<\/pre>\n\n\n\n
    For TeaA, this reveals whether the highest-weight structures favor the open or closed conformation:
    <\/p>\n\n\n\n
    \"\"<\/a><\/figure>\n\n\n\n
    Top weighted structures from different fitting methods, colored by RMSD to the open state. Cyan indicates structures similar to the open conformation; yellow indicates structures closer to closed.<\/em><\/p>\n\n\n\n
    Visualizing Cluster Mean Structures with Per-Residue Deviation<\/h2>\n\n\n\n
    When you have clustered your ensemble (e.g., into Folded, PUF1, PUF2 states), showing the mean structure<\/strong> of each cluster with per-residue RMSD coloring<\/strong> communicates both the average conformation and its internal variability.<\/p>\n\n\n\n
    First, split your trajectory into individual frames:<\/p>\n\n\n\n
    <\/pre>\n\n\n\n
    Then compute per-residue RMSD and store in B-factors:
    <\/p>\n\n\n\n
    # Accumulate squared deviations for mean structure RMSD\n    sum_sq_dev = {}\n    n_frames = len(frame_objs)\n    \n    for obj in frame_objs:\n        model = cmd.get_model(f\"{obj} and name CA\")\n        \n        for atom in model.atom:\n            resi = atom.resi.strip()\n            curr_pos = np.array(atom.coord)\n            \n            if resi in ref_coords_dict:\n                ref_pos = ref_coords_dict[resi]\n                dist = np.linalg.norm(curr_pos - ref_pos)\n                \n                # Accumulate for ensemble RMSD\n                if resi not in sum_sq_dev:\n                    sum_sq_dev[resi] = 0.0\n                sum_sq_dev[resi] += dist ** 2\n        \n        # Apply deviation to this frame's B-factors\n        stored.dev_map = {resi: np.linalg.norm(\n            np.array(cmd.get_model(f\"{obj} and name CA and resi {resi}\").atom[0].coord) - \n            ref_coords_dict.get(resi, np.zeros(3))\n        ) for resi in ref_coords_dict}\n        \n        cmd.alter(obj, \"b=stored.dev_map.get(resi.strip(), 0.0)\")<\/pre>\n\n\n\n
    Finally, create the mean structure with ensemble RMSD coloring:
    <\/p>\n\n\n\n
    # Create mean structure object\n    mean_obj = f\"{obj_name}_mean\"\n    cmd.create(mean_obj, frame_objs[0])\n    \n    # Compute mean coordinates for ALL atoms\n    mean_coords = compute_all_atom_mean(obj_name)\n    cmd.load_coords(mean_coords, mean_obj, state=1)\n    \n    # Calculate per-residue RMSD and apply to B-factors\n    rmsd_map = {resi: np.sqrt(sq_sum \/ n_frames) \n                for resi, sq_sum in sum_sq_dev.items()}\n    \n    stored.mean_rmsd = rmsd_map\n    cmd.alter(mean_obj, \"b=stored.mean_rmsd.get(resi.strip(), 0.0)\")\n    \n    # Visualize with putty representation\n    cmd.show(\"cartoon\", mean_obj)\n    cmd.cartoon(\"putty\", mean_obj)\n    cmd.spectrum(\"b\", \"paleyellow orange pink red black\", mean_obj,\n                minimum=1.0, maximum=4.0)\n    \n    # Putty settings scale tube width by B-factor\n    cmd.set(\"cartoon_putty_transform\", 7, mean_obj)\n    cmd.set(\"cartoon_putty_scale_min\", 0.5, mean_obj)\n    cmd.set(\"cartoon_putty_scale_max\", 3.5, mean_obj)<\/pre>\n\n\n\n
    \"\"<\/a><\/figure>\n\n\n\n
    Mean structures for Folded, PUF1, and PUF2 clusters. Tube width and color indicate per-residue RMSD: thin\/yellow regions are stable; thick\/black regions show high variability.<\/em><\/p>\n\n\n\n
    Ray Tracing Modes for NMR Ensembles<\/h2>\n\n\n\n
    NMR structures with multiple models benefit from careful ray tracing. Mode 3 (quantized color + black outline) works beautifully for showing ensemble spread:<\/p>\n\n\n\n
    # Load NMR ensemble and show all states\ncmd.load(\"2L39_crop.pdb\", \"nmr_ref\")\ncmd.set(\"all_states\", 1, \"nmr_ref\")\n\n# Style settings\ncmd.show(\"cartoon\", \"nmr_ref\")\ncmd.cartoon(\"tube\", \"nmr_ref\")\ncmd.color(\"white\", \"nmr_ref\")\ncmd.set(\"cartoon_tube_radius\", 0.5, \"nmr_ref\")\ncmd.set(\"cartoon_transparency\", 0.3, \"nmr_ref\")\n\n# Ray trace settings for clean ensemble visualization\ncmd.set(\"orthoscopic\", 1)\ncmd.set(\"ray_trace_mode\", 3)  # Quantized + outline\ncmd.set(\"ray_shadows\", 0)\ncmd.set(\"antialias\", 2)\ncmd.bg_color(\"white\")\n\n# Render\ncmd.ray(1200, 1200)<\/pre>\n\n\n\n
    \"\"<\/a><\/figure>\n\n\n\n
    NMR ensemble (PDB: 2L39) rendered with ray_trace_mode 3. The black outlines help distinguish overlapping conformers.<\/em><\/p>\n\n\n\n
    Compare with mode 1 (standard) for a softer look, or experiment with ray_transparency_oblique<\/code> for depth effects:<\/p>\n\n\n\n
    # Alternative: softer rendering with oblique transparency\ncmd.set(\"ray_trace_mode\", 1)\ncmd.set(\"ray_transparency_oblique\", 1.0)\ncmd.set(\"transparency_mode\", 1)<\/pre>\n\n\n\n
    Bringing It Together: A Complete Visualization Pipeline<\/h2>\n\n\n\n
    Here’s how these pieces combine for a weighted ensemble analysis:<\/p>\n\n\n\n
    # 1. Set up environment\ncmd.bg_color(\"white\")\ncmd.set(\"orthoscopic\", 1)\n\n# 2. Load and style reference(s)\ntarget_coords = get_reference_target_coords(\"reference.pdb\", ALIGN_SELECTION)\n\ncmd.load(\"reference.pdb\", \"ref\")\ncmd.cartoon(\"dash\", \"ref\")\ncmd.set(\"cartoon_transparency\", 0.4, \"ref\")\n\n# 3. Load ensemble and trajectory\ncmd.load(\"topology.pdb\", \"ensemble\")\ncmd.load_traj(\"trajectory.xtc\", \"ensemble\")\n\n# 4. Align to reference\nalign_trajectory_to_coords(\"ensemble\", target_coords, ALIGN_SELECTION)\n\n# 5. Load weights and show top structures\nweights = np.load(\"weights.npz\")[\"weights\"]\nshow_top_structures_by_rmsd(\"ensemble\", weights, n=20, reference_obj=\"ref\")\n\n# 6. Final rendering settings\ncmd.set(\"ray_trace_mode\", 3)\ncmd.set(\"ray_shadows\", 0)\ncmd.set(\"antialias\", 2)\ncmd.zoom()<\/pre>\n\n\n\n
    \"\"<\/a><\/figure>\n\n\n\n
    Comparison of ensemble fitting methods for MoPrP. Top structures colored by weight (cyan=low, magenta\/dark=high, overlaid on the NMR reference (gray).<\/em><\/p>\n\n\n\n
    Summary<\/h2>\n\n\n\n
    Visualizing weighted structural ensembles requires moving beyond single-structure thinking:<\/p>\n\n\n\n
      \n
    • Use cartoon dash<\/strong> to distinguish references from ensemble members<\/li>\n\n\n\n
    • Align to frame-averaged coordinates<\/strong> when your reference is itself an ensemble<\/li>\n\n\n\n
    • Color by RMSD<\/strong> to reveal conformational preferences in top structures<\/li>\n\n\n\n
    • Show cluster means with putty representation<\/strong> to communicate both average structure and variability<\/li>\n\n\n\n
    • Experiment with ray trace modes<\/strong> for clean multi-model rendering<\/li>\n<\/ul>\n\n\n\n
      In Part 2, we’ll tackle generating SASA surfaces from weighted ensembles\u2014including efficient algorithms for handling the computational demands of large trajectory datasets.<\/p>\n\n\n\n
      <\/p>\n","protected":false},"excerpt":{"rendered":"
      When working with structural ensembles from molecular dynamics, AlphaFold2 subsampling, or ensemble reweighting against experimental data, you quickly run into visualization problems. Many of these problems standard PyMOL tutorials don’t address: what do you do when there’s no single reference structure? In this two-part series, I’ll share the PyMOL techniques I’ve developed for visualizing weighted […]<\/p>\n","protected":false},"author":115,"featured_media":0,"comment_status":"closed","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"nf_dc_page":"","wikipediapreview_detectlinks":true,"_monsterinsights_skip_tracking":false,"_monsterinsights_sitenote_active":false,"_monsterinsights_sitenote_note":"","_monsterinsights_sitenote_category":0,"ngg_post_thumbnail":0,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[647,351,221],"tags":[132],"ppma_author":[725],"class_list":["post-13926","post","type-post","status-publish","format-standard","hentry","category-molecular-dynamics","category-molecular-visualization","category-python","tag-pymol"],"jetpack_featured_media_url":"","jetpack_sharing_enabled":true,"authors":[{"term_id":725,"user_id":115,"is_guest":0,"slug":"alexi","display_name":"Alexi Hussain","avatar_url":"https:\/\/secure.gravatar.com\/avatar\/aa505acb767aa3aa80da3b68d4684bf5d5064ac6aed2a01206a61383a967535a?s=96&d=mm&r=g","0":null,"1":"","2":"","3":"","4":"","5":"","6":"","7":"","8":""}],"_links":{"self":[{"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/posts\/13926","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/users\/115"}],"replies":[{"embeddable":true,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/comments?post=13926"}],"version-history":[{"count":5,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/posts\/13926\/revisions"}],"predecessor-version":[{"id":13956,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/posts\/13926\/revisions\/13956"}],"wp:attachment":[{"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/media?parent=13926"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/categories?post=13926"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/tags?post=13926"},{"taxonomy":"author","embeddable":true,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/ppma_author?post=13926"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}