Hi, Thank you for your reply. the code I use to do the processing is as follow:
I setup the dataset for both harmony and scvi as follow:
import scanpy as sc
#Create adata and basic filtering
adata1 = sc.read('path/adata1.loom', cache=True, var_names='var_names', obs_names='obs_names', validate=False)
adata2 = sc.read('path/adata2.loom', cache=True, var_names='var_names', obs_names='obs_names', validate=False)
adata3 = sc.read('path/adata3.loom', cache=True, var_names='var_names', obs_names='obs_names', validate=False)
adata4 = sc.read('path/adata4.loom', cache=True, var_names='var_names', obs_names='obs_names', validate=False)
adata5 = sc.read('path/adata5.loom', cache=True, var_names='var_names', obs_names='obs_names', validate=False)
#the same with the others 5 samples..
adata1.obs["sample"]='glioma 1'
adata2.obs["sample"]='glioma 2'
adata3.obs["sample"]='glioma 3'
adata4.obs["sample"]='glioma 4'
adata5.obs["sample"]='glioma 5'
#the same with the others 5 samples..
adata = adata1.concatenate(adata2, adata3, adata4, adata5, adata6,
adata7,adata8,adata9,adata10, join='outer')
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
sc.pl.highest_expr_genes(adata, n_top=20, )
mito_genes = adata.var.index.str.startswith('mt-')
adata.obs['percent_mito'] = np.sum(adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1
adata.obs['n_counts'] = adata.X.sum(axis=1).A1
#Filter cells by QC
adata = adata[adata.obs.n_genes < 6000, :]
adata = adata[adata.obs.n_counts < 25000, :]
adata = adata[adata.obs.n_counts > 1000, :]
adata = adata[adata.obs.percent_mito < 0.06, :]
sc.pl.violin(adata, ['n_genes', 'n_counts', 'percent_mito'], jitter=0.4, multi_panel=True)
sc.pl.scatter(adata, x='n_counts', y='percent_mito')
sc.pl.scatter(adata, x='n_counts', y='n_genes')
#Normalize data
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
Then my code for scanpy + harmony is:
#Find variable genes, regress non relevant variants, calculate PCA
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=6, min_disp=0.10)
sc.pl.highly_variable_genes(adata)
sc.pp.regress_out(adata, ['n_counts', 'percent_mito'], n_jobs=56)
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, svd_solver='arpack', n_comps=100)
sc.pl.pca_variance_ratio(adata, log=True, n_pcs = 100)
sc.pp.neighbors(adata, n_neighbors=30, n_pcs=75)
#First view of all data set without integration
sc.tl.umap(adata)
sc.tl.leiden(adata)
sc.pl.umap(adata, color=['leiden'])
#run harmony and generate new UMAP
import harmonypy as hm
adata.obsm['X_pca_old'] = adata.obsm['X_pca']
adata.obsm['X_umap_old'] = adata.obsm['X_umap']
data_mat = adata.obsm['X_pca']
meta_data = adata.obs
vars_use = ['sample']
ho = hm.run_harmony(data_mat, meta_data, vars_use)
adjusted_pcs = pd.DataFrame(ho.Z_corr).T
adata.obsm['X_pca'] = adjusted_pcs.values
sc.pp.neighbors(adata, n_neighbors=30, n_pcs=75)
sc.tl.umap(adata) maxiter=300)
sc.tl.leiden(adata)
sc.pl.umap(adata, color=['leiden'], legend_loc='on data' )
Alternatively my code for scvi to reach this point is:
import scvi
sc.pp.highly_variable_genes(adata,flavor="seurat_v3",n_top_genes=2000,batch_key="sample", subset=True)
scvi.data.setup_anndata(adata batch_key="sample")
vae = scvi.model.SCVI(adata, n_layers=2, n_latent=30)
vae.train()
adata.obsm["X_scVI"] = vae.get_latent_representation()
sc.pp.neighbors(adata, use_rep="X_scVI")
sc.tl.leiden(adata, resolution=1)
sc.tl.umap(adata, maxiter=300)
sc.pl.umap(adata, color=["leiden"], legend_loc = 'on data')
**Finally I run scVelo as follow:**
import scvelo as scv
scv.set_figure_params()
cv.tl.recover_dynamics(adata, n_jobs=56)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding(adata, basis='umap', color='leiden')
scv.pl.velocity_embedding_grid(adata, basis='umap', color='leiden')
scv.pl.velocity_embedding_stream(adata, basis='umap', color='leiden')
scv.tl.latent_time(adata)
scv.pl.scatter(adata, color='latent_time', color_map='gnuplot', size=80)
#Latent Time with harmony
Regions close to 1 are by gene markers differentiated states of the tumors