数据来源
单细胞数据来自2023年的Cell人胎脑，"Spatiotemporal transcriptome atlas reveals the regional specification of the developing human brain"，选取了其中14PCW的部分进行迁移整合；
空间转录组数据来自华大Stereo-seq测序的下机GEM/gef原始文件，选取了其中一张1*1的小芯片进行测试。所有的数据都需要转成AnnData格式，即适用于scanpy分析的h5ad格式。格式转换方法可以参考Stereopy官方教程(Input & Output - Stereopy)，这里不再详细说明。
环境配置
首先需要配置好环境，包的下载来自GitHub - mexchy1000/CellDART: domain adaptation of spatial and single-cell transcriptome
pip install git+https://github.com/mexchy1000/CellDART.gitGithub项目里也有相应的教程，但直接套用Stereo-seq数据会出现一些问题，但可以作为参考
CellDART/CellDART_example_mousebrain_markers.ipynb at master · mexchy1000/CellDART · GitHub
数据和包导入
配置好环境后具体如下，首先导入需要的包
import scanpy as sc
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from CellDART import da_cellfraction
from CellDART.utils import random_mix
from sklearn.manifold import TSNE分别导入单细胞数据adata和空间转录组数据adata_spatial
adata = sc.read('/data/work/Cell_fetal_brain_singleCell/14.h5ad')
adataAnnData object with n_obs × n_vars = 35175 × 36601
obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'percent.mt', 'region', 'week', 'celltype'
var: 'vst.mean', 'vst.variance', 'vst.variance.expected', 'vst.variance.standardized', 'vst.variable'
adata_spatial = sc.read('/data/work/scanpy/D03657A1_scanpy_out.h5ad')
adata_spatialAnnData object with n_obs × n_vars = 116233 × 47627
obs: 'total_counts', 'n_genes_by_counts', 'pct_counts_mt', 'orig.ident', 'x', 'y'
var: 'n_cells', 'n_counts', 'mean_umi'
uns: 'bin_size', 'bin_type', 'resolution', 'sn'
obsm: 'spatial'
简单查看一下数据的基本结构。这里文献的单细胞数据没有进行太多的预处理，很多时候可以和文章作者申请中间数据获得，这里由于空间组数据的采样时间，选择性地提取了第14孕周的数据，即文中的14.h5ad。空间转录组数据使用stereo包（地址见上文）进行格式转换，输出为适合scanpy分析的AnnData格式，即*_scanpy_out.h5ad。
数据预处理
下面对单细胞数据进行简单的预处理
#Preprocessing
adata.var['mt'] = adata.var_names.str.startswith('Mt-') # annotate the group of mitochondrial genes as 'mt'
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
sc.pp.normalize_total(adata)之后对单细胞数据进行降维聚类，文献中的聚类结果储存在obs的'celltype'列中
#PCA and clustering : Known markers with 'celltype'
sc.tl.pca(adata, svd_solver='arpack')
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution = 1)
sc.pl.umap(adata, color=['leiden','celltype'])
可视化了单细胞数据自聚类（leiden）和文献中的聚类结果。
一个小问题
可以查看一下聚类结果对应的具体数目
adata.obs['celltype'].value_counts()GABA_Neuron 8391 NPC 7646 Gluta_Neuron 5098 Microglia 4308 Cere_GC 3462 Endothelial 1847 UBC/eCN 1402 Pericyte 1303 OPC 1010 Astrocyte 402 VLMC 116 neutrophils 73 NK 60 Erythrocyte 33 Neuroectoderm 18 Purkinje 5 NeuralCrest 1 Name: celltype, dtype: int64
这里发现一个问题，名称为‘NeuralCrest’的类只包含1个细胞，肯定无法进行后续的rank_gene分析，这里我为了保证数据完整的同时顺利进行后续分析，选择了将这里的1个细胞合并进一个神经大类，使其不影响其他大类的迁移整合。方法如下：
groups_to_merge = {'NeuralCrest': 'NPC'}
adata.obs['celltype'] = adata.obs['celltype'].replace(groups_to_merge)即将NeuralCrest类的一个值合并进NPC类。
特征基因提取
紧接着进行后续高表达基因的选取（这样就不会因为那个类报错了）
sc.tl.rank_genes_groups(adata, 'celltype', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=20, sharey=False)
将选取的高表达基因转换成数据框格式检查一下前5行
genelists=adata.uns['rank_genes_groups']['names']
df_genelists = pd.DataFrame.from_records(genelists)
df_genelists.head(5)
创建一个新的数据集res_gene存放单细胞的特征基因
num_markers=10
res_genes = []
for column in df_genelists.head(num_markers):
res_genes.extend(df_genelists.head(num_markers)[column].tolist())
res_genes_ = list(set(res_genes))紧接着可以在我们的空间转录组数据中选取相同的特征基因，并用特征基因替换现有的基因集，用来进行后续的迁移训练
adata_spatial.var_names_make_unique()
inter_genes = [val for val in res_genes_ if val in adata_spatial.var.index]
print('Selected Feature Gene number',len(inter_genes))
adata = adata[:,inter_genes]
adata_spatial = adata_spatial[:,inter_genes]格式转换
正式进行训练前还需要将数据格式转成数组和字典
mat_sc = adata.X.toarray()
mat_sp = adata_spatial.X.todense()
df_sc = adata.obs
lab_sc_sub = df_sc.celltype
sc_sub_dict = dict(zip(range(len(set(lab_sc_sub))), set(lab_sc_sub)))
sc_sub_dict2 = dict((y,x) for x,y in sc_sub_dict.items())
lab_sc_num = [sc_sub_dict2[ii] for ii in lab_sc_sub]
lab_sc_num = np.asarray(lab_sc_num, dtype='int')又一个小问题
这里我在尝试用官网代码转换时，默认将我的数据集转成了稀疏矩阵
mat_sc = adata_cortex.X #adata_cortex是单细胞数据
mat_sp = adata_spatial_anterior.X.todense() #adata_spatial_anterior是空间转录组数据稀疏矩阵是一种优化的数据结构，用于表示大多数元素为零的矩阵（这里猜测可能是我们选取的数据集体量比较小）。在转录组数据中，通常有很多基因在很多细胞中的表达量为零。如果使用密集数组来存储这些数据，会占据大量的内存空间。为了节省内存空间，稀疏矩阵只存储非零元素的数据和它们的索引。
但是在后续random_mix()函数中需要通过len()统计矩阵的大小，而稀疏矩阵不能通过len()，只能通过getnnz() 或者 shape[0]获得。所以这里我们需要手动将稀疏矩阵转换成密集矩阵，即上文给出的两行代码。这样我们使用len()验证一下
len(mat_sc)35175
是可以顺利输出矩阵大小的
生成随机混合训练集，选取的样本数为5000
sc_mix, lab_mix = random_mix(mat_sc, lab_sc_num, nmix=5, n_samples=5000)
def log_minmaxscale(arr):
arrd = len(arr)
arr = np.log1p(arr)
return (arr-np.reshape(np.min(arr,axis=0), (1, arr.shape[1])))/np.reshape((np.max(arr, axis=0)-np.min(arr,axis=0)),(1, arr.shape[1]))
sc_mix_s = log_minmaxscale(sc_mix)
mat_sp_s = log_minmaxscale(mat_sp)
mat_sc_s = log_minmaxscale(mat_sc)打印看一下训练集的大小
print(mat_sp_s.shape, mat_sc_s.shape, sc_mix_s.shape)(116233, 148) (35175, 148) (5000, 148)
训练过程
利用上文生成的训练集训练模型
embs, clssmodel = da_cellfraction.train(sc_mix_s, lab_mix, mat_sp_s,
alpha=1, alpha_lr=5, emb_dim = 64, batch_size = 512,
n_iterations = 2000,
initial_train=True,
initial_train_epochs=10)参数含义参考官方教程
预测过程
用训练好的模型预测我们的空间转录组数据，打印上文生成的单细胞聚类结果字典
pred_sp = clssmodel.predict(mat_sp_s)
print(sc_sub_dict){0: 'Gluta_Neuron', 1: 'GABA_Neuron', 2: 'UBC/eCN', 3: 'Astrocyte', 4: 'Pericyte', 5: 'Endothelial', 6: 'Purkinje', 7: 'NPC', 8: 'VLMC', 9: 'neutrophils', 10: 'Neuroectoderm', 11: 'NK', 12: 'Microglia', 13: 'Erythrocyte', 14: 'Cere_GC', 15: 'OPC'}
将预测结果放在空间转录组的一个obs中 'Pred_label'可视化，numlist中放的是想要看的聚类簇，对应关系参考上一条的字典
def plot_cellfraction(visnum):
adata_spatial.obs['Pred_label'] = pred_sp[:,visnum]
sc.pl.spatial(
adata_spatial,
img_key="hires",
color='Pred_label',
palette='Set1',
spot_size=50,
legend_loc=None,
title = sc_sub_dict[visnum])
numlist = [0,1,7]
for num in numlist:
plot_cellfraction(num)又又一个小问题
这里官网plot_cellfraction()函数中用的是size参数而非spot_size，和我一样报错而百思不得其解的友友可以试着替换一下这里（抹泪.jpg）以及仍是因为数据本身原因，我需要将spot_size调得很大才能看见图（指大家多试试）（再次抹泪）
最后终于能看见一个比较OK的可视化结果，如下
当然后续还可以进行很多其他的分析，记得保存结果，想要单独可视化可以参考
adata_spatial.write_h5ad("/data/work/CellDART/D03657A1_sc2sp.h5ad")再次导入因为前期训练中间文件无了，所以只能利用字典单独调用结果
import scanpy as sc
adata_spatial = sc.read('/data/work/CellDART/D03657A1_sc2sp.h5ad')
sc_sub_dict = {0: 'Gluta_Neuron', 1: 'GABA_Neuron', 2: 'UBC/eCN', 3: 'Astrocyte', 4: 'Pericyte', 5: 'Endothelial', 6: 'Purkinje', 7: 'NPC', 8: 'VLMC', 9: 'neutrophils', 10: 'Neuroectoderm', 11: 'NK', 12: 'Microglia', 13: 'Erythrocyte', 14: 'Cere_GC', 15: 'OPC'}
sc.pl.spatial(
adata_spatial,
img_key="hires",
color='Pred_label',
palette='Set1',
spot_size=50,
legend_loc=None,
title = sc_sub_dict[0])祝大家好运！
746
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
