| name | harmonization-qc |
| description | Visual quality control for ComBat-family harmonization results. Generate PCA/UMAP plots colored by site before and after harmonization, heatmap comparisons, variance explained by site, and a summary report. Invoke when the user asks for "harmonization QC", "check harmonization", "visualize ComBat results", "PCA before after harmonization", "evaluate site correction", "does harmonization work", or "validate batch correction". Also invoke after any harmonization completes to verify the results are meaningful.
|
Harmonization Quality Control
Generate visual and quantitative assessments of how well harmonization
removed site effects while preserving biological signal.
Workflow
1. Load data
Read the pre- and post-harmonization feature matrices and metadata:
import numpy as np
import pandas as pd
orig = np.load('harmonization/input_edges.npy')
harm = np.load('harmonization/harmonized_edges.npy')
meta = pd.read_csv('harmonization/subject_metadata.csv')
sites = meta['site'].values
2. Quantitative metrics
Compute these metrics for both pre- and post-harmonization:
Site variance explained (lower = better):
from sklearn.decomposition import PCA
from scipy.stats import f_oneway
pca = PCA(n_components=min(10, n_subjects-1))
pca.fit(features)
for i in range(pca.n_components_):
f_stat, p_val = f_oneway(*[pca.transform(features)[sites==s, i] for s in np.unique(sites)])
print(f"PC{i+1}: F={f_stat:.2f}, p={p_val:.4f}")
Silhouette score by site (lower = better, sites should intermingle):
from sklearn.metrics import silhouette_score
score = silhouette_score(features, sites)
Correlation of site with each edge (count edges with |r| > 0.3 as "site-biased"):
site_encoded = pd.get_dummies(sites)
site_corrs = np.array([np.abs(np.corrcoef(site_encoded[s], features[:, i])[0,1])
for i in range(features.shape[1])])
n_biased = np.sum(site_corrs > 0.3)
3. Generate figures
Create a multi-panel QC figure (harmonization_qc.png):
Panel 1: PCA scatter (top row)
- Left: Pre-harmonization PC1 vs PC2, colored by site
- Right: Post-harmonization PC1 vs PC2, colored by site
- Sites should visually intermingle after harmonization
Panel 2: UMAP embedding (optional, middle row)
- Same layout as PCA but using UMAP for nonlinear structure
Panel 3: Edge-wise variance comparison (bottom left)
- Scatter plot: pre-harmonization SD per edge vs post-harmonization SD per edge
- Diagonal line for reference
- Points below the line = variance reduced
Panel 4: Bar chart of site variance explained by PC (bottom right)
- Grouped bars: pre vs post for top 5 PCs
- Y-axis: proportion of variance explained by site (from ANOVA F)
4. Generate report
Produce a text summary:
QC Report: ComBat Harmonization
================================
Method: ComBat-GAM
Date: YYYY-MM-DD
N subjects: XX
N sites: XX
N features: XX
Site variance explained:
Before: PC1=F=XX.X (p<.001), PC2=F=XX.X (p<.001)
After: PC1=F=X.X (p>.05), PC2=F=X.X (p>.05)
Site-biased edges:
Before: XX edges (XX%)
After: XX edges (XX%)
Silhouette by site:
Before: X.XXX
After: X.XXX (lower is better)
Assessment: PASS/FLAG
5. Writing the QC script
Write analysis_<N>_harmonization_qc.py with:
- All imports (numpy, pandas, matplotlib, sklearn, scipy)
- Figure generation (mpl.figure with subplots)
- Text report output to stdout and to
harmonization_qc_report.txt
- Save figure to
harmonization_qc.png at 300 DPI
6. Assessment criteria
PASS (harmonization likely succeeded):
- PCA site separation visibly reduced
- Silhouette score decreased
- Site-biased edge count decreased
- No NaN/Inf values in harmonized data
FLAG (needs investigation):
- PCA shows new/unexpected clusters after harmonization
- Many edges show increased variance
- Site-biased edges increased
- NaN values present
Constraints
- Always use both quantitative and visual assessment
- Never pass/fail based on visualization alone — compute metrics
- If UMAP is not installed (
pip install umap-learn), skip that panel rather than failing
- Set matplotlib style for publication-quality figures
