10X空间转录组Visium && HD && Xenium分析
2024-06-02 本文已影响0人
单细胞空间交响乐
作者，Evil Genius
哎呀， Xenium和HD的空间项目进来了，又要有很多新的工作要做了。
很多分析呢，大家如果只是自己玩一玩，那随便做做就也可以了，但如果像我一样作为公司的分析人员，尤其是大公司的核心分析人员，生信经理之类的，那一定要对自己的分析结果负责，对要求就会高得多,因为分析的水平代表了公司水平。
还有很多人我觉得态度上就有问题，认为官方代码都是现成的，完全不需要公司分析，这种偏见从我做生信就一直有，其实就跟我们上学一样，上高中大家用的一样的教材，都是看的一样的内容，但是高考的结果却是千差万别。cell、nature等高分杂志的文章代码很多都是公开的，但是大家用这些代码能发cell、nature么？再比如单细胞平台，10X单细胞的平台原理大家都知道，但是国内做的水平就是比国外要差，为什么？
其实我们很多都缺乏一种精神，这种精神是在西方世界、包括日本很推崇并且深入骨髓的，那就是工匠精神。
关于Xenium，主推squidpy，文章在高精度空间转录组分析之squidpy和空间网络图,主要分析就是降维聚类和邻域富集分析。
Xenium中，Seurat也出了示例教程，在Analysis of Image-based Spatial Data in Seurat
关于HD ，之前分享过STEP，文章在10Xvisium HD高精度平台探索, 主要分析精度为8um 和 16um
高精度HD目前Seurat也更新的教程，在Analysis, visualization, and integration of Visium HD spatial datasets with Seurat • Seurat (satijalab.org)
总体来讲，python版本的分析更加优秀，这一篇对visium、HD、Xenium
基础分析进行全系列的脚本更新(因为这是公司的要求)。
我们首先来看visium
import CellScopes as cs
sham = cs.read_visium("/home/data/Jia/shamA_visium/MGI0805_A3_147_msham/outs/")
Normalization, dimension reduction and cell clustering
sham = cs.normalize_object(sham; scale_factor = 10000)
sham = cs.scale_object(sham)
sham = cs.find_variable_genes(sham; nFeatures = 1000)
sham = cs.run_pca(sham; method=:svd, pratio = 1, maxoutdim = 10)
sham = cs.run_umap(sham; min_dist=0.2, n_neighbors=20)
sham = cs.run_clustering(sham; res=0.01, n_neighbors=20)
可视化
cs.dim_plot(sham; dim_type ="umap", marker_size=8)
cs.sp_dim_plot(sham, "cluster";
marker_size = 8, width=600, height=500,
do_label=false)
cs.sp_dim_plot(sham, "cluster";
marker_size = 8, width=600, height=500,
do_label=false, alpha=0.3)
图片剪切
cs.sp_dim_plot(sham, "cluster"; do_label = false, do_legend = true, img_res = "low",
marker_size = 8, width=600, height=500, adjust_contrast = 1,
adjust_brightness= 0.1, alpha =0.4, x_lims=(200, 440),
y_lims=(200, 400))
基因可视化
cs.sp_feature_plot(sham, ["Aqp2","Slc7a13", "Umod", "Slc12a1","Slc12a3","Slc34a1"];
marker_size = 8, color_keys=["gray90", "lemonchiffon" ,"red"],
adjust_contrast=1, adjust_brightness = 0.3, alpha=1)
Plot genes on the selected region
cs.sp_feature_plot(sham, ["Aqp2","Slc7a13", "Umod", "Slc12a1","Slc12a3","Slc34a1"];
marker_size = 8, color_keys=["gray90", "lemonchiffon" ,"red"],
adjust_contrast=1, adjust_brightness = 0.3, alpha=1, x_lims=(200, 440), y_lims=(200, 400))
接下来我们来看HD，示例数据在 https://www.10xgenomics.com/datasets/visium-hd-cytassist-gene-expression-libraries-of-mouse-brain-he，其中我们需要注意的是HD可以有多个精度，Data from each bin size is stored in separate layers, which allow users to select data at different resolutions to analyze.，多个精度的分析结果我们都要保留(2um、4um、8um、16um精度)。
import CellScopes as cs
hd_dir = "/mnt/sdb/visiumHD/hd_output/outs/"
hd = cs.read_visiumHD(hd_dir)
1. loading 2um binned data...
Formatting cell polygons...
Progress: 100%[==================================================] Time: 0:00:5039m
1. 2um binned data loaded!
2. loading 8um binned data...
Formatting cell polygons...
Progress: 100%[==================================================] Time: 0:00:08
2. 8um binned data loaded!
3. loading 16um binned data...
Formatting cell polygons...
3. 16um binned data loaded!
VisiumHDObject in CellScopes.jl
Genes x Cells = 18988 x 397258
Available data:
- layerData
- rawCount
- normCount
- scaleCount
- metaData
- spmetaData
- varGene
- dimReduction
- clustData
- imageData
- alterImgData
- polygonData
- defaultData
All fields:
- layerData
- rawCount
- normCount
- scaleCount
- metaData
- spmetaData
- varGene
- dimReduction
- clustData
- imageData
- alterImgData
- polygonData
- defaultData
Layer selection选择，我们需要选择需要分析的精度，当然，也可以是多个精度
hd = cs.set_default_layer(hd; layer_slot = "2_um")
## or
hd = cs.set_default_layer(hd; layer_slot = "16_um")
(Optional) Clustering and dimensional reduction,这一步是可选的，因为Space Ranger已经为8微米和16微米的桶尺寸提供了cluster和UMAP结果。
从8µm bin大小重新聚集细胞可能会很耗时，因为spot数量通常超过300,000。
可视化
hd = cs.set_default_layer(hd; layer_slot = "8_um")
cs.sp_dim_plot(hd, "cluster"; width=1300, height=1000,
do_legend=true, legend_size=40, legend_fontsize=40)
要使方向与原始组织方向对齐，可以使用convert_image_data函数。下面是调整方向的方法:
@time hd = cs.convert_image_data(hd)
转换图像数据后，所有图形将使用调整后的坐标进行可视化。
hd = cs.set_default_layer(hd; layer_slot = "8_um")
cs.sp_dim_plot(hd, "cluster"; width=1300, height=1000,
do_legend=true, legend_size=40, legend_fontsize=40)
还可以在同一图上突出一个特定的cluster
hd = cs.set_default_layer(hd; layer_slot = "16_um")
hd = cs.convert_image_data(hd; layer_slot = "16_um")
cs.sp_dim_plot(hd, "cluster"; width=1000, height=1200,
anno_color = Dict("13" => "green2") ,do_legend=true, img_res = "high", stroke_width=0.2,
cell_highlight = "13",legend_size=30, adjust_contrast =1, adjust_brightness=0.0, alpha=0.4)
Select a region of interest for detailed visualization
cs.plot_fov(hd, 20, 20; marker_size = 0, width=2000, height=2400)
df = cs.subset_fov(hd, [193,197,333,337], 20,20);
xlim, ylim = cs.get_subset_cood(df)
#((9518.974171506527, 17960.0047340106), (14384.796062010153, 20199.649061752563))
####这种初始裁剪可能不会精确地与所需的区域对齐。根据获得的坐标，您可以稍微调整限制以更好地适合您打算关注的区域:
xlim = (9418, 17960)
ylim=(15500, 19799)
####This method allows you to fine-tune the coordinates manually for an optimal close-up view of your region of interest.
cs.sp_dim_plot(hd, "cluster"; width=1800, height=1000, y_lims=ylim , x_lims=xlim,
do_legend=true, img_res = "high", stroke_width=0.4,
anno_color = Dict("11" => "slateblue1", "2"=>"green2","13"=>"yellow1"),
cell_highlight =["11", "2","13"],legend_size=40,
adjust_contrast =1, adjust_brightness=0.0, alpha=1)
Compare the spatial resolution of different bin sizes
xlim = (11818, 16100)
ylim=(15500, 17599)
hd = cs.set_default_layer(hd; layer_slot = "8_um")
cs.sp_dim_plot(hd, "cluster"; width=1500, height=900, y_lims=ylim , x_lims=xlim,
do_legend=true, img_res = "high", stroke_width=0.4,
anno_color = Dict("11" => "slateblue1", "2"=>"green2"),
cell_highlight =["11", "2"],legend_size=40,
adjust_contrast =1, adjust_brightness=0.0, alpha=0.7)
hd = cs.set_default_layer(hd; layer_slot = "16_um")
cs.sp_dim_plot(hd, "cluster"; width=1500, height=900, y_lims=ylim , x_lims=xlim,
do_legend=true, img_res = "high", stroke_width=0.4,
anno_color = Dict("13" => "slateblue1", "16"=>"green2"),
cell_highlight =["13", "16"],legend_size=40,
adjust_contrast =1, adjust_brightness=0.0, alpha=0.7)
使用sp_feature_plot在不同分辨率下可视化基因表达。
xlim = (11818, 16100)
ylim=(15500, 17599)
hd = cs.set_default_layer(hd; layer_slot="2_um")
hd = cs.normalize_object(hd; scale_factor = 10000)
hd = cs.convert_image_data(hd; layer_slot = "2_um")
cs.sp_feature_plot(hd, ["Pcp4"]; color_keys=["gray94", "cyan", "blue", "darkblue"],
width=800, height=500, x_lims=xlim , y_lims=ylim, img_res = "high", alpha=1)
hd = cs.set_default_layer(hd; layer_slot="8_um")
hd = cs.normalize_object(hd; scale_factor = 10000)
cs.sp_feature_plot(hd, ["Pcp4"]; color_keys=["gray94", "cyan", "blue", "darkblue"],
width=800, height=500, x_lims=xlim , y_lims=ylim, img_res = "high", alpha=1)
hd = cs.set_default_layer(hd; layer_slot="16_um")
hd = cs.normalize_object(hd; scale_factor = 10000)
cs.sp_feature_plot(hd, ["Pcp4"]; color_keys=["gray94", "cyan", "blue", "darkblue"],
width=800, height=500, x_lims=xlim , y_lims=ylim, img_res = "high", alpha=1)
写出数据
cs.save(hd; filename="visiumHD.jld2")
最后来看Xenium
import CellScopes as cs
xenium_dir = "/mnt/sdc/new_analysis_cellscopes/xenium_mouse_brain/"
@time brain = cs.read_xenium(xenium_dir; min_gene = 0, min_cell = 0, prefix = "brain")
Normalization, dimension reduction and cell clustering
brain = cs.normalize_object(brain)
brain = cs.scale_object(brain)
brain = cs.find_variable_genes(brain; nFeatures=200)
brain = cs.run_pca(brain; method=:svd, pratio = 1, maxoutdim = 20)
(可选)使用Baysor生成细胞多边形，并将基因表达分配给多边形(细胞分割)
Xenium对象已经包含由10x Analyzer生成的细胞多边形数据。如果您想使用其他工具(如Baysor)绘制细胞多边形，请按照以下步骤操作。在生成细胞多边形后，通过运行polygons_cell_mapping将每个多边形与最近的细胞基于欧几里得距离关联起来。最后，需要使用generate_polygon_counts将基因表达值分配给细胞多边形。这一步完全是可选的，可以省略而不影响任何后续的分析。
import Baysor as B
scale=20
min_pixels_per_cell = 15
grid_step = scale / min_pixels_per_cell
bandwidth= scale / 10
count_molecules = deepcopy(brain.spmetaData.molecule)
count_molecules.cell[count_molecules.cell.==-1] .= 0
count_molecules.cell = string.(count_molecules.cell)
cell_baysor = [replace(x, "brain_" => "") for x in count_molecules.cell]
cell_baysor = parse.(Int64, cell_baysor)
polygons = B.boundary_polygons(count_molecules, cell_baysor, grid_step=grid_step, bandwidth=bandwidth)
brain.polygonData = polygons
brain = cs.polygons_cell_mapping(brain)
brain = cs.generate_polygon_counts(brain)
可视化
可视化细胞和基因在降维空间和空间空间上的表达。
可视化细胞注释
colors =["#4f8c9d", "#8dd2d8", "#0a4f4e", "yellow", "#229743",
"#8dd2d8", "#0a4f4e","#a7d64e", "#788c3b", "#57e24c",
"#683c00", "#f1bb99", "blue", "#9f1845",
"#f87197", "#ff0087", "#a57a6a", "#fabd2a",
"#374475", "#628df2", "#691b9e", "slateblue1",
"#d5d0fa", "black", "#ab7b05","#4f8c9d",
"#8dd2d8", "#0a4f4e", "#4aeeb6", "#229743",
"#a7d64e", "#788c3b", "#57e24c", "#683c00",
"#f1bb99", "#db3c18", "cyan", "#f87197",
"#ff0087", "#a57a6a", "#fabd2a", "#374475",
"#628df2", "#691b9e", "#b25aed", "red",
"fuchsia", "#ab7b05","orangered3", "#0a4f4e"]
celltypes = string.(unique(brain.metaData.cluster))
anno_color=Dict(celltypes .=> colors)
cs.sp_dim_plot(brain, "cluster"; anno_color = anno_color,
do_label = false,marker_size = 5,
width=2500, height=2000, do_legend=false
)
Select a field of view to visualize the cell annotation
cs.plot_fov(brain, 20, 20)
上面的网格图允许我们选择感兴趣的特定区域进行可视化。
hippo = cs.subset_fov(brain, [93, 96, 173, 176],20,20)
x1 = minimum(hippo.x)
x2 = maximum(hippo.x)
y1 = minimum(hippo.y)
y2 = maximum(hippo.y)
cs.sp_dim_plot(brain, "cluster"; anno_color = anno_color,
do_label = false,marker_size = 7, x_lims = (x1, x2),
y_lims = (y1, y2),
width=700, height=600, do_legend=false
)
cs.plot_cell_polygons(brain, "cluster"; x_lims = (x1, x2),
y_lims = (y1, y2), cell_colors= colors, stroke_color="black",
width = 800, height = 550)
基因可视化
cs.sp_feature_plot(brain, "Bcl11b"; color_keys=["gray94", "lemonchiffon", "red"], height=3000, width=3000, marker_size = 4)
cs.plot_gene_polygons(brain, ["Bcl11b"]; x_lims = (x1, x2),y_lims = (y1, y2),
width = 800, height = 550, color_keys=["#440154", "#440154","#3b528b","#ffff67"], bg_color="black")
Gene imputation
提供包装函数运行基因插入使用三种不同的工具:SpaGE, gimVI和tangram。
data_path = "/mnt/sdc/new_analysis_cellscopes/brain_sc/"
SpaGE_path = "/mnt/sdc/new_analysis_cellscopes/SpaGE/"
brain = cs.run_spaGE(brain, data_path, SpaGE_path)
cs.sp_feature_plot(brain, ["Adcy5","Ccdc3", "Prox1"]; color_keys=["gray90", "lemonchiffon", "red"], use_imputed = true)
cs.sp_feature_plot(brain2, ["Ccdc3"];
color_keys=["gray80", "lemonchiffon", "red3"], order=true,
use_imputed = true, x_lims = (x1, x2), marker_size = 10,
y_lims = (y1, y2),width = 800, height = 550)
