細(xì)胞不會(huì)跳舞我不知道，但似乎nature biotech的文章代碼會(huì)有毒异吻。。。

CellDancer是最新發(fā)表在nature biotechnology上的一個(gè)類似scVelo的計(jì)算RNA Velocity的程序题涨，原文詳見：
https://www.nature.com/articles/s41587-023-01728-5

文本中說本文的創(chuàng)新之處在于可以克服scVelo不能分析的一些數(shù)據(jù)，如含有多個(gè)cell lineage時(shí)总滩，或者有基因表達(dá)不符合前提假設(shè)纲堵，有boost時(shí)，Velocity分析會(huì)出錯(cuò)闰渔。Celldancer 號(hào)稱可以解決這個(gè)問題席函。恰好有自己的數(shù)據(jù)一組，scVelo 分析會(huì)出現(xiàn)明顯的反向冈涧，決定嘗試一下茂附，看看虛實(shí)。

首先就發(fā)現(xiàn)這個(gè)包不能直接用于我的數(shù)據(jù)炕舵，celldancer需要數(shù)據(jù)以csv的形式讀入何之，雖然也提供了一個(gè)從anndata里提取splicing數(shù)據(jù)的功能：adata_to_df_with_embed，實(shí)測(cè)這個(gè)功能 暫時(shí)（截止至：04-04-2023） 有點(diǎn)雞肋咽筋，暫時(shí)不能滿足實(shí)際需要溶推。

一開始以為是小毛病，看完后發(fā)現(xiàn)必須得改奸攻，因?yàn)榱钊艘庀氩坏降氖撬馕＃瓉淼牟糠执a是這樣寫的：

    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)

竟然要求每循環(huán)一個(gè)基因，就寫盤一次睹耐，如果你幾萬個(gè)基因辐赞，就要寫盤幾萬次？硝训？响委？這不是TMD神經(jīng)病嗎。多跑幾次程序你的電腦就廢了窖梁。

到這里隱約感覺這可能又是一個(gè)坑赘风，代碼寫得這么爛，不太相信效果會(huì)好纵刘。邀窃。。假哎。這年頭發(fā)文章不吹牛都會(huì)死瞬捕。鞍历。。

另外非常不喜歡那個(gè)一行一行的打印處理進(jìn)度肪虎，打印了幾萬行劣砍，除了刷屏有毛用？用tqdm改了個(gè)進(jìn)度條扇救。

改動(dòng)后如下：

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)

下面開始搞
首先秆剪，導(dǎo)入所需要的庫，也不知道真的需不需要這么多

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

讀入數(shù)據(jù)

adata = sc.read_h5ad('/home/Documents/integrated_20L_with_splicing.h5ad')

用魔改后的adata_to_df_with_embed將我的保存于anndata的剪切數(shù)據(jù)提取出來爵政，此處不輸出csv 文件仅讽， 直接從內(nèi)存調(diào)用，你也可選擇將其保存成csv钾挟，而不是原來的強(qiáng)制保存洁灵。

cell_type_u_s =cdutil.adata_to_df_with_embed(adata,\
                              us_para=['unspliced','spliced'],\
                              cell_type_para='seurat_clusters',\
                              embed_para='X_umap') # 

該過程比較慢，取決于基因數(shù)目掺出，你也可以加入一個(gè)感興趣基因的列表徽千，單獨(dú)處理這些基因而非全部基因.

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

這一步貌似也有點(diǎn)問題，容易出錯(cuò)汤锨，似乎減少基因數(shù)目和線程數(shù) 有幫助双抽。問題似乎有點(diǎn)復(fù)雜，暫且跳過闲礼。

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()

結(jié)果似乎不咋地牍汹，估計(jì)如果挑幾個(gè)重要的基因，可能會(huì)有所改善柬泽。有待進(jìn)一步研究慎菲。

Figure_4.png