ixxmu
commented
4 months ago 高级替身！拟时序神器sctour，出图神似RNA velocity，兼顾时序和表达量给出降维 by 生信作曲家
前言
拟时序分析的方法已经很多了，我们之前向大家介绍了VECTOR算法，实际上可以类似于RNA速率分析，做出向量，且不需要bam文件。
今天我们介绍一种又能够做出更类似于RNA速率分析的拟时序算法，sctour !
sctour算法是一种用于稳健推断和精确预测细胞动态的深度学习架构，其优势在于，可以兼顾拟时间和表达了进行重降维聚类（一般拟时序不会给新的降维方法），并且，它出图真的非常类似RNA速率！
不管心不心动，小编已经非常激动啦！
代码实操
# pip install sctour
import sctour as sct
import scanpy as sc
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
adata = sc.read('G:/scouter/scRNA2.h5ad')
sc.pl.umap(adata, color=['celltype', 'Sample batch'], legend_loc='on data')
## 恶性细胞已经分好群了
# 计算
import os
os.environ['NUMBA_DISABLE_INTEL_SVML'] = '1'
import skmisc
#计算每个细胞中检测到的基因数量。这一步在scTour模型训练之前是必需的。
sc.pp.calculate_qc_metrics(adata, percent_top=None, log1p=False, inplace=True)
# Select highly variable genes. 你可能需要 pip install --user scikit-misc==0.2.0
sc.pp.highly_variable_genes(adata, flavor='seurat_v3', n_top_genes=1000, subset=True)
#训练scTour模型。默认的损失模式是负二项条件似然（nb），
# 默认情况下，当细胞总数少于10,000时，用于训练模型的细胞比例设置为0.9，而当细胞总数大于10,000时，该比例设置为0.2
# 用户可以通过使用参数percent来调整比例（例如，percent=0.6）。
#较小的alpha_recon_lec往往会得出基于伪时间排列细胞的潜在表示，因此不能很好地分离具有相似发育顺序的细胞类型。
# 这两个参数的默认值都是0.5。用户可以根据数据集调整它们。
tnode = sct.train.Trainer(adata, loss_mode='nb', alpha_recon_lec=0.5, alpha_recon_lode=0.5)
tnode.train()
#伪时间
#基于训练好的模型推断发育伪时间。
adata.obs['ptime'] = tnode.get_time()
#推断细胞位置，相当于替换掉umap
#两个参数alpha_z和alpha_predz用于调整从变分推断和ODE求解器获得的潜在空间的权重。
# 较大的alpha_z会使潜在空间偏向于内在的转录组结构，而较大的alpha_predz则更能代表外在的伪时间排序。
# 用户可以根据自己的目的调整这两个参数。
#zs represents the latent z from variational inference, and pred_zs represents the latent z from ODE solver
#mix_zs represents the weighted combination of the two, which is used for downstream analysis
mix_zs, zs, pred_zs = tnode.get_latentsp(alpha_z=0.5, alpha_predz=0.5)
adata.obsm['X_TNODE'] = mix_zs
## 推断向量
adata.obsm['X_VF'] = tnode.get_vector_field(adata.obs['ptime'].values, adata.obsm['X_TNODE'])
# 可以保存
#adata.write('G:/scouter/sctour.h5ad')
#adata = sc.read('G:/scouter/sctour.h5ad')
#可视化，我们看到，sctour给出了新的聚类！
adata = adata[np.argsort(adata.obs['ptime'].values), :]
sc.pp.neighbors(adata, use_rep='X_TNODE', n_neighbors=15)
sc.tl.umap(adata, min_dist=0.1)
fig, axs = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))
sc.pl.umap(adata, color='celltype', ax=axs[0, 0], legend_loc='on data', show=False, frameon=False)
sc.pl.umap(adata, color='Sample batch', ax=axs[0, 1], show=False, frameon=False)
sc.pl.umap(adata, color='ptime', ax=axs[1, 0], show=False, frameon=False)
sct.vf.plot_vector_field(adata, zs_key='X_TNODE', vf_key='X_VF', use_rep_neigh='X_TNODE', color='celltype', show=False, ax=axs[1, 1], legend_loc='none', frameon=False, size=100, alpha=0.2)
plt.show()
版本及报错指南
往往同学们运行时会遭遇各种报错，我们下面列出来一些常见的报错及解决策略。
0. 示例数据
购买后管理员处获取
https://mp.weixin.qq.com/s/CBdljX2zfSoz-VZIWbcoRw