1.4_process_and_compute_velocity.md

jupyter	jupytext kernelspec formats text_representation ipynb,md extension format_name format_version jupytext_version .md markdown 1.1 1.2.4 display_name language name (Fertig) Python 3.7 python fertig_python_3_7

from IPython.core.display import display, HTML
display(HTML("<style>.container { width:90% !important; }</style>"))
%matplotlib inline

import scanpy as sc
import scvelo as scv
import scipy as scp
import pandas as pd
import anndata as ad
import numpy as np
scv.logging.print_version()

scv.settings.verbosity = 3  # show errors(0), warnings(1), info(2), hints(3)
scv.settings.set_figure_params('scvelo')  # for beautified visualization

adata = scv.read('finalConcatanatedData.h5ad')
scv.utils.show_proportions(adata)
adata

# remove duplicated values that have different conditions
phenoDat = pd.read_csv('phenoSCCall.csv')
phenoDat.drop('Unnamed: 0',axis=1, inplace=True)
print(phenoDat.shape)
dup = phenoDat.duplicated('barcodes',keep=False)
dupPhenoDat = phenoDat.loc[dup,:]
dup2 = dupPhenoDat.duplicated(['barcodes','treatment'],keep=False)
dupKeep = dupPhenoDat.loc[dup2,:].sort_values('barcodes')
remove = [x for x in list(dupPhenoDat.index) if x not in list(dupKeep.index)]
pD = phenoDat.drop(remove,axis = 0)
print(pD.shape)
len(set(pD['barcodes']))

remove = []
for i,obsName in enumerate(adata.obs.index):
    if obsName not in list(pD['barcodes']):
        remove.append(i)
remove

#remove these indices
X = adata.X.toarray()
X = np.delete(X,remove,axis=0)
X = scp.sparse.csr_matrix(X)
S = adata.layers['spliced'].toarray()
S = np.delete(S,remove,axis=0)
S = scp.sparse.csr_matrix(S)
U = adata.layers['unspliced'].toarray()
U = np.delete(U,remove,axis=0)
U = scp.sparse.csr_matrix(U)

obsNames = [x for i,x in enumerate(adata.obs.index) if i not in remove]
var_names = list(adata.var.index)
concatAdata = ad.AnnData(X,
              {'obs_names': obsNames},
              {'var_names': var_names},
               layers={'spliced':S,
                      'unspliced':U})
filename = 'finalCleanConcatanatedData.h5ad'
print('final adata is:')
print(concatAdata)
print('writing file ' + filename)
concatAdata.write(filename)

adata = ad.read_h5ad('finalCleanConcatanatedData.h5ad')
adata

https://scvelo-notebooks.readthedocs.io/DentateGyrus.html

t2g = pd.read_table('homo_sapiens/transcripts_to_genes.txt',header=None)
t2g.head()

geneNameDict = dict(zip(t2g[1].values,t2g[2].values))
geneNames = []
for geneId in adata.var.index:
    geneNames.append(geneNameDict[geneId])
geneNames[:5]

adata.var['gene'] = geneNames
adata.var.head()

cell = []
pD = pD.set_index('barcodes')

cellDict = dict(zip(pD.index.values,pD['cell'].values))
treatmentDict = dict(zip(pD.index.values,pD['treatment'].values))
replicateDict = dict(zip(pD.index.values,pD['replicate'].values))
cells = []
treatments = []
replicates = []
for barcode in adata.obs.index:
    cells.append(cellDict[barcode])
    treatments.append(treatmentDict[barcode])
    replicates.append(treatmentDict[barcode])

adata.obs['cell'] = cells
adata.obs['treatment'] = treatments
adata.obs['replicate'] = replicates
adata

scv.pp.filter_and_normalize(adata, min_shared_counts=30, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

adata

scv.tl.velocity(adata)

scv.tl.velocity_graph(adata)
adata

sc.tl.umap(adata)

adata.write_h5ad('finalCleanadataWithVelocity.h5ad')

scv.pl.velocity_embedding_stream(adata, basis='umap',color='cell',save='velocityStream.png')

scv.pl.velocity_embedding_stream(adata, basis='umap',color='treatment',save='velocityStreamTreatment.png',legend_loc='right margin')

scv.pl.velocity_embedding(adata, basis='umap', arrow_length=2, arrow_size=1.5, dpi=150,color='cell')

scv.tl.recover_dynamics(adata)

adata.write_h5ad('finalCleanWithDynamics.h5ad')

df = adata.var
df.head()
df['geneId'] = df.index
df = df.set_index('gene')
df.head()
adata.var = df

adata

scv.tl.recover_latent_time(adata)

scv.pl.scatter(adata, color='latent_time', fontsize=24, size=100,
               color_map='gnuplot', perc=[2, 98], colorbar=True, rescale_color=[0,1],save='latentTime.png')

top_genes = adata.var_names[adata.var.fit_likelihood.argsort()[::-1]][:300]
scv.pl.heatmap(adata, var_names=top_genes, tkey='latent_time', n_convolve=100, col_color='cell',save='topGenesColorCell.png')

scv.pl.scatter(adata, basis=top_genes[5:10], legend_loc='none',
               size=80, frameon=True, ncols=5, fontsize=20,color='treatment',save='top5-10genesColByTreatment.png')

scv.pl.scatter(adata, basis=top_genes[5:10], legend_loc='right margin',
               size=80, frameon=True, ncols=5, fontsize=20,color='cell',save='top5-10GenesColByCell.png')

scv.pl.velocity_embedding_stream(adata, basis='umap',color='latent_time',save='velocityStreamLatentTime.png',legend_loc='right margin')

scv.pl.velocity_embedding_stream(adata, basis='umap',color='velocity_pseudotime',save='velocityStreamVelocityPseudoTime.png',legend_loc='right margin')

adata.var.fit_likelihood.argsort()[::-1][:300]