ixxmu
commented
1 year ago 整合单细胞数据和Bulk数据的多种方法（三）：Python包Bulk2Space by 生信随笔
接着前两个帖子：
除了scAB和Scissor这两个工具以外，衔接单细胞数据和Bulk测序数据的方法还有scPrognosis 和DEGAS 等工具。scPrognosis是2020年发表在PLoS Comput Biol上的一个R包，专门针对乳腺癌的bulk和单细胞数据，目前引用了11次。然而作者的github（https://github.com/XiaomeiLi1/scPrognosis）并没有给出一个示例文档，我也没有那么多精力去研究里面的源代码，所以我只能放弃这个R包。。。此外，DEGAS算法是在2022年发表于Genome Med上的一个python工具，基于深度学习算法，依赖GPU。然而，我在安装这个算法的python依赖包functools时，一直报错，尽管严格按照教程的软件版本，我换了三台服务器依然是报错，所以我也放弃了这个算法。。。本以为本系列推文到此就结束了，我突然想到前段时间发表在NC上的一个工具Bulk2Space，大致目的是把bulk数据反卷积映射到空间转录组数据上。
Bulk2Space发表在2022年10月份的Nature communications 上，题为《De novo analysis of bulk RNA-seq data at spatially resolved single-cell resolution》。Bulk2Space ( https://github.com/ZJUFanLab/bulk2space)是一种基于深度学习框架的空间反卷积算法（python环境），可以使用现有的单细胞和空间转录组学参考同时揭示Bulk RNA-seq 数据的空间和细胞异质性。
不得不吐槽一句，这个算法需要基于Bulk数据+单细胞数据+空转数据，然后去推算Bulk数据的单细胞矩阵和空间位置信息，从这篇文章出来之初，我就在想这个算法有啥用。。。到今天重新读了这篇文章及其代码，我还是没搞明白这套组合拳能够去探讨一个怎么样的生物学问题？
也可能是我的认识有限，因此非常欢迎有想法的朋友在文末进行留言讨论！
一. Bulk2Space算法的工作流程
简言之，Bulk2Space算法需要输入Bulk数据+单细胞数据（表达矩阵+细胞类型信息）+空转数据（表达矩阵+细胞空间位置分布信息），具体格式和内容见下文实战部分。而且，从示例数据看来，每次只能针对一个样本/患者的Bulk数据进行推断，推断逻辑具体是：
第一步反卷积（Deconvolution）：利用单细胞参考数据对Bulk数据进行反卷积从而构建一个模拟的单细胞矩阵。注意这里的反卷积不同于CIBERSORTx、MuSiC 和贝叶斯棱镜算法（BayesPrism），这些算法结合bulk数据和单细胞数据，然后直接给我们bulk数据里的细胞类型的比例情况。而Bulk2Space是直接输出一个基于Bulk数据推断的“单细胞表达矩阵”。 第二步空间映射（Spatial mapping）：在假定的“单细胞表达矩”的基础上，将每个细胞映射到空间位置上。本质上，这和传统的单细胞数据映射到空转数据是一个逻辑。 具体的算法逻辑见b,c和d。
二. Bulk2Space的代码实现
本文使用PDAC示例数据集进行展示，示例数据在谷歌云盘https://drive.google.com/file/d/1xB-Gk_KLxQA320-tycJp4CFHA66zF3LE/view。考虑到很多人无法使用谷歌，我也把示例数据上传到了百度云盘：https://pan.baidu.com/s/1Z2FoobwGmKgdk6Ijky0n9Q?pwd=1234 提取码：1234
Step0. 环境部署
Bulk2Space算法依赖python环境，需要在Linux和python IDE下进行分析：
conda create -n bulk2space python=3.8
conda activate bulk2space
pip install torch --extra-index-url https://download.pytorch.org/whl/cu116
github上的R包下载至服务器，然后：
cd bulk2space-main
pip install -r requirements.txt
python setup.py build
python setup.py install
快速启动
要使用Bulk2Space，我们需要5个格式化的“.csv”文件作为输入(即由pandas读取)。示例数据在tutorial/data/example_data文件夹中（https://github.com/ZJUFanLab/bulk2space/blob/main/tutorial/data/example_data）。
如果您选择基于spot的数据(如10x Genomics、ST或Slide-seq等)作为空间参考数据集，可以查阅教程:
Demonstration of Bulk2Space on demo1 dataset
如果您选择基于图像的数据(如MERFISH、SeqFISH或STARmap等)作为空间参考，请参考:
Demonstration of Bulk2Space on demo2 dataset
gitHub上给出了两个示例教程：
Integrating spatial gene expression and histomorphology in pancreatic ductal adenocarcinoma (PDAC) Spatially resolved analysis of mouse hypothalamus by Bulk2Space
Step1. 加载依赖包
import scanpy
import pandas as pd
import numpy as np
import matplotlib.colors as clr
import matplotlib.pyplot as plt
from scipy.stats import pearsonr
import bulk2space
from bulk2space import Bulk2Space
model = Bulk2Space()
这里我们检查一下示例数据的构成。示例数据包括Bulk数据+单细胞数据（表达矩阵+细胞类型信息）+空转数据（表达矩阵+细胞空间位置分布信息），共5个csv文件，具体如下：
pdac_bulk.csv bulk数据，单样本的一个表达矩阵
pdac_bulk = pd.read_csv('tutorial/data/pdac/pdac_bulk.csv')
print(pdac_bulk)
sc_data.csv 单细胞表达矩阵，单样本的一个表达矩阵，示例数据包含了1927个单细胞
sc_data = pd.read_csv('tutorial/data/pdac/sc_data.csv')
print(sc_data)
sc_meta.csv 单细胞表型数据，包含单细胞的细胞类型
sc_meta = pd.read_csv('tutorial/data/pdac/sc_meta.csv')
print(sc_meta)
st_data.csv 空转表达矩阵，示例数据包含429个spot。
st_data = pd.read_csv('tutorial/data/pdac/st_data.csv')
print(st_data)
st_data.csv 空转表型数据，包含spot的空间位置坐标。
st_meta = pd.read_csv('tutorial/data/pdac/st_meta.csv')
print(st_meta)
Step2. 将bulk-seq数据分解为scRNA-seq数据
generate_sc_meta, generate_sc_data = model.train_vae_and_generate(
input_bulk_path='tutorial/data/pdac/pdac_bulk.csv',
input_sc_data_path='tutorial/data/pdac/sc_data.csv',
input_sc_meta_path='tutorial/data/pdac/sc_meta.csv',
input_st_data_path='tutorial/data/pdac/st_data.csv',
input_st_meta_path='tutorial/data/pdac/st_meta.csv',
ratio_num=1,
top_marker_num=200,
gpu=0,
batch_size=512,
learning_rate=1e-4,
hidden_size=256,
epoch_num=3500,
vae_save_dir='tutorial/data/pdac/predata/save_model',
vae_save_name='pdac_vae',
generate_save_dir='tutorial/data/pdac/predata/output',
generate_save_name='pdac')
其实这个示例数据单细胞和空转数据并不大，但是我在运行的时候，这个步骤花了好久。。。
2.1 绘制生成的scRNA-seq数据的细胞类型比例
运行成功之后，生成的generate_sc_meta和generate_sc_data分别代表从bulk数据中反卷积获得的假定的单细胞表型数据（即细胞类型）和表达矩阵，可以可视化一下生成的细胞类型数量：
# The number of cells per cell type in deconvoluted bulk-seq data
ct_stat = pd.DataFrame(generate_sc_meta['Cell_type'].value_counts())
ct_name = list(ct_stat.index)
ct_num = list(ct_stat['Cell_type'])
color = ["#14B8A6", "#EF4444", '#F97316', "#7C3AED", '#64748B', '#3B82F6', "#0EA5E9",
'#F59E0B', '#10B981', '#EC4899', "#7DD3FC", "#F472B6", "#D946EF"]
plt.bar(ct_name, ct_num, color=color)
plt.xticks(ct_name, ct_name, rotation=90)
plt.title("The number of cells per cell type in bulk-seq data")
plt.xlabel("Cell type")
plt.ylabel("Cell number")
plt.show()
2.2 计算scRNA-seq参考数据和生成的scRNA-seq数据之间的标记基因表达相关性
此外，我们也可以比较真实的单细胞数据的细胞类型marker和bulk反卷积得到的单细胞细胞类型的marker的相似度：
# load input sc data
input_data = bulk2space.utils.load_data(
input_bulk_path='tutorial/data/pdac/pdac_bulk.csv',
input_sc_data_path='tutorial/data/pdac/sc_data.csv',
input_sc_meta_path='tutorial/data/pdac/sc_meta.csv',
input_st_data_path='tutorial/data/pdac/st_data.csv',
input_st_meta_path='tutorial/data/pdac/st_meta.csv'
)
# Calculate 200 marker genes of each cell type
sc = scanpy.AnnData(input_data['input_sc_data'].T)
sc.obs = input_data['input_sc_meta'][['Cell_type']]
scanpy.tl.rank_genes_groups(sc, 'Cell_type', method='wilcoxon')
marker_df = pd.DataFrame(sc.uns['rank_genes_groups']['names']).head(200)
marker = list(np.unique(np.ravel(np.array(marker_df))))
# the mean expression of 200 marker genes of input sc data
sc_marker = input_data['input_sc_data'].loc[marker, :].T
sc_marker['Cell_type'] = input_data['input_sc_meta']['Cell_type']
sc_marker_mean = sc_marker.groupby('Cell_type')[marker].mean()
# the mean expression of 200 marker genes of deconvoluted bulk-seq data
generate_sc_meta.index = list(generate_sc_meta['Cell'])
generate_sc_data_new = generate_sc_data.T
generate_sc_data_new['Cell_type'] = generate_sc_meta['Cell_type']
generate_sc_marker_mean = generate_sc_data_new.groupby('Cell_type')[marker].mean()
intersect_cell = list(set(sc_marker_mean.index).intersection(set(generate_sc_marker_mean.index)))
generate_sc_marker_mean= generate_sc_marker_mean.loc[intersect_cell]
sc_marker_mean= sc_marker_mean.loc[intersect_cell]
# calculate correlation
sc_marker_mean = sc_marker_mean.T
generate_sc_marker_mean = generate_sc_marker_mean.T
coeffmat = np.zeros((sc_marker_mean.shape[1], generate_sc_marker_mean.shape[1]))
for i in range(sc_marker_mean.shape[1]):
for j in range(generate_sc_marker_mean.shape[1]):
corrtest = pearsonr(sc_marker_mean[sc_marker_mean.columns[i]],
generate_sc_marker_mean[generate_sc_marker_mean.columns[j]])
coeffmat[i,j] = corrtest[0]
rf_ct = list(sc_marker_mean.columns)
generate_ct = list(generate_sc_marker_mean.columns)
# plot
fig, ax = plt.subplots()
im = ax.imshow(coeffmat, cmap='RdBu_r')
ax.set_xticks(np.arange(len(rf_ct)))
ax.set_xticklabels(rf_ct)
ax.set_yticks(np.arange(len(generate_ct)))
ax.set_yticklabels(generate_ct)
plt.xlabel("scRNA-seq reference")
plt.ylabel("deconvoluted bulk-seq")
plt.setp(ax.get_xticklabels(), rotation=90, ha="right", rotation_mode="anchor")
plt.colorbar(im)
ax.set_title("Expression correlation")
fig.tight_layout()
plt.show()
Step3. 将ST数据分解为空间分辨的单细胞转录组数据
df_meta, df_data = model.train_df_and_spatial_deconvolution(
generate_sc_meta,
generate_sc_data,
input_st_data_path='tutorial/data/pdac/st_data.csv',
input_st_meta_path='tutorial/data/pdac/st_meta.csv',
spot_num=500,
cell_num=10,
df_save_dir='tutorial/data/pdac/predata/save_model/',
df_save_name='pdac_df',
map_save_dir='tutorial/data/pdac/result',
map_save_name='pdac',
top_marker_num=200,
marker_used=True,
k=10)
利用深度学习模型，将空转数据分解。
Step4. 用单细胞分辨率绘制空间反褶积结果
我们首先研究了单细胞分辨率的空间反褶积结果。结果表明腺泡细胞、癌细胞和导管细胞具有明显的空间分布规律，与组织形态一致。
# Spatial mapping with single cell resolution
ct_type = list(df_meta['Cell_type'].unique())
color = ["#EF4444", "#10B981", '#14B8A6', "#7C3AED", '#3B82F6', '#64748B', "#D946EF",
'#F59E0B', '#F97316', '#7DD3FC', "#F472B6", "#EC4899"]
fig, ax = plt.subplots(figsize=(8,8))
for i in range(len(ct_type)):
ax.scatter(df_meta.loc[df_meta.Cell_type == ct_type[i], 'Cell_xcoord'],
df_meta.loc[df_meta.Cell_type == ct_type[i], 'Cell_ycoord'],
color = color[i], label = ct_type[i], s = 15)
plt.title("Spatial mapping with single cell resolution")
plt.xlabel("Cell_xcoord")
plt.ylabel("Cell_ycoord")
plt.legend(bbox_to_anchor=(1, 0.2), loc=3, borderaxespad=0, frameon=False, fontsize=15)
plt.show()
Step5. 用单细胞分辨率绘制标记基因的空间表达谱
标记基因的空间表达模式也与细胞类型的空间分布一致。
# gene expression
plot = df_meta
plot.index = list(plot['Cell'])
# 'REG1A': Acinar cells; 'CLDN1': Cancer clone A; 'KRT16': Cancer clone B; 'MUC5B': Ductal
gene = ['REG1A', 'CLDN1', 'KRT16', 'MUC5B']
plot[gene] = df_data.T[gene]
xcoord = np.array(plot['Cell_xcoord'])
ycoord = np.array(plot['Cell_ycoord'])
gene1 = np.array(plot['REG1A'])
gene2 = np.array(plot['CLDN1'])
gene3 = np.array(plot['KRT16'])
gene4 = np.array(plot['MUC5B'])
fig, axs = plt.subplots(1, 4, figsize=(12, 3))
a = axs[0].scatter(xcoord, ycoord, s=2.5, cmap='viridis', c=gene1)
fig.colorbar(a, ax=axs[0])
axs[0].set_title('REG1A')
b = axs[1].scatter(xcoord, ycoord, s=2.5, cmap='viridis', c=gene2)
fig.colorbar(b, ax=axs[1])
axs[1].set_title('CLDN1')
c = axs[2].scatter(xcoord, ycoord, s=2.5, cmap='viridis', c=gene3)
fig.colorbar(c, ax=axs[2])
axs[2].set_title('KRT16')
d = axs[3].scatter(xcoord, ycoord, s=2.5, cmap='viridis', c=gene4)
fig.colorbar(d, ax=axs[3])
axs[3].set_title('MUC5B')
plt.tight_layout()
plt.show()
Step6. 在spot分辨率下探讨反卷积结果
# spot-level
# calculate cell type proportion per spot
prop = df_meta[['Cell', 'Cell_type', 'Spot']].pivot_table(index=['Spot'], columns=['Cell_type'], aggfunc='count',values = 'Cell', fill_value=0)
prop = prop.div(prop.sum(axis=1), axis=0)
prop.columns = pd.Index(list(prop.columns))
prop['Spot_xcoord'] = np.array(df_meta.pivot_table(index=['Spot'])['Spot_xcoord'])
prop['Spot_ycoord'] = np.array(df_meta.pivot_table(index=['Spot'])['Spot_ycoord'])
# aggregate gene expression per spot
pred_spot_new = df_data.T
genes = pred_spot_new.columns
pred_spot_new['Spot'] = df_meta['Spot']
pred_spot_mean = pred_spot_new.groupby('Spot')[genes].mean()
6.1 绘制每个点的单元格类型比例
# plot cell type proportion
# Acinar cells; Cancer clone A; Cancer clone B; Ductal
spot_x = np.array(prop['Spot_xcoord'])
spot_y = np.array(prop['Spot_ycoord'])
cell1 = np.array(prop['Acinar cells'])
cell2 = np.array(prop['Cancer clone A'])
cell3 = np.array(prop['Cancer clone B'])
cell4 = np.array(prop['Ductal'])
fig, axs = plt.subplots(1, 4, figsize=(12, 3))
a = axs[0].scatter(spot_x, spot_y, s=10, cmap='viridis', c=cell1)
fig.colorbar(a, ax=axs[0])
axs[0].set_title('Acinar cells')
b=axs[1].scatter(spot_x, spot_y, s=10, cmap='viridis', c=cell2)
fig.colorbar(b, ax=axs[1])
axs[1].set_title('Cancer clone A')
c=axs[2].scatter(spot_x, spot_y, s=10, cmap='viridis', c=cell3)
fig.colorbar(c, ax=axs[2])
axs[2].set_title('Cancer clone B')
d=axs[3].scatter(spot_x, spot_y, s=10, cmap='viridis', c=cell4)
fig.colorbar(d, ax=axs[3])
axs[3].set_title('Ductal')
plt.tight_layout()
plt.show()
6.2 绘制每个spot的基因表达量
# plot gene expression proportion
# 'REG1A': Acinar cells; 'CLDN1': Cancer clone A; 'KRT16': Cancer clone B; 'MUC5B': Ductal
spot_x = np.array(prop['Spot_xcoord'])
spot_y = np.array(prop['Spot_ycoord'])
g1 = np.array(pred_spot_mean['REG1A'])
g2 = np.array(pred_spot_mean['CLDN1'])
g3 = np.array(pred_spot_mean['KRT16'])
g4 = np.array(pred_spot_mean['MUC5B'])
fig, axs = plt.subplots(1, 4, figsize=(12, 3))
a = axs[0].scatter(spot_x, spot_y, s=10, cmap='viridis', c=g1)
fig.colorbar(a, ax=axs[0])
axs[0].set_title('REG1A')
b=axs[1].scatter(spot_x, spot_y, s=10, cmap='viridis', c=g2)
fig.colorbar(b, ax=axs[1])
axs[1].set_title('CLDN1')
c=axs[2].scatter(spot_x, spot_y, s=10, cmap='viridis', c=g3)
fig.colorbar(c, ax=axs[2])
axs[2].set_title('KRT16')
d=axs[3].scatter(spot_x, spot_y, s=10, cmap='viridis', c=g4)
fig.colorbar(d, ax=axs[3])
axs[3].set_title('MUC5B')
plt.tight_layout()
plt.show()
三. 写在最后
关于Bulk2Space算法的实际用途，我唯一想到的是，针对稀有样本如果只有Bulk数据，同时存在单细胞分辨率下的数据，且存在空转级别的数据，然后利用这个算法去推断bulk数据，扩大单细胞/空转的样本量。此外，还可以用这个算法，基于bulk数据，来验证单细胞和空转的结果（有一点鸡生蛋，蛋生鸡的悖论感。。。）。除此之外，大家觉得这个算法还能讲一些什么故事呢？欢迎在文末进行留言。
- END -
https://mp.weixin.qq.com/s/GZAoS1ACW6IsN8kuAVZzew