|
| 1 | +import os |
| 2 | +import random |
| 3 | +import warnings |
| 4 | + |
| 5 | +# stlearn imports tensorflow, which logs its device setup to stderr on import |
| 6 | +os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' |
| 7 | + |
| 8 | +import anndata as ad |
| 9 | +import numpy as np |
| 10 | +import pandas as pd |
| 11 | +import scanpy as sc |
| 12 | +import scipy.sparse as sp |
| 13 | +import stlearn as st |
| 14 | + |
| 15 | +## VIASH START |
| 16 | +par = { |
| 17 | + 'input': 'resources_test/task_spatial_trajectory_inference/dlpfc_151673/dataset.h5ad', |
| 18 | + 'output': 'output.h5ad', |
| 19 | + 'n_comps': 50, |
| 20 | + 'n_neighbors': 200, |
| 21 | + 'resolution': 0.8, |
| 22 | + 'eps': 1500.0, |
| 23 | + 'seed': 0, |
| 24 | +} |
| 25 | +meta = { |
| 26 | + 'name': 'stlearn' |
| 27 | +} |
| 28 | +## VIASH END |
| 29 | + |
| 30 | +# warnings raised by stlearn internals, not actionable from here |
| 31 | +warnings.simplefilter('ignore', FutureWarning) |
| 32 | +warnings.simplefilter('ignore', ad.ImplicitModificationWarning) |
| 33 | + |
| 34 | +seed = par['seed'] |
| 35 | +np.random.seed(seed) |
| 36 | +random.seed(seed) |
| 37 | + |
| 38 | + |
| 39 | +def build_cell_type_int(adata): |
| 40 | + """Map cell_type strings to integer labels, as stLearn expects numeric cluster labels.""" |
| 41 | + adata.obs['cell_type'] = adata.obs['cell_type'].astype(str).astype('category') |
| 42 | + unique_cell_types = adata.obs['cell_type'].cat.categories |
| 43 | + |
| 44 | + ct_to_num = {str(ct): str(i) for i, ct in enumerate(unique_cell_types)} |
| 45 | + num_to_ct = {str(i): str(ct) for i, ct in enumerate(unique_cell_types)} |
| 46 | + |
| 47 | + adata.obs['cell_type_int'] = ( |
| 48 | + adata.obs['cell_type'].astype(str).map(ct_to_num).astype('category') |
| 49 | + ) |
| 50 | + return num_to_ct |
| 51 | + |
| 52 | + |
| 53 | +def select_root(adata): |
| 54 | + """Return the least differentiated cluster, i.e. the one whose cells express the largest |
| 55 | + number of genes on average. This is the same proxy stLearn uses to pick the root cell.""" |
| 56 | + n_expressed = np.asarray((adata.layers['counts'] > 0).sum(axis=1)).reshape(-1) |
| 57 | + scores = ( |
| 58 | + pd.DataFrame({'n_expressed': n_expressed, 'cluster': adata.obs['cell_type_int'].values}) |
| 59 | + .groupby('cluster', observed=True)['n_expressed'] |
| 60 | + .mean() |
| 61 | + ) |
| 62 | + return str(scores.idxmax()) |
| 63 | + |
| 64 | + |
| 65 | +def filter_branches(available_paths, root): |
| 66 | + """Keep the paths starting from the root and drop branches contained in a longer one.""" |
| 67 | + valid = [path for path in available_paths.values() if path[0] == root] |
| 68 | + valid.sort(key=len, reverse=True) |
| 69 | + |
| 70 | + unique = [] |
| 71 | + for branch in valid: |
| 72 | + if not any(set(branch).issubset(set(kept)) for kept in unique): |
| 73 | + unique.append(branch) |
| 74 | + return unique |
| 75 | + |
| 76 | + |
| 77 | +print('Reading input files', flush=True) |
| 78 | +adata = ad.read_h5ad(par['input']) |
| 79 | + |
| 80 | +adata.X = adata.layers['normalized'] |
| 81 | +if sp.issparse(adata.X): |
| 82 | + adata.X = adata.X.toarray() |
| 83 | + |
| 84 | +# stLearn reads the spatial coordinates from these slots |
| 85 | +adata.obsm['spatial'] = adata.obsm['X_spatial'] |
| 86 | +adata.obs['imagerow'] = adata.obsm['X_spatial'][:, 1] |
| 87 | +adata.obs['imagecol'] = adata.obsm['X_spatial'][:, 0] |
| 88 | + |
| 89 | +print('Embed and cluster the cells', flush=True) |
| 90 | +st.em.run_pca(adata, n_comps=par['n_comps']) |
| 91 | +sc.pp.neighbors(adata, n_neighbors=par['n_neighbors'], use_rep='X_pca') |
| 92 | +st.tl.clustering.leiden(adata, resolution=par['resolution'], random_state=seed) |
| 93 | + |
| 94 | +num_to_ct = build_cell_type_int(adata) |
| 95 | + |
| 96 | +print('Determine root cell', flush=True) |
| 97 | +root = select_root(adata) |
| 98 | +print(f'Root cluster: {num_to_ct[root]}', flush=True) |
| 99 | + |
| 100 | +# use_raw is False because the dataset does not carry a raw layer |
| 101 | +adata.uns['iroot'] = st.spatial.trajectory.set_root( |
| 102 | + adata, |
| 103 | + use_label='cell_type_int', |
| 104 | + cluster=int(root), |
| 105 | + use_raw=False, |
| 106 | +) |
| 107 | + |
| 108 | +print('Calculate the pseudotime and the trajectory branches', flush=True) |
| 109 | +st.spatial.trajectory.pseudotime( |
| 110 | + adata, |
| 111 | + eps=par['eps'], |
| 112 | + use_rep='X_pca', |
| 113 | + use_label='cell_type_int', |
| 114 | +) |
| 115 | +branches = filter_branches(adata.uns.get('available_paths', {}), int(root)) |
| 116 | + |
| 117 | +if not branches: |
| 118 | + raise RuntimeError( |
| 119 | + f'No valid branches found from root cluster {root}. Try adjusting --resolution or --eps.' |
| 120 | + ) |
| 121 | + |
| 122 | +print('Generate predictions', flush=True) |
| 123 | +# cells not covered by any branch keep a NaN pseudotime |
| 124 | +adata.obs['pseudotime_inferred'] = np.nan |
| 125 | + |
| 126 | +for branch in branches: |
| 127 | + branch_labels = [str(node) for node in branch] |
| 128 | + try: |
| 129 | + st.spatial.trajectory.pseudotimespace_global( |
| 130 | + adata, |
| 131 | + use_label='cell_type_int', |
| 132 | + list_clusters=branch_labels, |
| 133 | + ) |
| 134 | + except Exception as exc: |
| 135 | + print(f'Skipping branch {branch}: {exc}', flush=True) |
| 136 | + continue |
| 137 | + |
| 138 | + # the first branch a cell belongs to takes precedence |
| 139 | + unfilled = ( |
| 140 | + adata.obs['cell_type_int'].isin(branch_labels) |
| 141 | + & adata.obs['pseudotime_inferred'].isna() |
| 142 | + ) |
| 143 | + adata.obs.loc[unfilled, 'pseudotime_inferred'] = adata.obs.loc[unfilled, 'dpt_pseudotime'] |
| 144 | + |
| 145 | +n_assigned = adata.obs['pseudotime_inferred'].notna().sum() |
| 146 | +print(f'Pseudotime assigned to {n_assigned}/{adata.n_obs} cells', flush=True) |
| 147 | + |
| 148 | +print('Write output AnnData to file', flush=True) |
| 149 | +output = ad.AnnData( |
| 150 | + obs=adata.obs[['pseudotime_inferred']], |
| 151 | + uns={ |
| 152 | + 'dataset_id': adata.uns['dataset_id'], |
| 153 | + 'normalization_id': adata.uns['normalization_id'], |
| 154 | + 'method_id': meta['name'], |
| 155 | + }, |
| 156 | +) |
| 157 | +output.write_h5ad(par['output'], compression='gzip') |
0 commit comments