# Clinical Modeling Workflows > Classify variants via ACMG/AMP + CADD/REVEL/SpliceAI; dock ligands with AutoDock Vina; set up GROMACS MD; run Scanpy scRNA-seq QC/clustering. Use for variant classification, docking, MD setup, or scRNA-seq. - Skill: `pavel-kravchenko/clinical-modeling-workflows` (Agent Skill) - Install (CLI): `npx skillmds@latest add pavel-kravchenko/clinical-modeling-workflows` - Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/clinical-modeling-workflows/raw - Safety review: pending - Works with: Claude Code, Claude.ai, OpenAI Codex - Category: Coding & Dev Tools - Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko) - Updated: 2026-09-17 - Page: https://skillmd.com/skills/pavel-kravchenko/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)` | Strength | Pathogenic codes | Key triggers | |----------|-----------------|-------------| | Very Strong | PVS1 | Null variant (nonsense/frameshift/splice) in a LoF-mechanism disease gene | | Strong | PS1-PS4 | Same AA change as established pathogenic; confirmed de novo; damaging functional assay; prevalence in affected >> controls | | Moderate | PM1-PM6 | Mutational hotspot; absent from population databases; in-frame indel; assumed de novo; novel missense at a known pathogenic residue | | Supporting | PP1-PP5 | Co-segregation with disease; PP3 computational evidence; phenotype specific for the gene | | Stand-alone benign | BA1 | MAF > 5% in any general population | | Strong benign | BS1-BS4 | Frequency greater than expected for disorder; observed in healthy adult; benign functional study; lack of segregation | | Supporting benign | BP1-BP7 | Missense in gene where only truncating cause disease; in silico benign (BP4); synonymous with no predicted splice effect (BP7) | ### In Silico Predictors (for PP3 / BP4) | Tool | Range | Damaging threshold | |------|-------|--------------------| | SIFT | 0-1 (lower = damaging) | < 0.05 | | PolyPhen-2 | 0-1 (higher = damaging) | > 0.908 probably damaging; > 0.446 possibly damaging | | CADD (Phred) | higher = worse | ≥ 20 (top 1%); ≥ 25 (top 0.3%) | | REVEL | 0-1 | > 0.75 likely pathogenic | | AlphaMissense | 0-1 | > 0.564 likely pathogenic; < 0.34 likely benign | | SpliceAI (delta score) | 0-1 | > 0.2 suggestive; > 0.5 high confidence | 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). ```python 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) ```bash 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. ```python 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 ```bash 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. ```python 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