官網(wǎng)的代碼演示差不多已經(jīng)能夠完成我的部分需求，但是有部分函數(shù)他介紹的不太完善疆导，我進(jìn)行了一些修改，并對(duì)腳本進(jìn)行了封裝葛躏。

step1 數(shù)據(jù)讀取

1.讀取loom文件澈段，將所有的loom文件都放在一個(gè)名為loom_file的文件夾里，并進(jìn)行讀取舰攒。

loom_data = {}
index_sample = []

# 遍歷 loom 文件列表
for file in loom_file:
    # 檢查文件是否為 loom 文件
    if '.loom' in file:
        # 打印文件路徑
        print(f"{loom_dir}/{file}")
        
        # 讀取 loom 文件數(shù)據(jù)
        loom_temp_data = sc.read(f"{loom_dir}/{file}", cache=True)
        
        # 確卑芨唬基因名稱(chēng)唯一
        loom_temp_data.var_names_make_unique()
        
        # 替換觀(guān)測(cè)名稱(chēng)中的特殊字符,這段不固定，需要根據(jù)你的loom里面的barcode和單細(xì)胞數(shù)據(jù)里面的barcode進(jìn)行統(tǒng)一
        loom_temp_data.obs_names = loom_temp_data.obs_names.str.replace("x", "")
        
        # 將樣本索引添加到列表中
        index_sample += list(loom_temp_data.obs_names)
        
        # 將數(shù)據(jù)存儲(chǔ)在 loom 數(shù)據(jù)字典中
        loom_data[file] = loom_temp_data

2.讀取rds或者h(yuǎn)5ad文件摩窃，rds文件需要提前轉(zhuǎn)換成h5ad文件

sc_data = scv.read('test.h5ad')
loom_data_merge = loom_data[list(loom_data.keys())[0]]

# 遍歷 loom_data 中的其他 loom 數(shù)據(jù)
for i in range(1, len(loom_data)):
    # 獲取下一個(gè) loom 數(shù)據(jù)
    next_data = loom_data[list(loom_data.keys())[i]]
    
    # 將當(dāng)前 loom 數(shù)據(jù)與下一個(gè) loom 數(shù)據(jù)連接
    loom_data_merge = loom_data_merge.concatenate(next_data)

loom_data_merge.obs_names=index_sample

adata = scv.utils.merge(sc_data, loom_data_merge)

3.讀取一下每個(gè)cluster的剪切和未剪切的比例,也可以觀(guān)察其他的信息兽叮，比如sample芬骄，group，celltype等

groupby='clusters'
output='scvelocity_outdir'
scv.pl.proportions(adata, groupby=groupby,
save=f'{output}/spliced_proportions_groupby_{groupby}.pdf')

step2 數(shù)據(jù)處理

1.數(shù)據(jù)預(yù)處理鹦聪，有些人可能之前已經(jīng)用scanpy做過(guò)一些數(shù)據(jù)的預(yù)處理账阻，會(huì)疑問(wèn)這一步是否還有必要進(jìn)行，作者給出的答案是就算之前已經(jīng)做過(guò)了預(yù)處理泽本，但是也沒(méi)有對(duì)spliced和unspliced進(jìn)行處理淘太，還是要運(yùn)行這一步，該函數(shù)會(huì)自動(dòng)判斷數(shù)據(jù)是否進(jìn)行過(guò)標(biāo)準(zhǔn)化规丽，如今已經(jīng)進(jìn)行了蒲牧，則不會(huì)再log處理。

scv.pp.filter_and_normalize(adata赌莺，min_shared_counts=20, n_top_genes=2000)

scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

step3 速度計(jì)算

scvelo支持動(dòng)態(tài)模型和穩(wěn)態(tài)模型冰抢，使用動(dòng)態(tài)模型比穩(wěn)態(tài)模型能夠獲得更多的信息。
通過(guò)代表基因推算細(xì)胞的未來(lái)分化方向艘狭，并投影到低緯度上挎扰。

scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata,vkey='dynamical_velocity', mode='dynamical',n_jobs=5)
scv.tl.velocity_graph(adata,vkey='dynamical_velocity',n_jobs=5)

step4 結(jié)果可視化

作者提供三種速度投影到umap的展示方式
細(xì)胞水平 scv.pl.velocity_embedding，
網(wǎng)格線(xiàn) scv.pl.velocity_embedding_grid
流線(xiàn)型 scv.pl.velocity_embedding_stream

#第一張
scv.pl.velocity_embedding_stream(adata, basis='umap',key='dynamical_velocity',
title=f'velocity_embedding_stream_groupby_{groupby}',color=groupby,
legend_fontsize=18, legend_loc='right margin',colorbar=TRUE,dpi=300,
save=f'{output}/velocity_embedding_stream_groupby{groupby}.svg')

image.png

#第二張
scv.pl.velocity_embedding(adata, basis='umap',vkey='dynamical_velocity',
title=f'velocity_embedding_groupby_{groupby}',
arrow_length=2,color=groupby,arrow_size=1, legend_loc='right margin',
dpi=300,save=f'{output}/velocity_embedding_groupby_{groupby}.svg')

image.png

#第三張
scv.pl.velocity_embedding_grid(adata, basis='umap',vkey='dynamical_velocity',
color=groupby,dpi=300,title=f'velocity_embedding_grid_groupby_{groupby}',
legend_loc='right margin',save=f'{output}/velocity_embedding_grid_groupby_{groupby}.pdf')

image.png

step5 latent time

    #計(jì)算潛伏時(shí)間
    scv.tl.recover_latent_time(adata, vkey='dynamical_velocity')
    scv.pl.scatter(adata,basis='umap',
                   color='latent_time',
                   color_map='gnuplot',
                   colorbar=True,
                   perc=[2, 98], #指定連續(xù)著色的百分比缓升。
                   rescale_color=[0,1], #顏色重新縮放的邊界鼓鲁，設(shè)置顏色條的最小/最大值。
                   dpi = 300,
                   save=f'{output}/latent_time_based_dynamical.pdf')

image.png

    #取出Top300的likelihood genes信息港谊，并保存骇吭，繪制熱圖
    TopGenes = adata.var['fit_likelihood'].sort_values(ascending=False).index[:300]
    df = pd.DataFrame(TopGenes)
    df=df.transpose()
    df.to_csv(f'{output}/Top300_likelihood_genes_data.csv')
   #繪制熱圖
    scv.pl.heatmap(adata, var_names=TopGenes, 
                   sortby='latent_time', 
                   col_color=groupby,#選擇作為注釋條的列名
                   n_convolve=100,#平滑曲線(xiàn)的核數(shù)
                   dpi=400, #這個(gè)圖片很大，不要保存太大的dpi
                   save = f'{output}/latent_time_heatmap_groupby_{groupby}.pdf')

image.png

step6 Top-likelihood genes

動(dòng)態(tài)模型中的高似然基因可以認(rèn)為是驅(qū)動(dòng)基因歧寺。這些基因展現(xiàn)出明顯的動(dòng)態(tài)特征燥狰，并且可以明顯的地被檢測(cè)到

scv.pl.scatter(adata2,
        basis=top_genes[:15],
        ncols=5,
        use_raw=False,
        frameon=False,
        legend_loc='right margin',
        save=f'{output}/Top-likelihood_genes_groupby_{groupby}.pdf')