细胞不会跳舞我不知道，但似乎nature biotech的文章代码会有毒。。。
CellDancer是最新发表在nature biotechnology上的一个类似scVelo的计算RNA Velocity的程序，原文详见：
https://www.nature.com/articles/s41587-023-01728-5
文本中说本文的创新之处在于可以克服scVelo不能分析的一些数据，如含有多个cell lineage时，或者有基因表达不符合前提假设，有boost时，Velocity分析会出错。Celldancer 号称可以解决这个问题。恰好有自己的数据一组，scVelo 分析会出现明显的反向，决定尝试一下，看看虚实。
首先就发现这个包不能直接用于我的数据，celldancer需要数据以csv的形式读入，虽然也提供了一个从anndata里提取splicing数据的功能：adata_to_df_with_embed，实测这个功能 暂时（截止至：04-04-2023） 有点鸡肋，暂时不能满足实际需要。
一开始以为是小毛病，看完后发现必须得改，因为令人意想不到的是，原来的部分代码是这样写的：
for i,gene in enumerate(gene_list):
data_onegene = adata_to_raw_one_gene(adata, us_para=us_para, gene=gene)
if i==0:
data_onegene.to_csv(save_path,header=True,index=False)
else:
data_onegene.to_csv(save_path,mode='a',header=False,index=False)
竟然要求每循环一个基因，就写盘一次，如果你几万个基因，就要写盘几万次？？？这不是TMD神经病吗。多跑几次程序你的电脑就废了。
到这里隐约感觉这可能又是一个坑，代码写得这么烂，不太相信效果会好。。。。这年头发文章不吹牛都会死。。。
另外非常不喜欢那个一行一行的打印处理进度，打印了几万行，除了刷屏有毛用？用tqdm改了个进度条。
改动后如下：
def adata_to_df_with_embed(adata,
us_para=['Mu', 'Ms'],
cell_type_para='celltype',
embed_para='X_umap',
save_path='cell_type_u_s_sample_df.csv',
gene_list=None):
def adata_to_raw_one_gene(data, us_para, gene):
data2 = data[:,data.var.index.isin([gene])].copy()
n=len(data2)
u0 = data2.layers[us_para[0]][:,0].copy().astype(np.float32)
u0 = scipy.sparse.csr_matrix.todense(u0).tolist()
s0 = data2.layers[us_para[1]][:,0].copy().astype(np.float32)
s0 = scipy.sparse.csr_matrix.todense(s0).tolist()
raw_data = pd.DataFrame({'gene_name':[gene]*n, 'unsplice':u0, 'splice':s0})
raw_data.splice = [i[0] for i in raw_data.splice]
raw_data.unsplice = [i[0] for i in raw_data.unsplice]
return(raw_data)
if gene_list is None: gene_list=adata.var.index
dfs = []
for gene in tqdm(gene_list):
global combined
data_onegene = adata_to_raw_one_gene(adata, us_para=us_para, gene=gene)
data_onegene.sort_index(inplace=True)
dfs.append(data_onegene)
combined = pd.concat(dfs).reset_index(drop=True)
# cell info
gene_num=len(gene_list)
cellID=pd.DataFrame({'cellID':adata.obs.index})
celltype_meta=adata.obs[cell_type_para].reset_index(drop=True)
celltype=pd.DataFrame({'clusters':celltype_meta})#
embed_map=pd.DataFrame({'embedding1':adata.obsm[embed_para][:,0],'embedding2':adata.obsm[embed_para][:,1]})
# embed_info_df = pd.concat([embed_info]*gene_num)
embed_info=pd.concat([cellID,celltype,embed_map],axis=1)
embed_raw=pd.concat([embed_info]*gene_num)
embed_raw=embed_raw.reset_index(drop=True)
raw_data=pd.concat([combined,embed_raw],axis=1)
return(raw_data)
下面开始搞
首先，导入所需要的库，也不知道真的需不需要这么多
import pandas as pd
import celldancer.utilities as cdutil
import scanpy as sc
import os
import sys
import glob
import math
import matplotlib.pyplot as plt
import celldancer as cd
import celldancer.cdplt as cdplt
from celldancer.cdplt import colormap
读入数据
adata = sc.read_h5ad('/home/Documents/integrated_20L_with_splicing.h5ad')
用魔改后的adata_to_df_with_embed将我的保存于anndata的剪切数据提取出来，此处不输出csv 文件， 直接从内存调用，你也可选择将其保存成csv，而不是原来的强制保存。
cell_type_u_s =cdutil.adata_to_df_with_embed(adata,\
us_para=['unspliced','spliced'],\
cell_type_para='seurat_clusters',\
embed_para='X_umap') #
该过程比较慢，取决于基因数目，你也可以加入一个感兴趣基因的列表，单独处理这些基因而非全部基因.
gene_list=['Rora','Fgfr3','Wee1','Cux1','Slc6a6','Hlf','Myof','Dmkn','Sema3d']
loss_df, cellDancer_df=cd.velocity(cell_type_u_s,\
gene_list=gene_list,\
permutation_ratio=0.125,\
n_jobs=2)
cellDancer_df
这一步貌似也有点问题，容易出错，似乎减少基因数目和线程数 有帮助。问题似乎有点复杂，暂且跳过。
import seaborn as sns
colevels=adata.obs.seurat_clusters.unique()
cellcolor=dict(zip(colevels,sns.color_palette("husl", len(colevels)).as_hex()))
# compute cell velocity
cellDancer_df=cd.compute_cell_velocity(cellDancer_df=cellDancer_df, projection_neighbor_choice='gene', expression_scale='power10', projection_neighbor_size=10, speed_up=(100,100))
# plot cell velocity
fig, ax = plt.subplots(figsize=(10,10))
cdplt.scatter_cell(ax,
cellDancer_df,
colors=cellcolor,
alpha=0.5,
s=10,
velocity=True,
legend='on',
min_mass=15,
arrow_grid=(20,20),
custom_xlim=[-6,13],
custom_ylim=[2,16], )
ax.axis('off')
plt.show()
import random
# set parameters
dt = 0.05
t_total = {dt:int(10/dt)}
n_repeats = 10
# estimate pseudotime
cellDancer_df = cd.pseudo_time(cellDancer_df=cellDancer_df,
grid=(30,30),
dt=dt,
t_total=t_total[dt],
n_repeats=n_repeats,
speed_up=(100,100),
n_paths = 3,
plot_long_trajs=True,
psrng_seeds_diffusion=[i for i in range(n_repeats)],
n_jobs=2)
# plot pseudotime
fig, ax = plt.subplots(figsize=(6,6))
im=cdplt.scatter_cell(ax,cellDancer_df, colors='pseudotime', alpha=0.5, velocity=False, custom_xlim=(-5,11), custom_ylim=(4,18))
ax.axis('off')
plt.show()
结果似乎不咋地，估计如果挑几个重要的基因，可能会有所改善。有待进一步研究。
Figure_4.png
