目录
环境导入
关键函数定义
运行前设定
数据准备
正式运行与保存
可视化与概率调整
偶然发现的一个好用的transfer方法,计算效率相当高,解了我的燃眉之急hh
原方法来自由以色列耶路撒冷希伯来大学的Mor Nitzan、美国麻省理工学院-哈佛大学博德研究所Aviv Regev和Johanna Klughammer共同通讯发表在Nature Biotechnology的研究成果。
环境导入
所需要的安装条件已经详细列出https://github.com/simonwm/tacco
这里不再赘述,我们先导入需要的环境包
import os
import sys
import cv2
import seaborn as sns
import anndata as ad
import argparse
import pandas as pd
import numpy as np
import tacco as tc
import scanpy as sc
import matplotlib.pyplot as plt
from matplotlib.transforms import Affine2D
关键函数定义
一些后续需要使用到的可视化以及修改函数定义(直接使用可跳过)
def cluster_small_multiples(adata, clust_key, size=10, frameon=False, **kwargs):fig, axes = plt.subplots(figsize=(7*len(celltype_order), 3))tmp = adata.copy()for i, clust in enumerate(adata.obs[clust_key].cat.categories):tmp.obs[clust] = adata.obs[clust_key].isin([clust]).astype('category')tmp.uns[clust+'_colors'] = ['#d3d3d3', adata.uns[clust_key+'_colors'][i]]sc.pl.scatter(tmp, groups=tmp.obs[clust].cat.categories[1:].values, color=adata.obs[clust_key].cat.categories.tolist(), size=size, frameon=frameon, **kwargs)def markers_plot(adata, clust_key, df, **kwargs):adata = adata.copy()# Normalizing to median total countssc.pp.normalize_total(adata)# Logarithmize the datasc.pp.log1p(adata)markers = {i['celltype']: [j for j in i['marker'].split(',') if j in adata.var_names] for (_, i) in df.iterrows()}sc.pl.dotplot(adata, markers, groupby=clust_key, **kwargs)def update_celltype_proportions(df, modifications, ct_lowres='celltype_lowres', ct_hires="celltype_hires"):# 计算每个粗注释的细胞比例celltype_proportions = df[ct_lowres].value_counts(normalize=True)# 对每个粗注释分组,然后计算每个粗注释中其细分群的细胞比例hires_proportions = df.groupby(ct_lowres)[ct_hires].value_counts(normalize=True)grouped_proportions = hires_proportions.mul(celltype_proportions, level=0)if modifications is None:return celltype_proportions, grouped_proportions# 计算未被修改的细胞类型的当前总比例current_remaining_proportion = celltype_proportions.drop(modifications.keys()).sum()# 计算修改后剩余的比例remaining_proportion = 1 - sum(modifications.values())# 更新未被修改的细胞类型的比例celltype_proportions.update(celltype_proportions.drop(modifications.keys()) * remaining_proportion / current_remaining_proportion)# 更新指定的粗注释细胞比例celltype_proportions.update(pd.Series(modifications))# 调整每个细分群的细胞比例grouped_proportions = hires_proportions.mul(celltype_proportions, level=0)return celltype_proportions, grouped_proportions
运行前设定
这里假设我要批量用一个注释好的单细胞转录组数据,去映射一组空间数据
sc_path = "/data/your_annotated_sc_data.h5ad"
sp_dir = "/data/your_sp_data_path/"
output_dir = "/data/output/"
os.makedirs(output_dir, exist_ok=True) # 创建输出目录,如果不存在
数据准备
读取单细胞参考和所有空间数据,这里prob默认按照原单细胞注释中细胞比例来设定,但后续可以手动调整这个比例,从而调整映射效果
# 读取单细胞参考数据
reference = sc.read_h5ad(sc_path)
reference = reference[reference.obs['annotation'] != 'Other']
table = reference.obs["annotation"].value_counts()
prob = reference.obs["annotation"].value_counts() / reference.obs["annotation"].size# 获取空间转录组数据目录下的所有文件
sp_files = [f for f in os.listdir(sp_dir) if f.endswith('.h5ad')]
正式运行与保存
for sp_file in sp_files:print(f"Processing {sp_file}...")sp_path = os.path.join(sp_dir, sp_file)adata = sc.read_h5ad(sp_path)adata.X = adata.X.astype('float32')# 使用 Tacco 进行细胞类型注释tc.tl.annotate(adata, reference, "annotation",result_key='pred_celltype', annotation_prior=prob, verbose=False)# 提取预测结果prediction = adata.obsm['pred_celltype'].copy()prediction['pred_celltype'] = prediction.idxmax(1)prediction['max_score'] = prediction.iloc[:, :-1].max(1)pred_prob = prediction['pred_celltype'].value_counts(normalize=True)mse = ((pred_prob - prob) ** 2).mean()remove_celltype = pred_prob[pred_prob == 0].index.tolist()# 更新 adata 的注释信息adata.obs['pred_celltype'] = prediction['pred_celltype']adata.obs["pred_celltype"] = adata.obs["pred_celltype"].cat.remove_categories(remove_celltype)# 保存预测后的 h5ad 文件output_h5ad_path = os.path.join(output_dir, f"{sp_file.split('.')[0]}_predicted.h5ad")adata.write(output_h5ad_path)print("All files processed successfully.")
这里就得到了所有预测标签,放置在‘pred_celltype’标签下
可视化与概率调整
这里可以轻松地通过注释标签可视化预测的标签空间分布情况
adata = sc.read_h5ad("/data/your_predicted_sp.h5ad")
sc.pl.spatial(adata,color='pred_celltype',spot_size=100,title=f'Tacco Predicted Cell Types',legend_loc='right margin',palette='tab20',show=True,
# save=output_png_path)
同时你也可以输出上述默认的“prob”检查预测的细胞比例,如果有需要调整的,或者有异常值的,可以修改之后重新运行如下
# 修改概率
prob["label1"] = 0.5
prob["label2"] = 0.1# 重新运行
tc.tl.annotate(adata, reference, "annotation",result_key='pred_celltype', annotation_prior=prob, verbose=False)# 提取预测结果
prediction = adata.obsm['pred_celltype'].copy()
prediction['pred_celltype'] = prediction.idxmax(1)
prediction['max_score'] = prediction.iloc[:, :-1].max(1)
pred_prob = prediction['pred_celltype'].value_counts(normalize=True)
mse = ((pred_prob - prob) ** 2).mean()
remove_celltype = pred_prob[pred_prob == 0].index.tolist()# 更新 adata 的注释信息
adata.obs['pred_celltype'] = prediction['pred_celltype']
adata.obs["pred_celltype"] = adata.obs["pred_celltype"].cat.remove_categories(remove_celltype)
所有的迁移方法都是服务于我们的科学问题,本身并无好坏之分!主要这个方法真的计算很快,有相关需求的朋友大可以尝试一下~