ixxmu
commented
1 year ago Seurat空间转录组分析（二）标准分析流程 by 生信随笔
上期推文 Seurat空间转录组分析（一）数据读入 介绍了三种方法将空转数据读入R语言中。本期推文是跟着Seurat官网对空转数据进行标准化分析（https://satijalab.org/seurat/articles/spatial_vignette.html）。值得注意的是，本文基于Seurat v4进行分析。在撰写本文时，Seurat官网进行了更新，目前最新版本是Seurat v5测试版。Seurat v5相比于v4，有更高的计算性能和更低的硬件要求。在Seurat v5中，作者引入了新方法可以更好的分析大样本的单细胞数据集，此外还极大地完善了单细胞多组学数据处理的需求，特别是对单细胞空间转录组数据的处理上进行了极大的提升。参见：
需要注意的是，目前v5仍存在诸多bug，有待开发者去解决，所以建议初学者谨慎尝试。总的来说，尽管v5将至，但是其分析的标准流程仍然和v4差不多，所以了解和熟悉Seurat v4的空间转录组标准流程，对后续v5版本的学习具有指导意义。
一. 读入数据
Step1. R包加载及安装
library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
Step2. 加载Seurat官网的示例数据
示例数据在https://support.10xgenomics.com/spatial-gene-expression/datasets获得，并使用Load10X_Spatial()函数将其加载到Seurat。这将读取spaceranger管道的输出，并返回Seurat对象，该对象包含spot级别的表达数据以及相关的组织切片图像。
本文的示例数据是10X平台的，内置于SeuratDataR包，加载方式如下：
InstallData("stxBrain")
brain <- LoadData("stxBrain", type = "anterior1")
brain
# An object of class Seurat
# 31053 features across 2696 samples within 1 assay
# Active assay: Spatial (31053 features, 0 variable features)
## 检查示例数据
head(brain@meta.data)
# orig.ident nCount_Spatial nFeature_Spatial slice region
# AAACAAGTATCTCCCA-1 anterior1 13069 4242 1 anterior
# AAACACCAATAACTGC-1 anterior1 37448 7860 1 anterior
# AAACAGAGCGACTCCT-1 anterior1 28475 6332 1 anterior
# AAACAGCTTTCAGAAG-1 anterior1 39718 7957 1 anterior
# AAACAGGGTCTATATT-1 anterior1 33392 7791 1 anterior
# AAACATGGTGAGAGGA-1 anterior1 20955 6291 1 anterior
SpatialDimPlot(brain,alpha = 0)
#更多参数参考
?SpatialDimPlot()
二. 数据预处理
在spot中基因表达数据进行初始的预处理步骤与典型的scRNA-seq相似。我们首先需要对数据进行规范化，以考虑数据点之间测序深度的差异。我们注意到，对于空间数据集来说，分子数/点的差异可能是巨大的，特别是在整个组织的细胞密度存在差异的情况下。我们在这里看到大量的异质性，这需要有效的标准化。
plot1 <- VlnPlot(brain, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(brain, features = "nCount_Spatial") + theme(legend.position = "right")
p = wrap_plots(plot1, plot2)
#ggsave(p,filename = "../Output/test.pdf")
p3 <- VlnPlot(brain, features = "nFeature_Spatial",pt.size = 0.1,cols = "tomato",group.by = "orig.ident") + NoLegend()
p4 <- FeatureScatter(brain, feature1 = "nCount_Spatial",feature2= "nFeature_Spatial",group.by = "orig.ident")+ NoLegend()
wrap_plots(p3, p4)
这些图表明，分子计数（molecular counts）在点间的差异不仅是技术性的，而且还依赖于组织解剖。例如，组织中神经元消耗的区域(如皮层白质)，可再生地显示出较低的分子计数。因此，标准方法(如LogNormalize函数)可能会有问题，因为它会强制每个数据点在标准化之后具有相同的底层“大小”。
作为一种替代方法，我们推荐使用sctransform (Hafemeister和Satija，已出版)，它构建了基因表达的正则化负二项模型，以便在保留生物差异的同时考虑技术因素。有关sctransform的更多信息，请参见 here(https://www.biorxiv.org/content/10.1101/576827v2)的预印和here(https://satijalab.org/seurat/sctransform_vignette.html)的Seurat教程。sctransform将数据归一化，检测高方差特征，并将数据存储在SCT分析中:
brain <- SCTransform(brain, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)
#提取表达量变化最高的10个基因
top10<- head(VariableFeatures(brain),10)
top10
##[1] "Ttr" "Hbb-bs" "Plp1" "Mbp" "Hba-a1" "Ptgds" "Penk" "Hba-a2" "S100a5" "Ppp1r1b"
#类似常规转录组，我们也可以展示高变基因；
p5 <- VariableFeaturePlot(brain,cols = c( "gray60", "red"))+NoLegend()
p6 <- LabelPoints(plot = p5,points = top10, repel = TRUE,xnudge=0,ynudge=0)
p5 | p6
Q1：空转数据需不需要进行过滤？
A：参见追风少年【关于10X空间转录组质控分析的讨论】https://www.jianshu.com/p/13fafd6567a4。据我的了解，基于10x Visium测序技术的空转数据，一般不需要对单个spot进行过滤，除非有一些质量极低的spot（例如极低的基因数量）。
另外，不同于单细胞数据分析，空转的数据强烈建议使用SCT转换进行标准化和归一化。
三. 基因表达可视化
在Seurat，我们有功能来探索和互动空间数据固有的视觉特性。Seurat中的SpatialFeaturePlot()函数扩展了FeaturePlot()，可以将分子数据覆盖在组织组织学上。例如，在这个小鼠大脑的数据集中，Hpca基因是强海马marker，Ttr基因是脉络膜丛marker。
SpatialFeaturePlot(brain, features = c("Hpca", "Ttr"))
Seurat的默认参数强调分子数据的可视化。但是，你也可以调整斑点的大小(和它们的透明度)来改善组织学图像的可视化，通过改变以下参数:
pt.size。因素-这将比例大小的斑点。默认为1.6 alpha -最小和最大透明度。默认是c(1,1) 尝试设置为alpha c(0.1, 1)，以降低表达式较低的点的透明度 stroke：设置spots的描边粗细,默认stroke = 0.25
p7 <- SpatialFeaturePlot(brain, features = "Ttr", pt.size.factor = 1)
p8 <- SpatialFeaturePlot(brain, features = "Ttr", alpha = c(0.1, 1))
p7 + p8
四. 降维聚类及可视化
然后，我们可以继续对RNA表达数据进行降维和聚类，使用与我们用于scRNA-seq分析相同的工作流程：
brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)
brain <- FindNeighbors(brain, reduction = "pca", dims = 1:30)
brain <- FindClusters(brain, verbose = FALSE)
brain <- RunUMAP(brain, reduction = "pca", dims = 1:30)
我们可以在UMAP空间，使用DimPlot或使用SpatialDimPlot将聚类的结果显示在图像上：
p9 <- DimPlot(brain, reduction = "umap", label = TRUE)
p10 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)
p9 + p10
因为有许多颜色，所以可视化哪个颜色属于哪个簇是很有挑战性的。我们有一些策略来帮助解决这个问题。通过设置label参数，可以在每个集群的中间位置放置一个彩色框（参见上面的图）以及do.hover。SpatialDimPlot的悬停参数允许交互式查看当前的spot标识。
SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(2, 1, 4, 3,
5, 8)), facet.highlight = TRUE, ncol = 3)
## 或者
cluster <- CellsByIdentities(object = brain, idents = c(1, 2, 3, 4,5,8))
str(cluster)
##List of 6
##$ 1: chr [1:291] "AAAGGCTCTCGCGCCG-1" "AAATCGTGTACCACAA-1" "AAATGATTCGATCAGC-1" "AAATTAACGGGTAGCT-1" ...
##$ 2: chr [1:288] "AAACAAGTATCTCCCA-1" "AAACTCGTGATATAAG-1" "AAATAGCTTAGACTTT-1" "AAATCTAGCCCTGCTA-1" ...
##$ 3: chr [1:215] "AAACAGAGCGACTCCT-1" "AAACCGGGTAGGTACC-1" "AAATAACCATACGGGA-1" "AAATTTGCGGGTGTGG-1" ...
##$ 4: chr [1:211] "AAACCGTTCGTCCAGG-1" "AAACGGTTGCGAACTG-1" "AAATCGCGGAAGGAGT-1" "AAATTACACGACTCTG-1" ...
##$ 5: chr [1:208] "AAACACCAATAACTGC-1" "AAACATGGTGAGAGGA-1" "AAACGAAGATGGAGTA-1" "AAACTTAATTGCACGC-1" ...
##$ 8: chr [1:155] "AAACATTTCCCGGATT-1" "AAACGGGCGTACGGGT-1" "AAGACATACGTGGTTT-1" "AAGAGGCCCTTTGGAA-1" ...
SpatialDimPlot(brain, cells.highlight = cluster, facet.highlight = TRUE,cols.highlight = c( "green", "red"), ncol = 3)
五. 交互式绘图可视化
SpatialDimPlot(brain, interactive = TRUE)
SpatialFeaturePlot(brain, features = "Ttr", interactive = TRUE)
LinkedDimPlot(brain)
六. 空间变量特征的识别
Seurat提供了两个工作流程来识别与组织空间位置相关的分子特征。
方法1：
第一种是根据组织内预先标注的解剖区域进行差异表达，这种差异表达可以通过非监督聚类或先验知识来确定。这种策略在这种情况下有效，因为上面的clusters显示出明显的空间差异：
de_markers <- FindMarkers(brain, ident.1 = 5, ident.2 = 6)
head(de_markers)
| p_val | avg_log2FC | pct.1 | pct.2 | p_val_adj | |
|---|---|---|---|---|---|
| Calb2 | 6.427214e-69 | 3.336874 | 1.000 | 0.537 | 1.135560e-64 |
| Camk2n1 | 1.519204e-68 | -2.450388 | 1.000 | 1.000 | 2.684130e-64 |
| Nrgn | 1.573095e-68 | -3.229826 | 0.971 | 1.000 | 2.779344e-64 |
| Stx1a | 2.104119e-68 | -2.238353 | 0.784 | 1.000 | 3.717558e-64 |
| Nptxr | 9.820482e-68 | -1.918293 | 0.942 | 1.000 | 1.735083e-63 |
| Lingo1 | 2.124777e-67 | -1.913281 | 0.880 | 1.000 | 3.754056e-63 |
可视化，绘制前3个差异基因的表达映射图：
SpatialFeaturePlot(object = brain, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3)
方法2：
在FindSpatiallyVariables中实现的另一种方法是在没有预注释的情况下搜索显示空间模式的特性。默认的方法(method = 'markvariogram')受到Trendsceek的启发，后者将空间转录组学数据建模为标记点过程，并计算一个variogram，它识别其表达水平取决于其空间位置的基因。更具体地说，这个过程计算伽玛（r）值，测量两个点之间一定的“r”距离的相关性。默认情况下，我们在这些分析中使用的r值为‘5’，并且只计算可变基因的这些值（其中的变异是独立于空间位置计算的），以节省时间。
## 3000个高变基因
length(VariableFeatures(brain))
## 从前1000个高变基因中查找差异基因
variable_top1000 = VariableFeatures(brain)[1:1000]
brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT",
features = variable_top1000,
selection.method = "markvariogram")
## 得到差异基因
DEGs<- SpatiallyVariableFeatures(brain, selection.method = "markvariogram")
length(DEGs) # 1000
## 可视化top差异化的高变基因
top.features <- head(DEGs, 6)
SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1)))
七. 可视化解剖区域的子集
与单细胞对象一样，您可以对该对象进行子集设置，以将重点放在数据的子集上。在这里，我们大致取到了额叶皮质的子集。这个过程也促进了这些数据与下一节的皮层scRNA-seq数据集的整合。首先，我们取clusters的一个子集，然后根据精确的坐标信息进一步“精修”。设置好亚组后，我们可以在完整图像或裁剪图像上看到皮质细胞。
Step1. 粗略的取到皮层得到子集
cortex <- subset(brain, idents = c(1, 2, 3, 4, 6, 7))
## 可视化1 2 3 4 6 7 clusters的子集结果
p11 <- SpatialDimPlot(cortex, crop = TRUE, label = TRUE, label.size = 3)
p11
Step2. 过滤，取皮层细胞
可以看到需要删除其他的细胞，这可以用SpatialDimPlot可视化，根据切片图像坐标（前缀需和Key保持一致），进一步去除多余的spots数据:
## 2.1 获取切片图像坐标和Key：
Key(object =cortex@images$anterior1)
img<- GetTissueCoordinates(cortex)
head(img)
## imagerow imagecol
##AAACAAGTATCTCCCA-1 385.9725 438.9501
##AAACAGAGCGACTCCT-1 163.3735 410.4991
##AAACCGGGTAGGTACC-1 336.5060 175.9208
##AAACCGTTCGTCCAGG-1 398.3649 225.6971
##AAACGGTTGCGAACTG-1 491.1016 286.1102
## 2.2 分群信息；
ident<- Idents(cortex)
head(ident)
## ident
##AAACAAGTATCTCCCA-1 2
##AAACAGAGCGACTCCT-1 3
##AAACCGGGTAGGTACC-1 3
##AAACCGTTCGTCCAGG-1 4
##AAACGGTTGCGAACTG-1 4
##AAACTCGTGATATAAG-1 2
## 2.3 合并切片图像坐标、分群信息；
imgId<- data.frame(img,ident)
head(imgId)
## imagerow imagecol ident
##AAACAAGTATCTCCCA-1 385.9725 438.9501 2
##AAACAGAGCGACTCCT-1 163.3735 410.4991 3
##AAACCGGGTAGGTACC-1 336.5060 175.9208 3
##AAACCGTTCGTCCAGG-1 398.3649 225.6971 4
##AAACGGTTGCGAACTG-1 491.1016 286.1102 4
##AAACTCGTGATATAAG-1 219.0361 478.0379 2
显微照片的像素坐标原点在左上角，而图表的原点在左下角，因此需要做垂直翻转y=540-imagerow：
#spots坐标信息的探索：
p12 <- ggplot(img, aes(x=imagecol, y=540-imagerow,color = ident)) +
geom_point(size = 1.6)+theme_bw()
p11 +p12
不垂直翻转的话：
p13 <- ggplot(img, aes(x=imagecol, y=imagerow,color = ident)) +
geom_point(size = 1.6)+theme_bw()
p11+p_13
根据上图切片图像坐标（前缀需和Key保持一致，即为anterior1）进一步去除多余的spots数据：
cortex <- subset(cortex, anterior1_imagerow > 400 | anterior1_imagecol < 150, invert = TRUE)
SpatialDimPlot(cortex, crop = TRUE, label = TRUE, label.size = 5)
cortex <- subset(cortex, anterior1_imagerow > 275 & anterior1_imagecol > 370, invert = TRUE)
SpatialDimPlot(cortex, crop = TRUE, label = TRUE, label.size = 5)
cortex <- subset(cortex, anterior1_imagerow > 250 & anterior1_imagecol > 440, invert = TRUE)
SpatialDimPlot(cortex, crop = TRUE, label = TRUE, label.size = 5)
最后，对切片特定区域的数据进行局部和整体可视化；
p14<- SpatialDimPlot(cortex, crop = TRUE, label = TRUE, label.size = 3)
p15<- SpatialDimPlot(cortex, crop = FALSE, label = TRUE,
pt.size.factor = 1, label.size = 3)
p14+p15
如需进一步分析的话，需要再次进行SCTransform转换，详见下一章节：
## After subsetting, we renormalize cortex
#cortex <- SCTransform(cortex, assay = "Spatial", verbose = FALSE) %>%
# RunPCA(verbose = FALSE)
## the annotation is stored in the 'subclass' column of object metadata
#DimPlot(allen_reference, group.by = "subclass", label = TRUE)
八. 整合空间转录组与单细胞数据
在~50um时，visium检测的斑点将包含多个细胞的表达谱。对于越来越多可用scRNA-seq数据的系统，用户可能有兴趣“反卷积”每个空间体素，以预测单元类型的底层组成。在准备这篇文章的过程中，我们测试了各种各样的脱卵方法和整合方法（decovonlution and integration methods），使用的是来自Allen研究所的参考scRNA-seq数据集（reference scRNA-seq dataset），其中包含了大约14000只成年小鼠的皮层细胞分类，并使用SMART-Seq2 protocol 生成。
作者认为，使用集成方法（ integration methods）（与反褶积方法deconvolution methods相反）可以获得更好的性能，这可能是因为空间和单细胞数据集的噪声模型本质上是不同的，而集成方法的特殊设计是为了对这些差异具有鲁棒性。因此我们应用“锚”的集成工作流一如最近在 Seurat v3介绍的，使注释的概率从一个reference 数据query 数据集转移。因此,我们遵循这里介绍的标签转移流程,利用sctransform正常化,但预测新方法被开发来完成这项任务。
我们首先加载数据（ https://www.dropbox.com/s/cuowvm4vrf65pvq/allen_cortex.rds?dl=1提供下载，990MB大小），预处理scRNA-seq 参考数据集，然后执行标签传输。该过程为每个spot输出每个scRNA-seq派生类的概率分类。我们将这些预测添加到Seurat对象中作为新的试验。
8.1 下载并读入单细胞数据
数据来源于《Adult mouse cortical cell taxonomy revealed by single cell transcriptomics》
allen_reference <- readRDS("../Rawdata/allen_cortex.rds")
allen_reference
##An object of class Seurat
##34617 features across 14249 samples within 1 assay
##Active assay: RNA (34617 features, 0 variable features)
8.2 数据预处理
首先对单细胞数据进行 SCTransform转换。注意，设置ncells=3000将使完整的数据集normalizes，但在3k时学习噪声模型，这大大提高了SCTransform的速度，并且没有性能损失：
# note that setting ncells=3000 normalizes the full dataset but learns noise models on 3k
# cells this speeds up SCTransform dramatically with no loss in performance
start_time <- Sys.time()
allen_reference <- SCTransform(allen_reference, ncells = 3000, verbose = FALSE) %>%
RunPCA(verbose = FALSE) %>%
RunUMAP(dims = 1:30)
end_time <- Sys.time()
end_time - start_time
allen_reference
##An object of class Seurat
##69225 features across 14249 samples within 2 assays
##Active assay: SCT (34608 features, 3000 variable features)
##1 other assay present: RNA
##2 dimensional reductions calculated: pca, umap
可视化，"subclass"是原文的注释结果：
# the annotation is stored in the 'subclass' column of object metadata
DimPlot(allen_reference, group.by = "subclass", label = TRUE)
对于取子集后的空间数据，再做一次SCTransform转换：
# After subsetting, we renormalize cortex
start_time <- Sys.time()
cortex <- SCTransform(cortex, assay = "Spatial", verbose = FALSE) %>%
RunPCA(verbose = FALSE)
end_time <- Sys.time()
end_time - start_time
##Time difference of 1.195518 mins
8.3 整合空间数据和单细胞数据
start_time <- Sys.time()
anchors <- FindTransferAnchors(reference = allen_reference, query = cortex, normalization.method = "SCT")
end_time <- Sys.time()
end_time - start_time
##Time difference of 1.0687 mins
predictions.assay <- TransferData(anchorset = anchors, refdata = allen_reference$subclass, prediction.assay = TRUE, weight.reduction = cortex[["pca"]], dims = 1:30)
cortex[["predictions"]] <- predictions.assay
现在，我们得到了每个class每个spot的预测分数。在额叶皮层区域特别有趣的是层状兴奋性神经元。在这里，我们可以区分这些神经元亚型的不同顺序层，例如:
#DefaultAssay(cortex) <- "SCT"
rownames(cortex) #基因名names
DefaultAssay(cortex) <- "predictions"
rownames(cortex) #预测的单细胞clusters
SpatialFeaturePlot(cortex, features = c("L2/3 IT", "L4"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)
根据这些预测分数，我们还可以预测位置受到空间限制的细胞类型。我们使用基于标记点过程的相同方法来定义空间变量特征，但使用细胞类型预测评分作为“标记”，而不是使用基因表达。
## 利用预测的细胞cluster，寻找高变特征（类似常规的高变基因 ）
cortex <- FindSpatiallyVariableFeatures(cortex, assay = "predictions", selection.method = "markvariogram",
features = rownames(cortex), r.metric = 5, slot = "data")
DECs = SpatiallyVariableFeatures(cortex)
top.clusters <- head(DECs, 4)
SpatialPlot(object = cortex, features = top.clusters, ncol = 2)
最后，作者证明Seurat的整合程序能够恢复已知的神经元和非神经元亚群的空间定位模式，包括层兴奋性亚群、第1层星形胶质细胞和皮层灰质：
SpatialFeaturePlot(cortex, features = c("Astro", "L2/3 IT", "L4", "L5 PT", "L5 IT", "L6 CT", "L6 IT",
"L6b", "Oligo"), pt.size.factor = 1, ncol = 3, crop = FALSE, alpha = c(0.1, 1))
九. 处理Seurat中的多个slices
这个小鼠大脑的数据集包含了另一个与大脑另一半相对应的切片。这里我们读取它并执行相同的初始规范化。
brain2 <- LoadData("stxBrain", type = "posterior1")
# An object of class Seurat
# 31053 features across 3353 samples within 1 assay
# Active assay: Spatial (31053 features, 0 variable features)
brain2 <- SCTransform(brain2, assay = "Spatial", verbose = FALSE)
brain2
# An object of class Seurat
# 48721 features across 3353 samples within 2 assays
# Active assay: SCT (17668 features, 3000 variable features)
# 1 other assay present: Spatial
## 这里读入的数据已经设置好了，如果没设置的话，手动设置一下
brain2@project.name
Idents(brain2)
brain2$orig.ident
##"Posterior1"
为了在同一个Seurat对象中处理多个切片，官方提供了merge函数：
start_time <- Sys.time()
brain.merge <- merge(brain, brain2)
end_time <- Sys.time()
end_time - start_time
##Time difference of 54.28318 secs
brain.merge
##An object of class Seurat
##49293 features across 6049 samples within 2 assays
##Active assay: SCT (18240 features, 0 variable features)
##1 other assay present: Spatial
## 可以保存一下再继续
##save(brain.merge,file = "brain.merge_sct.Rdata")
这样就可以对底层的RNA表达数据进行联合降维和聚类：
DefaultAssay(brain.merge) <- "SCT"
VariableFeatures(brain.merge) <- c(VariableFeatures(brain), VariableFeatures(brain2))
brain.merge <- RunPCA(brain.merge, verbose = FALSE)
brain.merge <- FindNeighbors(brain.merge, dims = 1:30)
brain.merge <- FindClusters(brain.merge, verbose = FALSE)
brain.merge <- RunUMAP(brain.merge, dims = 1:30)
brain.merge
##An object of class Seurat
##49293 features across 6049 samples within 2 assays
##Active assay: SCT (18240 features, 6000 variable features)
##1 other assay present: Spatial
##2 dimensional reductions calculated: pca, umap
最后，可以在单个UMAP图中共同可视化数据。SpatialDimPlot和SpatialFeaturePlot将默认将所有片绘制为列，将groupings/features 绘制为行。
DimPlot(brain.merge, reduction = "umap", group.by = c("ident", "orig.ident"))
接下来可视化比较两个切片的分群情况，默认使用两个slice的名称分组：
SpatialDimPlot(brain.merge)
#如果创建对象时没有对slice重命名，也可使用ident进行分组，绘图效果一样：
#SpatialDimPlot(brain.merge,group.by = c( "ident"))
#比较两个切片的特定基因表达水平
SpatialFeaturePlot(brain.merge, features = c("Hpca", "Plp1"))
十. 结语
关于整合空间转录组与单细胞数据整合方法汇总，详见2022年Nat Methods上的一篇Benchmarking文章《Benchmarking spatial and single-cell transcriptomics integration methods for transcript distribution prediction and cell type deconvolution》（DOI: 10.1038/s41592-022-01480-9），推荐大家使用cell2location或RCTD算法对空间转录组数据进行解卷积。值得注意的是，Seurat v5空转部分的教程使用了RCTD算法作为示例文档。
- END -
https://mp.weixin.qq.com/s/X8eP58iDEx4j7pSyGbze9Q