Skills › Knowledge Work Plugins › scvi-tools › validate_adata.py
validate_adata.py
397 lines · bio-research/skills/scvi-tools/scripts/validate_adata.py at 2574e29
1#!/usr/bin/env python3
2"""
3Validation utilities for checking AnnData compatibility with scvi-tools.
5Usage:
6 python validate_adata.py data.h5ad
8 # Or import as module
9 from validate_adata import validate_for_scvi, ValidationResult
10 result = validate_for_scvi(adata)
11 print(result.summary())
12"""
14import argparse
15import sys
16from dataclasses import dataclass, field
17from typing import Optional, List, Dict, Any
18import warnings
21@dataclass
22class ValidationResult:
23 """Results from AnnData validation."""
25 is_valid: bool = True
26 errors: List[str] = field(default_factory=list)
27 warnings: List[str] = field(default_factory=list)
28 info: Dict[str, Any] = field(default_factory=dict)
29 recommendations: List[str] = field(default_factory=list)
31 def add_error(self, msg: str):
32 """Add an error (makes validation fail)."""
33 self.errors.append(msg)
34 self.is_valid = False
36 def add_warning(self, msg: str):
37 """Add a warning (doesn't fail validation)."""
38 self.warnings.append(msg)
40 def add_recommendation(self, msg: str):
41 """Add a recommendation for improvement."""
42 self.recommendations.append(msg)
44 def summary(self) -> str:
45 """Generate summary report."""
46 lines = []
47 lines.append("=" * 60)
48 lines.append("scvi-tools AnnData Validation Report")
49 lines.append("=" * 60)
51 # Status
52 status = "PASSED" if self.is_valid else "FAILED"
53 lines.append(f"\nStatus: {status}")
55 # Info
56 if self.info:
57 lines.append("\n--- Data Summary ---")
58 for key, value in self.info.items():
59 lines.append(f" {key}: {value}")
61 # Errors
62 if self.errors:
63 lines.append(f"\n--- Errors ({len(self.errors)}) ---")
64 for i, err in enumerate(self.errors, 1):
65 lines.append(f" {i}. {err}")
67 # Warnings
68 if self.warnings:
69 lines.append(f"\n--- Warnings ({len(self.warnings)}) ---")
70 for i, warn in enumerate(self.warnings, 1):
71 lines.append(f" {i}. {warn}")
73 # Recommendations
74 if self.recommendations:
75 lines.append(f"\n--- Recommendations ({len(self.recommendations)}) ---")
76 for i, rec in enumerate(self.recommendations, 1):
77 lines.append(f" {i}. {rec}")
79 lines.append("\n" + "=" * 60)
80 return "\n".join(lines)
83def validate_for_scvi(
84 adata,
85 layer: Optional[str] = None,
86 batch_key: Optional[str] = None,
87 labels_key: Optional[str] = None,
88 check_hvg: bool = True
89) -> ValidationResult:
90 """
91 Validate AnnData for scvi-tools compatibility.
93 Parameters
94 ----------
95 adata : AnnData
96 Data to validate
97 layer : str, optional
98 Layer containing counts (if None, checks X)
99 batch_key : str, optional
100 Expected batch column in obs
101 labels_key : str, optional
102 Expected labels column in obs
103 check_hvg : bool
104 Check for highly variable genes
106 Returns
107 -------
108 ValidationResult with errors, warnings, and recommendations
109 """
110 import numpy as np
111 from scipy.sparse import issparse
113 result = ValidationResult()
115 # Basic info
116 result.info["shape"] = f"{adata.n_obs} cells x {adata.n_vars} genes"
117 result.info["layers"] = list(adata.layers.keys()) if adata.layers else "None"
119 # Get data matrix to check
120 if layer is not None:
121 if layer not in adata.layers:
122 result.add_error(f"Layer '{layer}' not found. Available: {list(adata.layers.keys())}")
123 return result
124 X = adata.layers[layer]
125 result.info["checking"] = f"layer '{layer}'"
126 else:
127 X = adata.X
128 result.info["checking"] = "adata.X"
130 # Check for None or empty
131 if X is None:
132 result.add_error("Data matrix is None")
133 return result
135 if X.shape[0] == 0 or X.shape[1] == 0:
136 result.add_error(f"Data matrix is empty: shape {X.shape}")
137 return result
139 # Convert to array for checking
140 if issparse(X):
141 result.info["sparse"] = True
142 X_check = X.data # Just check non-zero values
143 else:
144 result.info["sparse"] = False
145 X_check = X.flatten()
147 # Check for raw counts (integers)
148 if len(X_check) > 0:
149 is_integer = np.allclose(X_check, X_check.astype(int))
150 result.info["contains_integers"] = is_integer
152 if not is_integer:
153 result.add_error(
154 "Data does not contain integers (raw counts required). "
155 "Found float values - data may be normalized."
156 )
157 result.add_recommendation(
158 "Use adata.raw.to_adata() to recover raw counts, "
159 "or specify a layer with raw counts"
160 )
162 # Check for negative values
163 min_val = X.min()
164 if min_val < 0:
165 result.add_error(f"Data contains negative values (min={min_val})")
167 # Check for NaN/Inf
168 if issparse(X):
169 has_nan = np.isnan(X.data).any()
170 has_inf = np.isinf(X.data).any()
171 else:
172 has_nan = np.isnan(X).any()
173 has_inf = np.isinf(X).any()
175 if has_nan:
176 result.add_error("Data contains NaN values")
177 if has_inf:
178 result.add_error("Data contains Inf values")
180 # Check data range
181 max_val = X.max()
182 result.info["value_range"] = f"[{min_val}, {max_val}]"
184 if max_val < 10:
185 result.add_warning(
186 f"Maximum value is {max_val}, which is very low. "
187 "Data may be log-transformed or normalized."
188 )
190 # Check sparsity
191 if issparse(X):
192 sparsity = 1 - (X.nnz / (X.shape[0] * X.shape[1]))
193 result.info["sparsity"] = f"{sparsity:.1%}"
195 if sparsity < 0.5:
196 result.add_warning(
197 f"Data is only {sparsity:.1%} sparse. "
198 "Consider if this is expected for your data type."
199 )
201 # Check batch key
202 if batch_key is not None:
203 if batch_key not in adata.obs.columns:
204 result.add_error(
205 f"batch_key '{batch_key}' not found in obs. "
206 f"Available columns: {list(adata.obs.columns)}"
207 )
208 else:
209 n_batches = adata.obs[batch_key].nunique()
210 result.info["n_batches"] = n_batches
212 if n_batches == 1:
213 result.add_warning(
214 "Only 1 batch found. Batch correction may not be needed."
215 )
217 # Check for small batches
218 batch_counts = adata.obs[batch_key].value_counts()
219 small_batches = batch_counts[batch_counts < 50]
220 if len(small_batches) > 0:
221 result.add_warning(
222 f"{len(small_batches)} batches have fewer than 50 cells. "
223 "Consider merging small batches."
224 )
226 # Check labels key
227 if labels_key is not None:
228 if labels_key not in adata.obs.columns:
229 result.add_error(
230 f"labels_key '{labels_key}' not found in obs. "
231 f"Available columns: {list(adata.obs.columns)}"
232 )
233 else:
234 n_labels = adata.obs[labels_key].nunique()
235 result.info["n_labels"] = n_labels
237 # Check for rare labels
238 label_counts = adata.obs[labels_key].value_counts()
239 rare_labels = label_counts[label_counts < 30]
240 if len(rare_labels) > 0:
241 result.add_warning(
242 f"{len(rare_labels)} cell types have fewer than 30 cells. "
243 "Rare types may not be well learned."
244 )
246 # Check HVG
247 if check_hvg:
248 if 'highly_variable' not in adata.var.columns:
249 result.add_recommendation(
250 "No highly variable genes found. Run sc.pp.highly_variable_genes() "
251 "and subset to HVGs for better performance."
252 )
253 else:
254 n_hvg = adata.var['highly_variable'].sum()
255 result.info["n_hvg"] = n_hvg
257 if n_hvg < 1000:
258 result.add_warning(
259 f"Only {n_hvg} HVGs selected. Consider using 2000-4000 for best results."
260 )
261 elif n_hvg > 5000:
262 result.add_warning(
263 f"{n_hvg} HVGs selected. Consider reducing to 2000-4000 "
264 "for efficiency."
265 )
267 # Check gene count
268 if adata.n_vars > 30000:
269 result.add_recommendation(
270 f"Dataset has {adata.n_vars} genes. Subset to HVGs (2000-4000) "
271 "for faster training and better results."
272 )
274 # Check cell count
275 if adata.n_obs < 1000:
276 result.add_warning(
277 f"Dataset has only {adata.n_obs} cells. "
278 "Deep learning models work best with >5000 cells."
279 )
281 # Check for counts layer
282 if layer is None and 'counts' not in adata.layers:
283 result.add_recommendation(
284 "Store raw counts in adata.layers['counts'] before any normalization. "
285 "This preserves the original data for scvi-tools."
286 )
288 # Check for raw attribute
289 if adata.raw is not None:
290 result.info["has_raw"] = True
291 result.add_recommendation(
292 "adata.raw exists. If X is normalized, use adata.raw.to_adata() "
293 "to recover raw counts."
294 )
295 else:
296 result.info["has_raw"] = False
298 return result
301def suggest_model(adata, result: ValidationResult) -> str:
302 """
303 Suggest appropriate scvi-tools model based on data.
305 Parameters
306 ----------
307 adata : AnnData
308 Data to analyze
309 result : ValidationResult
310 Validation result with info
312 Returns
313 -------
314 String with model suggestion
315 """
316 suggestions = []
318 # Check for multi-modal data
319 if 'protein_expression' in adata.obsm:
320 suggestions.append("totalVI: CITE-seq data detected (protein + RNA)")
322 if 'spliced' in adata.layers and 'unspliced' in adata.layers:
323 suggestions.append("veloVI: RNA velocity data detected (spliced + unspliced)")
325 # Check for labels
326 has_labels = result.info.get('n_labels', 0) > 0
327 has_batches = result.info.get('n_batches', 0) > 1
329 if has_batches:
330 if has_labels:
331 suggestions.append(
332 "scANVI: Integration with cell type labels (recommended for label transfer)"
333 )
334 else:
335 suggestions.append(
336 "scVI: Unsupervised batch integration"
337 )
338 else:
339 suggestions.append(
340 "scVI: Dimensionality reduction and differential expression"
341 )
343 if not suggestions:
344 suggestions.append("scVI: General-purpose single-cell analysis")
346 return "\n".join([f" - {s}" for s in suggestions])
349def main():
350 """Command-line interface."""
351 parser = argparse.ArgumentParser(
352 description="Validate AnnData for scvi-tools compatibility"
353 )
354 parser.add_argument("file", help="Path to h5ad file")
355 parser.add_argument("--layer", help="Layer to check (default: X)")
356 parser.add_argument("--batch-key", help="Batch column to check")
357 parser.add_argument("--labels-key", help="Labels column to check")
358 parser.add_argument("--suggest", action="store_true", help="Suggest model type")
360 args = parser.parse_args()
362 try:
363 import scanpy as sc
364 except ImportError:
365 print("Error: scanpy is required. Install with: pip install scanpy")
366 sys.exit(1)
368 # Load data
369 print(f"Loading {args.file}...")
370 try:
371 adata = sc.read_h5ad(args.file)
372 except Exception as e:
373 print(f"Error loading file: {e}")
374 sys.exit(1)
376 # Validate
377 result = validate_for_scvi(
378 adata,
379 layer=args.layer,
380 batch_key=args.batch_key,
381 labels_key=args.labels_key
382 )
384 # Print report
385 print(result.summary())
387 # Suggest model
388 if args.suggest:
389 print("\nSuggested models:")
390 print(suggest_model(adata, result))
392 # Exit code
393 sys.exit(0 if result.is_valid else 1)
396if __name__ == "__main__":
397 main()