作者碗殷，Evil Genius

哎呀， Xenium和HD的空間項目進(jìn)來了速缨，又要有很多新的工作要做了锌妻。

很多分析呢，大家如果只是自己玩一玩旬牲，那隨便做做就也可以了仿粹，但如果像我一樣作為公司的分析人員搁吓，尤其是大公司的核心分析人員，生信經(jīng)理之類的吭历，那一定要對自己的分析結(jié)果負(fù)責(zé)堕仔，對要求就會高得多,因為分析的水平代表了公司水平。

還有很多人我覺得態(tài)度上就有問題晌区，認(rèn)為官方代碼都是現(xiàn)成的摩骨，完全不需要公司分析，這種偏見從我做生信就一直有朗若，其實就跟我們上學(xué)一樣恼五，上高中大家用的一樣的教材，都是看的一樣的內(nèi)容哭懈，但是高考的結(jié)果卻是千差萬別灾馒。cell、nature等高分雜志的文章代碼很多都是公開的银伟，但是大家用這些代碼能發(fā)cell、nature么绘搞？再比如單細(xì)胞平臺彤避，10X單細(xì)胞的平臺原理大家都知道，但是國內(nèi)做的水平就是比國外要差夯辖，為什么琉预？

其實我們很多都缺乏一種精神，這種精神是在西方世界蒿褂、包括日本很推崇并且深入骨髓的圆米，那就是工匠精神。

關(guān)于Xenium啄栓，主推squidpy娄帖，文章在高精度空間轉(zhuǎn)錄組分析之squidpy和空間網(wǎng)絡(luò)圖,主要分析就是降維聚類和鄰域富集分析。

Xenium中昙楚，Seurat也出了示例教程近速，在Analysis of Image-based Spatial Data in Seurat

關(guān)于HD ，之前分享過STEP堪旧，文章在10Xvisium HD高精度平臺探索, 主要分析精度為8um 和 16um

高精度HD目前Seurat也更新的教程削葱，在Analysis, visualization, and integration of Visium HD spatial datasets with Seurat ? Seurat (satijalab.org)

總體來講，python版本的分析更加優(yōu)秀淳梦，這一篇對visium析砸、HD、Xenium

基礎(chǔ)分析進(jìn)行全系列的腳本更新(因為這是公司的要求)爆袍。

我們首先來看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首繁，示例數(shù)據(jù)在 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.蛮瞄，多個精度的分析結(jié)果我們都要保留(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選擇，我們需要選擇需要分析的精度闲先，當(dāng)然状土，也可以是多個精度

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已經(jīng)為8微米和16微米的桶尺寸提供了cluster和UMAP結(jié)果伺糠。

從8μm bin大小重新聚集細(xì)胞可能會很耗時蒙谓，因為spot數(shù)量通常超過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)

要使方向與原始組織方向?qū)R训桶，可以使用convert_image_data函數(shù)累驮。下面是調(diào)整方向的方法:

@time hd = cs.convert_image_data(hd)

轉(zhuǎn)換圖像數(shù)據(jù)后，所有圖形將使用調(diào)整后的坐標(biāo)進(jìn)行可視化舵揭。

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))
####這種初始裁剪可能不會精確地與所需的區(qū)域?qū)R谤专。根據(jù)獲得的坐標(biāo)，您可以稍微調(diào)整限制以更好地適合您打算關(guān)注的區(qū)域:
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在不同分辨率下可視化基因表達(dá)午绳。

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)

寫出數(shù)據(jù)

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生成細(xì)胞多邊形置侍，并將基因表達(dá)分配給多邊形(細(xì)胞分割)

Xenium對象已經(jīng)包含由10x Analyzer生成的細(xì)胞多邊形數(shù)據(jù)。如果您想使用其他工具(如Baysor)繪制細(xì)胞多邊形拦焚，請按照以下步驟操作蜡坊。在生成細(xì)胞多邊形后，通過運行polygons_cell_mapping將每個多邊形與最近的細(xì)胞基于歐幾里得距離關(guān)聯(lián)起來赎败。最后秕衙，需要使用generate_polygon_counts將基因表達(dá)值分配給細(xì)胞多邊形。這一步完全是可選的僵刮，可以省略而不影響任何后續(xù)的分析灾梦。

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)

可視化

可視化細(xì)胞和基因在降維空間和空間空間上的表達(dá)。

可視化細(xì)胞注釋

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)

上面的網(wǎng)格圖允許我們選擇感興趣的特定區(qū)域進(jìn)行可視化妓笙。

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

提供包裝函數(shù)運行基因插入使用三種不同的工具: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)

生活很好，有你更好