細胞軌跡分析方法有哪些？單細胞數(shù)據(jù)分析新思路(內(nèi)附代碼)

pip?install?-U?scvelo?#需要3.6以上版本python
pip?install?velocyto

– adata.X：數(shù)據(jù)矩陣
– adata.obs：細胞的注釋
–?adata.var：基因的注釋
– adata.nus：非結(jié)構(gòu)化注釋，例如圖
– adata.layers：附加層數(shù)據(jù)，如spliced/unspliced矩陣

import?seaborn?as?sns
import?scvelo?as?scv
import?pandas?as?pd
import?scanpy?as?sc
import?numpy?as?np
import?loompy

scv.settings.verbosity?=?3??#?show?errors(0),?warnings(1),?info(2),?hints(3)
scv.settings.presenter_view?=?True??#?set?max?width?size?for?presenter?view
scv.set_figure_params(‘scvelo’)??#?for?beautified?visualization

##?運行run10x
/software/python3.9.6/bin/velocyto?run10x??????????\
#?相應(yīng)的repeatmasker的重復(fù)序列注釋文件（可從USCU中獲?。?br> #?http://genome.ucsc.edu/cgi-bin/hgTables
-m?RMSK
cellranger_results??#?cellranger的輸出結(jié)果路徑
sample_name???????????#?樣本名
seurat_GTF??????????????#?用于cellranger定量的基因組注釋文件
##?輸出結(jié)果，在cellranger的結(jié)果中會創(chuàng)建單個樣本的velocyto文件夾，其中的loom文件即為最終結(jié)果
data/SC/Normal1-sc/velocyto:
Normal1-sc.loom
data/SC/Normal2-sc/velocyto:
Normal2-sc.loom
data/SC/Normal3-sc/velocyto:
Normal3-sc.loom
data/SC/Normal4-sc/velocyto:
Normal4-sc.loom

#?讀取單個樣本loom文件
Normal_loom?=?sc.read(‘./Normal1.loom’,?sparse=True,?cache=True)
#?合并同組樣本的loom文件
import?loompy
Normal_loom_list?=?[Normal1-sc.loom,
Normal2-sc.loom,
Normal3-sc.loom,
Normal4-sc.loom]
loompy.combine(Normal_loom_list,?output_file?=?‘Normal_combined.loom’)
#?讀取
Normal_loom?=?sc.read(‘./Normal_combined.loom’,?sparse=True,?cache=True)

library(Seurat)
library(tidyverse)
seurat_obj@meta.data?%>%?rownames_to_column(var?=?‘barcodes’)?%>%
write.table(‘seurat_meta_df.txt’,?sep?=?‘\t’,?quote?=?F)
seurat_obj@reductions$umap@cell.embeddings?%>%?rownames_to_column(var?=?‘barcodes’)?%>%
write.table(‘seurat_umap.txt’,?sep?=?‘\t’,?quote?=?F)
def?using_seurat_meta_umap(loom_file,?seurat_meta_df,?seurat_umap):
#?讀入seurat的meta.data以及umap降維結(jié)果
seurat_meta_data_raw?=?pd.read_csv(seurat_meta_df)
seurat_umap_raw?=?pd.read_csv(seurat_umap)
#?loom文件原始的index順序
anndata_index?=?loom_file.obs.index.values
#?合并
meta_data?=?loom_file.obs.reset_index().?\
merge(seurat_meta_data_raw,
left_on=”CellID”,
right_on=”barcodes”,?how=’left’).?\
set_index(‘CellID’).reindex(anndata_index)
meta_data[‘seurat_clusters’]?=?meta_data[‘seurat_clusters’].?\
apply(lambda?x:?“%d”?%?int(x)?if?~np.isnan(x)?else?x).astype(‘category’)
loom_file.obs?=?meta_data
loom_file?=?loom_file[loom_file.obs.dropna().index,?:]
loom_file.obsm[‘umap’]?=?seurat_umap_raw.set_index(‘barcodes’).?\
reindex(loom_file.obs.index).to_numpy().astype(‘float64’)
return?loom_file
Normal_loom?=?using_seurat_meta_umap(Normal_loom,
‘seurat_meta_df.txt’,
‘seurat_umap.txt’)
#?查看是否替換成功，與seurat結(jié)果進行比較
sc.pl.umap(Normal_loom,?color=’seurat_clusters’)

scv.pl.proportions(Normal_loom,?groupby=’seurat_clusters’)

#?預(yù)處理
scv.pp.filter_genes(Normal_loom,?min_shared_counts=20)
scv.pp.normalize_per_cell(Normal_loom)
scv.pp.filter_genes_dispersion(Normal_loom,?n_top_genes=2000)
scv.pp.log1p(Normal_loom)
scv.pp.filter_and_normalize(Normal_loom,
min_shared_counts=20,
n_top_genes=2000)
#?運行scVelo動態(tài)模型，如果算力不夠可以考慮靜態(tài)或者隨機方法
scv.tl.recover_dynamics(Normal_loom,?n_jobs=30)
scv.tl.velocity(Normal_loom,?mode=’dynamical’)
scv.tl.velocity_graph(Normal_loom,?n_jobs=n_jobs)
#?保存計算結(jié)果
Normal_loom.write(‘Normal_dy_model.h5ad’,?compression=’gzip’)
scv.pl.velocity_embedding_stream(Normal_loom,
basis=’umap’,
color=’seurat_clusters’)
scv.pl.velocity_embedding(Normal_loom,
arrow_length=3,
arrow_size=2,
show=False,
color=’seurat_clusters’)

scv.tl.latent_time(Normal_loom)
scv.pl.scatter(Normal_loom,?color=’latent_time’,
color_map=’gnuplot’,?size=80,?show=False)

top_genes?=?Normal_loom.var[‘fit_likelihood’].?\
sort_values(ascending=False).index[:300]

scv.pl.heatmap(Normal_loom,
var_names=top_genes,?sortby=’latent_time’,
col_color=’seurat_clusters’,?n_convolve=100)

scv.pl.scatter(Normal_loom,?basis=top_genes[:15],
ncols=5,?frameon=False,
color=’seurat_clusters’)

scv.tl.paga(Normal_loom,?groups=’seurat_clusters’)

scv.pl.paga(Normal_loom,?basis=’umap’,
size=50,?alpha=.1,
min_edge_width=2,?node_size_scale=1.5)