Setting the file. One moment. Spatial Deconvolution
22
Validate Data
Tech Debt
62
Recruiting Pipeline
71
Vendor Check
125
Zoom Meeting SDK Web
88
Vendor Review
181
Create An Asset
Video Sdk/web
scripts/prepare_data.py
scripts/prepare_data.py
Python·169 lines·5 KB
def
prepare_data
(
18 adata,
19 batch_key=None,
20 n_top_genes=2000,
21 min_genes=200,
22 max_genes=5000,
23 max_mito_pct=20.0,
24 min_cells=3,
25 skip_filter=False
26):
27 """
28 Prepare AnnData for scvi-tools.
29
30 Parameters
31 ----------
32 adata : AnnData
33 Raw count data
34 batch_key : str, optional
35 Batch column for batch-aware HVG selection
36 n_top_genes : int
37 Number of highly variable genes
38 min_genes : int
39 Minimum genes per cell
40 max_genes : int
41 Maximum genes per cell
42 max_mito_pct : float
43 Maximum mitochondrial percentage
44 min_cells : int
45 Minimum cells per gene
46 skip_filter : bool
47 Skip QC filtering (use if already filtered)
48
49 Returns
50 -------
51 AnnData prepared for scvi-tools
52 """
53 import scanpy as sc
54 import numpy as np
55 from model_utils import get_mito_genes
56
57 adata = adata.copy()
58 print(f"Input: {adata.shape[0]} cells, {adata.shape[1]} genes")
59
60 if not skip_filter:
61 # Calculate QC metrics
62 adata.var['mt'] = get_mito_genes(adata)
63 sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)
64
65 # Filter cells
66 n_before = adata.n_obs
67 adata = adata[adata.obs['n_genes_by_counts'] >= min_genes].copy()
68 adata = adata[adata.obs['n_genes_by_counts'] <= max_genes].copy()
69 adata = adata[adata.obs['pct_counts_mt'] < max_mito_pct].copy()
70 print(f"Filtered cells: {n_before} → {adata.n_obs}")
71
72 # Filter genes
73 n_genes_before = adata.n_vars
74 sc.pp.filter_genes(adata, min_cells=min_cells)
75 print(f"Filtered genes: {n_genes_before} → {adata.n_vars}")
76
77 # Store raw counts in layer
78 adata.layers["counts"] = adata.X.copy()
79
80 # HVG selection
81 if batch_key is not None and batch_key in adata.obs.columns:
82 print(f"Selecting {n_top_genes} HVGs (batch-aware: {batch_key})")
83 sc.pp.highly_variable_genes(
84 adata,
85 n_top_genes=n_top_genes,
86 flavor="seurat_v3",
87 batch_key=batch_key,
88 layer="counts"
89 )
90 else:
91 print(f"Selecting {n_top_genes} HVGs")
92 # Need to normalize for non-seurat_v3
93 sc.pp.normalize_total(adata, target_sum=1e4)
94 sc.pp.log1p(adata)
95 sc.pp.highly_variable_genes(adata, n_top_genes=n_top_genes)
96 # Restore counts to X
97 adata.X = adata.layers["counts"].copy()
98
99 # Subset to HVGs
100 n_hvg = adata.var['highly_variable'].sum()
101 adata = adata[:, adata.var['highly_variable']].copy()
102 print(f"Selected {n_hvg} highly variable genes")
103
104 print(f"Output: {adata.shape[0]} cells, {adata.shape[1]} genes")
105
106 return adata
107
108
109def main():
110 parser = argparse.ArgumentParser(
111 description="Prepare AnnData for scvi-tools",
112 formatter_class=argparse.RawDescriptionHelpFormatter,
113 epilog="""
114Examples:
115 # Basic preparation
116 python prepare_data.py raw.h5ad prepared.h5ad
117
118 # With batch-aware HVG selection
119 python prepare_data.py raw.h5ad prepared.h5ad --batch-key sample
120
121 # Custom parameters
122 python prepare_data.py raw.h5ad prepared.h5ad --n-hvgs 3000 --max-mito 15
123
124 # Skip filtering (data already QC'd)
125 python prepare_data.py filtered.h5ad prepared.h5ad --no-filter
126 """
127 )
128 parser.add_argument("input", help="Input h5ad file")
129 parser.add_argument("output", help="Output h5ad file")
130 parser.add_argument("--batch-key", help="Batch column for HVG selection")
131 parser.add_argument("--n-hvgs", type=int, default=2000, help="Number of HVGs (default: 2000)")
132 parser.add_argument("--min-genes", type=int, default=200, help="Min genes per cell (default: 200)")
133 parser.add_argument("--max-genes", type=int, default=5000, help="Max genes per cell (default: 5000)")
134 parser.add_argument("--max-mito", type=float, default=20.0, help="Max mito %% (default: 20)")
135 parser.add_argument("--min-cells", type=int, default=3, help="Min cells per gene (default: 3)")
136 parser.add_argument("--no-filter", action="store_true", help="Skip QC filtering")
137
138 args = parser.parse_args()
139
140 try:
141 import scanpy as sc
142 except ImportError:
143 print("Error: scanpy required. Install with: pip install scanpy")
144 sys.exit(1)
145
146 # Load data
147 print(f"Loading {args.input}...")
148 adata = sc.read_h5ad(args.input)
149
150 # Prepare
151 adata = prepare_data(
152 adata,
153 batch_key=args.batch_key,
154 n_top_genes=args.n_hvgs,
155 min_genes=args.min_genes,
156 max_genes=args.max_genes,
157 max_mito_pct=args.max_mito,
158 min_cells=args.min_cells,
159 skip_filter=args.no_filter
160 )
161
162 # Save
163 print(f"Saving to {args.output}...")
164 adata.write_h5ad(args.output)
165 print("Done!")
166
167
168if __name__ == "__main__":
169 main()