Clinical Modeling Workflows

Clinical Modeling Workflows

When to Use

• Classifying germline/somatic variants for clinical reporting under ACMG/AMP guidelines
• Deciding whether an in silico predictor score (CADD, REVEL, SpliceAI, AlphaMissense) counts as PP3/BP4 evidence
• Virtual screening or pose scoring with AutoDock Vina
• Setting up or troubleshooting a GROMACS MD simulation (solvation, equilibration, production)
• Single-cell RNA-seq QC, clustering, and cell-type annotation with Scanpy/AnnData

Version Compatibility

• Python ≥3.10, scanpy ≥1.10, anndata ≥0.10
• AutoDock Vina ≥1.2 (new scoring function; Vina 1.1 syntax differs slightly)
• GROMACS ≥2023
• ACMG/AMP rules per Richards et al. 2015 (Table 5); thresholds below reflect commonly used defaults, not a single canonical cutoff

Prerequisites

• pip install scanpy pandas numpy requests (requests only needed for live ClinVar/NCBI queries)
• AutoDock Vina + MGLTools (prepare_receptor4.py, prepare_ligand4.py) or Meeko for PDBQT prep
• GROMACS built with the force field you intend to use (e.g. AMBER99SB-ILDN)
• Familiarity with VCF/variant annotation (see bio-applied-variant-calling-and-snp-analysis) and PDB structure files (see bio-applied-structural-methods)

ACMG/AMP Variant Classification

Pathogenic (P) > Likely Pathogenic (LP) > VUS > Likely Benign (LB) > Benign (B)

In Silico Predictors (for PP3 / BP4)

Per current ClinGen SVI recommendations, require concordance across multiple orthogonal tools (not a single one) before applying PP3/BP4 — a common convention is ≥4 of 5 tools agreeing.

gnomAD constraint: pLI > 0.9 = haploinsufficient/LoF-intolerant; LOEUF < 0.35 = strong LoF constraint; missense o/e < 0.5 = missense-constrained. ClinVar review status: **** expert panel/practice guideline, *** multiple submitters no conflicts, * single submitter.

Goal: turn a list of ACMG criteria codes into a final classification. Approach: count codes by strength/direction, then apply the Richards 2015 Table 5 combining rules (benign rules checked first since BA1 short-circuits everything else).

from collections import Counter

CRITERIA_STRENGTH = {
    'PVS1': 'very_strong',
    'PS1': 'strong', 'PS2': 'strong', 'PS3': 'strong', 'PS4': 'strong',
    'PM1': 'moderate', 'PM2': 'moderate', 'PM3': 'moderate',
    'PM4': 'moderate', 'PM5': 'moderate', 'PM6': 'moderate',
    'PP1': 'supporting', 'PP2': 'supporting', 'PP3': 'supporting',
    'PP4': 'supporting', 'PP5': 'supporting',
    'BA1': 'stand_alone',
    'BS1': 'strong', 'BS2': 'strong', 'BS3': 'strong', 'BS4': 'strong',
    'BP1': 'supporting', 'BP2': 'supporting', 'BP3': 'supporting',
    'BP4': 'supporting', 'BP5': 'supporting', 'BP6': 'supporting', 'BP7': 'supporting',
}


def classify_variant(criteria: list[str]) -> str:
    """Classify a variant from ACMG/AMP criteria codes (Richards et al. 2015, Table 5).

    Args:
        criteria: e.g. ['PVS1', 'PM2', 'PP3']
    Returns:
        One of 'Pathogenic', 'Likely Pathogenic', 'VUS', 'Likely Benign', 'Benign'.
    """
    path = [c for c in criteria if c.startswith(('PVS', 'PS', 'PM', 'PP'))]
    benign = [c for c in criteria if c.startswith(('BA', 'BS', 'BP'))]
    strength = Counter(CRITERIA_STRENGTH.get(c) for c in path)
    bstrength = Counter(CRITERIA_STRENGTH.get(c) for c in benign)
    pvs, ps, pm, pp = (strength[s] for s in ('very_strong', 'strong', 'moderate', 'supporting'))
    ba, bs, bp = bstrength['stand_alone'], bstrength['strong'], bstrength['supporting']

    if ba >= 1 or bs >= 2:
        return 'Benign'
    if (bs >= 1 and bp >= 1) or bp >= 2:
        return 'Likely Benign'

    if pvs >= 1 and (ps >= 1 or pm >= 2 or (pm >= 1 and pp >= 1) or pp >= 2):
        return 'Pathogenic'
    if ps >= 2:
        return 'Pathogenic'
    if ps >= 1 and (pm >= 3 or (pm >= 2 and pp >= 2) or (pm >= 1 and pp >= 4)):
        return 'Pathogenic'

    if pvs >= 1 and pm >= 1:
        return 'Likely Pathogenic'
    if ps >= 1 and 1 <= pm <= 2:
        return 'Likely Pathogenic'
    if ps >= 1 and pp >= 2:
        return 'Likely Pathogenic'
    if pm >= 3 or (pm >= 2 and pp >= 2) or (pm >= 1 and pp >= 4):
        return 'Likely Pathogenic'
    if pvs >= 1:  # pragmatic convention (e.g. InterVar); not stated verbatim in Table 5
        return 'Likely Pathogenic'

    return 'VUS'


assert classify_variant(['PVS1', 'PM2']) == 'Likely Pathogenic'
assert classify_variant(['PVS1', 'PS1']) == 'Pathogenic'
assert classify_variant(['BA1']) == 'Benign'
assert classify_variant(['PM2', 'PP3']) == 'VUS'

Molecular Docking (AutoDock Vina)

prepare_receptor4.py -r protein.pdb -o receptor.pdbqt
prepare_ligand4.py -l ligand.mol2 -o ligand.pdbqt
vina --receptor receptor.pdbqt --ligand ligand.pdbqt \
     --config config.txt --out output.pdbqt
# config.txt: center_x/y/z, size_x/y/z define the search box (Å)
# Scores in kcal/mol -- more negative = stronger predicted binding

Force field energy: E_total = E_bond + E_angle + E_dihedral + E_electrostatic + E_vdW. Common FFs: AMBER (proteins/nucleic acids), CHARMM (proteins/lipids), OPLS-AA (small organics) — never mix parameters across FFs.

Goal: parse a Vina log/output table into a DataFrame for ranking and plotting. Approach: split the fixed-width text block, coerce columns, then sort by affinity.

import pandas as pd


def parse_vina_output(text: str) -> pd.DataFrame:
    """Parse AutoDock Vina docking output into a DataFrame.

    Args:
        text: raw stdout/log from `vina`, containing the 'mode | affinity | rmsd' table.
    Returns:
        DataFrame with columns mode, affinity_kcal_mol, rmsd_lb, rmsd_ub.
    """
    rows = []
    for line in text.strip().split('\n'):
        line = line.strip()
        if line and line[0].isdigit():
            parts = line.split()
            rows.append({
                'mode': int(parts[0]),
                'affinity_kcal_mol': float(parts[1]),
                'rmsd_lb': float(parts[2]),
                'rmsd_ub': float(parts[3]),
            })
    return pd.DataFrame(rows)


vina_output = """mode |   affinity | dist from best mode
     | (kcal/mol) | rmsd l.b.| rmsd u.b.
-----+------------+----------+----------
   1       -8.7      0.000      0.000
   2       -8.3      1.245      2.187
   3       -7.9      2.056      4.321
"""
hits = parse_vina_output(vina_output)
best = hits.loc[hits['affinity_kcal_mol'].idxmin()]
assert best['mode'] == 1

GROMACS MD Workflow

gmx pdb2gmx -f protein.pdb -o protein.gro -water tip3p -ff amber99sb-ildn
gmx editconf -f protein.gro -o box.gro -c -d 1.0 -bt dodecahedron
gmx solvate -cp box.gro -cs spc216.gro -o solvated.gro -p topol.top
# Then: add ions (gmx genion) -> energy minimization (em.mdp) ->
#       NVT equilibration -> NPT equilibration -> production MD (md.mdp)

Analysis: gmx rms (RMSD), gmx rmsf (per-residue fluctuation), gmx gyrate (radius of gyration), gmx hbond (hydrogen bonds).

Homology model quality: >50% sequence identity = reliable; 30-50% = reasonable; <30% = twilight zone. AlphaFold pLDDT: >90 high confidence; 70-90 moderate; <50 likely disordered — always energy-minimize a model before MD.

Scanpy Single-Cell Pipeline

Goal: go from a raw 10x count matrix to QC-filtered, clustered, annotated cell types. Approach: filter → normalize → find HVGs (saving raw first) → PCA/neighbors/UMAP → Leiden clustering → marker-gene-based annotation.

import scanpy as sc

adata = sc.read_10x_mtx('filtered_feature_bc_matrix/')  # or sc.datasets.pbmc3k()

# QC: flag mitochondrial genes, compute metrics, filter
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata = adata[adata.obs.pct_counts_mt < 20]

# Normalize -> log -> HVG -> snapshot raw -> scale
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata.raw = adata                              # save unscaled data before subsetting
adata = adata[:, adata.var.highly_variable]
sc.pp.scale(adata, max_value=10)

# Dimensionality reduction + clustering
sc.tl.pca(adata, svd_solver='arpack', n_comps=50)
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=40)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5, flavor='igraph', n_iterations=2)  # 0.2 coarse - 2.0 fine

# Marker genes + manual annotation
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
cluster_annotations = {'0': 'CD4+ T cells', '1': 'CD14+ Monocytes'}
adata.obs['cell_type'] = adata.obs['leiden'].map(cluster_annotations)

Pitfalls

• ACMG: BA1 alone is stand-alone benign — no pathogenic evidence overrides it. VUS is the correct default when evidence is insufficient or conflicting.
• Variant predictors: no single tool is authoritative for PP3/BP4; require concordance across multiple orthogonal tools.
• gnomAD frequency: use population-specific max MAF (popmax), not global AF; recessive disorders tolerate higher carrier frequency than dominant ones.
• Force fields: never mix parameters from different FFs; always energy-minimize before starting MD.
• Docking scores: reliable for ranking poses of the same ligand; not reliable for cross-ligand ranking without rescoring (e.g. MM-GBSA).
• Scanpy adata.raw: must be set before HVG subsetting — rank_genes_groups and DE tools read from the raw, unscaled data.
• Leiden resolution: default 1.0 tends to over-split PBMC-like data; 0.3-0.6 typically yields interpretable major cell types. Use flavor='igraph' (current scanpy default going forward) for speed and determinism.

See Also

• bio-applied-clinical-genomics — VCF annotation pipelines and ClinVar/gnomAD lookups feeding ACMG scoring
• bio-applied-variant-calling-and-snp-analysis — upstream variant calling before classification
• bio-applied-docking — deeper AutoDock Vina / docking workflow patterns
• bio-applied-single-cell-scanpy — advanced Scanpy/AnnData QC, integration, and annotation patterns