-
Notifications
You must be signed in to change notification settings - Fork 5
/
scVelo.py
82 lines (63 loc) · 2.45 KB
/
scVelo.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
###### load python libraries
import sys
import scvelo as scv
import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib as mpl
from matplotlib import cm
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
from matplotlib import pyplot as plt
from matplotlib import rcParams
##### Setup scVelo
sc.settings.verbosity = 3
sc.set_figure_params(dpi=80, color_map='viridis')
scv.settings.set_figure_params('scvelo')
##### Load Data
LoomFile=sys.argv[1]
OUTFILE_PREFIX=sys.argv[2]
Experiment = OUTFILE_PREFIX
adata = scv.read(LoomFile, cache=True)
adata.var_names_make_unique()
cellnames = adata.obs_names
df = pd.DataFrame(adata.obs_names)
df.to_csv('Samples.txt', index=False)
anno = pd.read_csv('Annotations.txt') ###################### Example of Annotations.txt is available on GitHub
adata.obs = anno
#### Basic Filtering
sc.pp.filter_cells(adata, min_genes=0)
sc.pp.filter_genes(adata, min_cells=0)
##### Normalization
sc.pp.normalize_per_cell(adata, counts_per_cell_after=1e4)
sc.pp.log1p(adata)
adata.raw = adata
####### PCA, Clustering and UMAP
########## set n_neighbors, n_pcs and n_comps based on the number of samples#######
sc.tl.pca(adata, svd_solver='arpack', n_comps=3)
sc.pp.neighbors(adata, n_neighbors=2, n_pcs=3)
sc.tl.umap(adata)
sc.tl.louvain(adata)
FigureOut = Experiment + "_UMAP_scVelo.png"
rcParams['figure.figsize'] = 15, 15
sc.pl.umap(adata, color=['louvain'], save=FigureOut)
FigureOut = Experiment + "_UMAP_SampleType_scVelo.png"
rcParams['figure.figsize'] = 15, 15
sc.pl.umap(adata, color=['SampleType'], save=FigureOut)
###### scVelo Velocity Estimation
scv.pp.moments(adata, n_pcs=3)
scv.tl.velocity(adata, mode='stochastic')
scv.tl.velocity_graph(adata)
scv.tl.velocity_embedding(adata, basis='umap')
##### scVelo Visualization
FigureOut = Experiment + "_VelocityEmbedding_scVelo.png"
rcParams['figure.figsize'] = 15, 15
scv.pl.velocity_embedding(adata, basis='umap', save=FigureOut)
FigureOut = Experiment + "_VelocityStream_scVelo.png"
rcParams['figure.figsize'] = 15, 15
scv.pl.velocity_embedding_stream(adata, basis='umap', save=FigureOut)
#### Variance Velocity of Genes of Interest
FigureOut = Experiment + "_Cited1-Wnt4-Hnf4a_Velocity_scVelo.png"
scv.pl.velocity(adata, var_names=['Notch2', 'Wnt4', 'Hnf4a']) ###### Change genes based on your interest
#### writing out scVelo object
results_file = Experiment + 'scVeloAnalysis.h5ad'
adata.write(results_file, compression='gzip')