Visium HD数据分析之Bin2Cell |生信开发实战

Python                  # In[1]: 导入模块                                   import matplotlib.pyplot as plt                  import scanpy as sc                  import numpy as np                  import os                  import bin2cell as b2c                                   # In[2]:地址设置                                   os.chdir("../2_bin2cell/Demo2/Output")                  os.makedirs("stardist", exist_ok=True)                  path = "../Input/binned_outputs/square_002um/"                  source_image_path = "../Input/Visium_HD_Mouse_Brain_tissue_image.tif"                  spaceranger_image_path = "../Input/spatial"      

Python# In[3]:导入HD表达数据（2um）和高清HE图片
adata = b2c.read_visium(path,                         source_image_path = source_image_path,                         spaceranger_image_path = spaceranger_image_path                       )adata.var_names_make_unique()adata
# In[4]:过滤低质量bin
sc.pp.filter_genes(adata, min_cells=3)sc.pp.filter_cells(adata, min_counts=1)adata
# In[5]:基于高清图像保存低分辨率图像并存入adata.obsm['spatial_cropped']
mpp = 0.5 # 每像素微米数
b2c.scaled_he_image(adata,                     mpp=mpp,                     save_path="stardist/he.tiff")  adata      

Python                  # In[6]: 校正bin的行列偏移                                   b2c.destripe(adata)                  adata                                   # 取子图像进行可视化                  mask = ((adata.obs['array_row'] >= 1450) &                          (adata.obs['array_row'] <= 1550) &                          (adata.obs['array_col'] >= 250) &                          (adata.obs['array_col'] <= 450)                         )                                           # In[7]: Picture1 使用adata.obsm['spatial_cropped']进行可视化                                   bdata = adata[mask]                  sc.pl.spatial(bdata, color=[None, "n_counts", "n_counts_adjusted"], img_key="0.5_mpp", basis="spatial_cropped")                  plt.savefig('destripe_mask.png', dpi=300, bbox_inches='tight')                  plt.close()                                   # In[8]: Picture2 使用默认图片进行可视化                                   sc.pl.spatial(bdata, color=[None, "n_counts", "n_counts_adjusted"])                  plt.savefig('destripe_mask_default.png', dpi=300, bbox_inches='tight')                  plt.close()

Python                  # In[9]: 细胞核分割（约13min）                                   # block_size StarDist 在 predict_instances_big() 函数中使用的输入参数，指定了作为单个瓦片处理的图像的正方形边长(4096)                  # min_overlap StarDist 在 predict_instances_big() 函数中使用的输入参数，指定了相邻瓦片之间的最小重叠区域(128)                  # context StarDist 在 predict_instances_big() 函数中使用的输入参数，指定了在每个瓦片的边界上要丢弃的图像上下文量                  # prob_thresh控制调用对象的严格性，确定像素被分类为前景（细胞）的概率，值越小分割的细胞越多                  # nms_thresh 非最大抑制中定义了两个或多个预测对象（putative objects）被视为重叠并因此被合并成一个单一标签的阈值                  # labels_npz_path代表输出的npz文件名（稀疏矩阵文件）                  # stardist_model选择StarDist 的 H&E 模型                  # nms_thresh                  b2c.stardist(image_path="stardist/he.tiff",                               labels_npz_path="stardist/he.npz",                               stardist_model="2D_versatile_he",                               prob_thresh=0.01                              )                                   # In[10]:转录本分配                                   b2c.insert_labels(adata,                                    labels_npz_path="stardist/he.npz",                                    basis="spatial",                                    spatial_key="spatial_cropped",                                    mpp=mpp,                                    labels_key="labels_he" # .obs key to store the labels under.                                   )                  adata                  # adata.uns多了'bin2cell'                  # adata.obs 多了'labels_he'                                   # In[11]: Picture3 使用adata.obsm['spatial_cropped']进行可视化nuclei                                   bdata = adata[mask]                                   #0 means unassigned                  bdata = bdata[bdata.obs['labels_he']>0]                  bdata.obs['labels_he'] = bdata.obs['labels_he'].astype(str)                                   sc.pl.spatial(bdata, color=[None, "labels_he"], img_key="0.5_mpp", basis="spatial_cropped")                  # sc.pl.spatial(bdata, color=[None, "labels_he"], img_key="0.5_mpp", basis="spatial_cropped",save='_he_labels.pdf')                  plt.savefig('he_nuclei_mask_test.png', dpi=300, bbox_inches='tight')                  plt.close()                                   # In[12]: Picture4 使用HE图片进行可视化nuclei                                   #the label viewer wants a crop of the processed image                  #get the corresponding coordinates spanning the subset object                  # 获取spatial_cropped的coordinate                  crop = b2c.get_crop(bdata,                                      basis="spatial",                                      spatial_key="spatial_cropped",                                      mpp=mpp)                                   #if this errors about missing get_cmap(), downgrade matplotlib below 3.9.0                  rendered = b2c.view_stardist_labels(image_path="stardist/he.tiff",                                                      labels_npz_path="stardist/he.npz",                                                      crop=crop                                                     )                  plt.imshow(rendered)                  plt.savefig('he_nuclei_mask2.png', dpi=300, bbox_inches='tight')                  # plt.savefig('he_labels_image.pdf')                  plt.close()                  

Python                  # In[13]: 扩展nuclei为cell（4 min）                                   # 如果一个 bin 与多个细胞核的距离相等，则根据其基因表达谱的相似性将其分配给一个 bin                  # max_bin_distance指定扩展核标签时最大的 bin 数量 (2)                  # volume_ratio如果 max_bin_distance 设置为 None，则根据ceil((volume_ratio**(1/3)-1) * sqrt(n_bins/pi))进行扩展（4）                  # k为每个未分配的 bin 寻找潜在最近邻的已分配空间坐标 bin 的数量 （4）                  # subset_pca 如果设置为 True，则只对于那些与多个细胞核的bins距离相等的未分配bins进行PCA表示和决断（True）                                   b2c.expand_labels(adata,                                    labels_key='labels_he',                                    expanded_labels_key="labels_he_expanded"                                   )                  # adata.obs 多了'labels_he_expanded'                                   # In[14]:Picture4 使用adata.obsm['spatial_cropped']进行可视化cell                                   bdata = adata[mask]                  bdata = bdata[bdata.obs['labels_he_expanded']>0]                  bdata.obs['labels_he_expanded'] = bdata.obs['labels_he_expanded'].astype(str)                                   sc.pl.spatial(bdata, color=[None, "labels_he_expanded"], img_key="0.5_mpp", basis="spatial_cropped")                  # sc.pl.spatial(bdata, color=[None, "labels_he_expanded"], img_key="0.5_mpp", basis="spatial_cropped")                  plt.savefig('he_cell_mask.png', dpi=300, bbox_inches='tight')                  plt.close()

Python                  # In[15]: 根据表达矩阵进行细胞分割                                   # 基于adata.obs['n_counts_adjusted']生成高斯模糊图像                  # log1p是否做log1p-transform（False）                  # sigma使用 skimage.filters.gaussian() 对最终图像进行高斯模糊处理，sigma 值决定模糊的程度(None)                  b2c.grid_image(adata, "n_counts_adjusted", mpp=mpp, sigma=5, save_path="stardist/gex.tiff")                                   # In[16]: 使用StarDist's fluorescence model进行细胞分割（20min）                                   # 由于网络原因会下载模型失败，多次尝试即可                  b2c.stardist(image_path="stardist/gex.tiff",                               labels_npz_path="stardist/gex.npz",                               stardist_model='2D_versatile_fluo',                               prob_thresh=0.05,                               nms_thresh=0.5                              )                                   # In[17]: 转录本分配                                   b2c.insert_labels(adata,                                    labels_npz_path="stardist/gex.npz",                                    basis="array",                                    mpp=mpp,                                    labels_key="labels_gex"                                   )                  adata                  # adata.obs 多了'labels_gex'                                   # In[18]: 使用adata.obsm['spatial_cropped']进行可视化cell                                   bdata = adata[mask]                  bdata = bdata[bdata.obs['labels_gex']>0]                  bdata.obs['labels_gex'] = bdata.obs['labels_gex'].astype(str)                                   sc.pl.spatial(bdata, color=[None, "labels_gex"], img_key="0.5_mpp", basis="spatial_cropped")                  plt.savefig('gex_cell_mask.png', dpi=300, bbox_inches='tight')                  plt.close()                                   # In[19]:使用gex.tiff进行可视化cell                                   #the label viewer wants a crop of the processed image                  #get the corresponding coordinates spanning the subset object                  crop = b2c.get_crop(bdata, basis="array", mpp=mpp)                                   #if this errors about missing get_cmap(), downgrade matplotlib below 3.9.0                  rendered = b2c.view_stardist_labels(image_path="stardist/gex.tiff",                                                      labels_npz_path="stardist/gex.npz",                                                      crop=crop                                                     )                  plt.imshow(rendered)                  plt.savefig('gex_cell_mask2.png', dpi=300, bbox_inches='tight')                  plt.close()    

Python                  # In[20]: 用GEX的结果填补HE结果中未识别的细胞                                   b2c.salvage_secondary_labels(adata,                                               primary_label="labels_he_expanded",                                               secondary_label="labels_gex",                                               labels_key="labels_joint"                                              )                  adata                  # adata.obs 多了'labels_joint', 'labels_joint_source'                                   # In[21]:使用adata.obsm['spatial_cropped']进行可视化cell                                   bdata = adata[mask]                  bdata = bdata[bdata.obs['labels_joint']>0]                  bdata.obs['labels_joint'] = bdata.obs['labels_joint'].astype(str)                                   sc.pl.spatial(bdata, color=[None, "labels_joint_source", "labels_joint"], img_key="0.5_mpp", basis="spatial_cropped")                  # sc.pl.spatial(bdata, color=[None, "labels_joint_source", "labels_joint"], img_key="0.5_mpp", basis="spatial_cropped",save='_salvage_labels.svg')                  plt.savefig('he_gex_mask_bin.png', dpi=300, bbox_inches='tight')                  plt.close()                                   # In[22]:将同一个cell中的bin进行合并                                   cdata = b2c.bin_to_cell(adata, labels_key="labels_joint", spatial_keys=["spatial", "spatial_cropped"])                  cdata                  # AnnData object with n_obs × n_vars = 69902 × 18823                  #     obs: 'bin_count', 'array_row', 'array_col', 'labels_joint_source'                  #     var: 'gene_ids', 'feature_types', 'genome', 'n_cells'                  #     uns: 'spatial'                  #     obsm: 'spatial', 'spatial_cropped'                                   # In[23]:使用adata.obsm['spatial_cropped']进行可视化cell                                   cell_mask = ((cdata.obs['array_row'] >= 1450) &                               (cdata.obs['array_row'] <= 1550) &                               (cdata.obs['array_col'] >= 250) &                               (cdata.obs['array_col'] <= 450)                              )                                   ddata = cdata[cell_mask]                  sc.pl.spatial(ddata, color=["bin_count", "labels_joint_source"], img_key="0.5_mpp", basis="spatial_cropped")                  # sc.pl.spatial(ddata, color=["bin_count", "labels_joint_source"], img_key="0.5_mpp", basis="spatial_cropped",s=4)                  plt.savefig('he_gex_mask_cell.png', dpi=300, bbox_inches='tight')                  plt.close()                                   # In[24]:保存数据                  cdata.write_h5ad('m_brain_b2c.h5ad')                  adata.write_h5ad('m_brain_2um.h5ad')

Python                  # In[1]: 导入模块                  

                  import scanpy as sc                  import os                  import bin2cell as b2c                  import celltypist                  from celltypist import models                  import numpy as np                  from matplotlib import rcParams                  from matplotlib import font_manager                  import matplotlib.pyplot as plt                  

                  # In[2]:地址设置                  

                  os.chdir("../2_bin2cell/Demo2/Output")                  os.makedirs("scanpy", exist_ok=True)                  os.makedirs("spatial", exist_ok=True)                  

                  path008 = "../Input/binned_outputs/square_008um/"                  path002 = "../Input/binned_outputs/square_002um/"                  source_image_path = "../Input/Visium_HD_Mouse_Brain_tissue_image.tif"                  prediction_model = "../CellTypist/models/Mouse_Whole_Brain.pkl"

Python                  # In[3]:导入HD表达数据（8um）                                    bdata = b2c.read_visium(path008, source_image_path = source_image_path)                  bdata.var_names_make_unique()                  bdata                                    # In[4]:                                   bdata.raw = bdata                                    # In[5]:过滤低质量bin                                   sc.pp.filter_genes(bdata, min_cells=3)                  sc.pp.filter_cells(bdata, min_genes=100)                  sc.pp.calculate_qc_metrics(bdata,inplace=True)                  bdata                                    # In[6]:寻找高变基因、标准化                  
                 sc.pp.highly_variable_genes(bdata,n_top_genes=5000,flavor="seurat_v3")                  sc.pp.normalize_total(bdata,target_sum=1e4)                  sc.pp.log1p(bdata)                                    # In[7]:基于celltypist模型进行细胞类型鉴定                  
                predictions_8bin = celltypist.annotate(bdata, model = prediction_model, majority_voting = False)                  bdata = predictions_8bin.to_adata()                  bdata                                    # In[8]:使用高变基因进行降维和可视化                                   bdata = bdata[:, bdata.var["highly_variable"]].copy()                  sc.pp.scale(bdata, max_value=10)                  sc.pp.pca(bdata, use_highly_variable=True)                  sc.pp.neighbors(bdata)                  sc.tl.umap(bdata)                                    # In[9]:表达量小提琴图                  # 每细胞基因数、log1p(x) = log(1 + x)、每细胞UMI数、log1p(x) = log(1 + x)、细胞中表达量最高的前50个基因、细胞中表达量最高的前100个基因                  sc.set_figure_params(dpi=50,fontsize=10,)                  sc.pl.violin(bdata, ['n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes'],                               jitter=0.1, multi_panel=True)                  plt.savefig('scanpy/8um_violin.png', dpi=300, bbox_inches='tight')                  plt.close()                                                # In[10]: UMAP                                    sc.set_figure_params(dpi=100,fontsize=10,)                  sc.tl.leiden(bdata,resolution=6,key_added='leiden')                  del bdata.uns["leiden_colors"]                  sc.pl.umap(bdata,color=['leiden'],size=2,wspace=0.25,frameon=False)                  plt.savefig('scanpy/8um_umap.png', dpi=300, bbox_inches='tight')                  plt.close()                                   # In[11]: save while recovering raw                                    bdata.raw.to_adata().write('mb_8um.h5ad')

Python                  # In[12]:导入HD细胞分割结果（2um）                  
                  cdata = sc.read_h5ad('../2_bin2cell/Demo2/Output/m_brain_b2c.h5ad')                  cdata                  # AnnData object with n_obs × n_vars = 69902 × 18823                  #     obs: 'bin_count', 'array_row', 'array_col', 'labels_joint_source'                  #     var: 'gene_ids', 'feature_types', 'genome', 'n_cells'                  #     uns: 'spatial'                  #     obsm: 'spatial', 'spatial_cropped'                                    # In[13]:过滤过小的细胞                                    cdata = cdata[cdata.obs['bin_count']>7] # min 8 bins                  cdata                  # need integers for seuratv3 hvgs                  cdata.X.data = np.round(cdata.X.data)                  cdata.raw = cdata.copy()                                    # In[14]:过滤低质量基因和细胞                                    sc.pp.filter_genes(cdata, min_cells=3)                  sc.pp.filter_cells(cdata,min_genes=100)                  sc.pp.calculate_qc_metrics(cdata,inplace=True)                  cdata                  # AnnData object with n_obs × n_vars = 60274 × 18552                  #     obs: 'bin_count', 'array_row', 'array_col', 'labels_joint_source', 'n_genes', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes'                  #     var: 'gene_ids', 'feature_types', 'genome', 'n_cells', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts'                  #     uns: 'spatial'                  #     obsm: 'spatial', 'spatial_cropped'                                    # In[15]:寻找高变基因、标准化                                    sc.pp.highly_variable_genes(cdata,n_top_genes=5000,flavor="seurat_v3")                  sc.pp.normalize_total(cdata,target_sum=1e4)                  sc.pp.log1p(cdata)                                    # In[16]:基于celltypist模型进行细胞类型鉴定                                    b2c_predictions = celltypist.annotate(cdata, model = prediction_model, majority_voting = False)                  cdata = b2c_predictions.to_adata()                  cdata                                    # In[17]:使用高变基因进行降维和可视化                                    cdata = cdata[:, cdata.var["highly_variable"]]                  sc.pp.scale(cdata, max_value=10)                  sc.pp.pca(cdata, use_highly_variable=True)                  sc.pp.neighbors(cdata)                  sc.tl.umap(cdata)                                    # In[18]:表达量小提琴图                                    sc.set_figure_params(dpi=50,fontsize=10,)                  sc.pl.violin(cdata, ['n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes','bin_count'],                               jitter=0.05, multi_panel=True)                  plt.savefig('scanpy/cell_violin.png', dpi=300, bbox_inches='tight')                  plt.close()                                    # In[19]:UMAP                  sc.set_figure_params(dpi=100,fontsize=10,)                  sc.tl.leiden(cdata,resolution=6,key_added='leiden')                  sc.pl.umap(cdata,color=['leiden'],size=2,wspace=0.25,frameon=False)                  plt.savefig('scanpy/cell_umap.png', dpi=300, bbox_inches='tight')                  plt.close()                                    # In[20]:save while recovering raw                                    cdata.raw.to_adata().write('mb_b2c.h5ad')