单细胞测序可以检测一块组织中每个细胞的转录组情况，但其实一块组织中存在着不同发育状态的细胞，因此通过基因的表达情况我们可以了解一些细胞的发育状况和细胞转化的过程，目前为止做拟时序的软件就有monocle，velocity 等等，在这里我们对目前比较主流的进行测试，后期有发现了好用的软件也会持续进行更新。
测试1：monocle2:Monocle 2然后在自动选择的一组数据质心上构造一棵生成树（DDRTree算法）。然后，该算法将细胞移动到它们最近的树的顶点，更新顶点的位置以适应细胞，学习新的生成树，并迭代地继续这个过程，直到树和细胞的位置已经收敛。在这个过程中，Monocle 2保持了高维空间和低维空间之间的可逆映射，从而既学习了轨迹，又降低了数据的维数。官网学习地址Monocle (cole-trapnell-lab.github.io)
测试2：monocle3：官网推荐的版本，相比于DDRTree算法，UMAP算法进行轨迹推断可能会更直观，算法上也进行了优化，Monocle3正式继承了Monocle2中迭代的灵魂算法，并舍弃了其中被诟病的轨迹图生成部分，引入了以PAGA为蓝本的基于图谱的轨迹推断流程，相比于2，速度和内存消耗都感人了很多，官网学习地址Monocle (cole-trapnell-lab.github.io)
测试3 ：velocity.R:刚转录的mRNA含有内含子与外显子，经过剪切后可形成成熟的mRNA，这个方法就是通过细胞中含有的前体RNA和成熟RNA的含量来预测谱系发生的方向的。
测试4：scvelo：感觉比在R里边操作，顺畅了很多，RNA Velocity Basics — scVelo 0.2.5.dev6+g1805ab4.d20210829 documentation
测试5 scanpy ：PAGA是一种基于图形的分析方法，先通过 Louvain algorithm算法对细胞进行降维，生成低纬度的聚类图，基于聚类图进一步分析不同细胞类群之间的关系。官网教程如下Trajectory inference for hematopoiesis in mouse — Scanpy documentation (scanpy-tutorials.readthedocs.io)

测试1：monocle2
monocle2安装,BiocManager真的好用，能用BiocManager安装的建议用BiocManager安装

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("monocle")

1 加载包

suppressMessages({
library(reshape2)
library(Seurat)
library(corrplot)
library(grid)
library(cowplot)
library(tidyverse)
library(monocle)
library(ggsci)
library(ggpubr)
})

2 读入Seurat数据，本次用的是上次Seurat分析的数据，但是monocle的计算速度太令人绝望了，所以一直在各种subset，subset完还剩17000+ cells，还是很占用计算资源。

PRO<-readRDS('/test1/PRO1.rds')
Idents(PRO)<-'orig.ident'
PRO<-subset(PRO,idents = 'Bong_marrow1')
Idents(PRO)<-'cluster'
PRO<-subset(PRO,idents = 'Neutrophils',invert = TRUE)

3 将Seurat对象转为monocle对象

PRO_matrix <- as(as.matrix(PRO@assays$RNA@counts), 'sparseMatrix')
pd <- new('AnnotatedDataFrame', data = PRO@meta.data) 
fd <- new('AnnotatedDataFrame', data = data.frame(gene_short_name = row.names(PRO),row.names = row.names(PRO)))
cds <- newCellDataSet(PRO_matrix,
                      phenoData = pd,
                      featureData = fd,
                      lowerDetectionLimit = 0.5,
                      expressionFamily = negbinomial.size())

4 数据均一化，这一步挺慢的

cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)

5 过滤基因

cds <- detectGenes(cds, min_expr = 0.1)#在fDate（cds）中添加一列num_cells_expressed
expressed_genes <- row.names(subset(fData(cds),num_cells_expressed >= 10))#过滤掉小于10个细胞中表达的基因
#如果Seurat过滤过，这一步可以不用过滤

6 选择高变基因，这一步比较关键，算法也比较多，不同的算法可能会导致不同的结果，大家可以了解其原理，选择适合自己的算法

##使用seurat选择的高变基因
#express_genes <- VariableFeatures(PRO)
#cds <- setOrderingFilter(cds, express_genes)
#plot_ordering_genes(cds)
##使用clusters差异表达基因
#deg.cluster <- FindAllMarkers(PRO)
#express_genes <- subset(deg.cluster,p_val_adj<0.05)$gene
#cds <- setOrderingFilter(cds, express_genes)
#plot_ordering_genes(cds)
##使用monocle选择的高变基因
disp_table <- dispersionTable(cds)
disp.genes <- subset(disp_table, mean_expression >= 0.1 )
cds <- setOrderingFilter(cds, disp.genes$gene_id)
plot_ordering_genes(cds)
#dpFeature，本来想用这个方法做测试，但是跑了很久没跑出来，就放弃了，使用了monocle选择的高变基因来做
#diff <- differentialGeneTest(cds[expressed_genes,],fullModelFormulaStr="~cluster",cores=1)
#deg <- subset(diff, qval < 0.01) #选出2829个基因
#deg <- deg[order(deg$qval,decreasing=F),]
#write.table(deg,file="train.monocle.DEG.xls",col.names=T,row.names=F,sep="\t",quote=F)
#ordergene <- rownames(deg) 
#cds <- setOrderingFilter(cds, ordergene)
#pdf("train.ordergenes.pdf")
#plot_ordering

image.png

7 降维分析，细胞排序

cds <- reduceDimension(cds, max_components = 2, method = 'DDRTree')
cds <- orderCells(cds)

8差异基因

ordering_genes=names(tail(sort(apply(cds@assayData$exprs,1,mad)),2000)) #用作后边的差异基因分析，也可以用前边得到的disp.genes，seurat得到的差异基因，自己指定基因啥的
diff_test_res <- differentialGeneTest(cds[ordering_genes,],fullModelFormulaStr = "~cluster") #这里可以选择用各种基因集，比如说seurat中发现的高变基因，monocle中的基因集等，这里用到了上一步中基因，很多教程中用到了expreesed gene ，但是这个基因集基因有点多，这一步跑的巨慢
deg <- subset(diff_test_res, qval < 0.01)
deg <- deg[order(deg$qval,decreasing=F),]
##差异基因的结果文件保存
write.table(deg,file="train.monocle.DEG.xls",col.names=T,row.names=F,sep="\t",quote=F)
ordergene <- rownames(deg)
Time_diff <- differentialGeneTest(cds[ordergene,], cores = 1, fullModelFormulaStr = "~sm.ns(Pseudotime)")
#Time_diff <- Time_diff[,c(5,2,3,4,1,6,7)] #调整列的顺序，可以不用跑
write.csv(Time_diff, "Time_diff_all.csv", row.names = F)
Time_genes <- Time_diff %>% pull(gene_short_name) %>% as.character()
p=plot_pseudotime_heatmap(cds[Time_genes,], num_clusters=4, show_rownames=T, return_heatmap=T)
ggsave("Time_heatmapAll.pdf", p, width = 5, height = 10)

9 绘制拟时序分析图

p1<-plot_cell_trajectory(cds, color_by = "cluster")
p2<-plot_cell_trajectory(cds,color_by="Pseudotime", size=1,show_backbone=TRUE) 
p <- p1|p2
ggsave("trajectory_plot.pdf", plot = plot3, width = 16, height = 8)
p3<-plot_cell_trajectory(cds, color_by = "cluster") + facet_wrap("~cluster", nrow = 1)#可以看每个cluster的分布情况
p3
p1<-plot_complex_cell_trajectory(cds, x = 1, y = 2,color_by = "cluster")#+scale_color_manual() +theme(legend.title = element_blank())
p2<-ggplot(pData(cds), aes(Pseudotime, colour = cluster, fill=cluster)) +geom_density(bw=0.5,size=1,alpha = 0.5)+theme_classic2()
p4<-p1|p2#绘制树状图和细胞密度图
ggsave("trajectory_2plot.pdf", plot = p4, width = 16, height = 8)

p

p3

p4

还可以选择不同
10分支点分析
选择自己感兴趣的细胞亚型和分化路径，然后可以看发生变化的基因和基因的变化趋势

BEAM_res=BEAM(cds[Time_genes,],branch_point = 3)#在这里我们看一下第三个节点的分化基因
BEAM_res <- BEAM_res[order(BEAM_res$qval),] #将得到的基因进行排序
write.csv(BEAM_res, "BEAM_res.csv", row.names = F)# 这一笔跑了好久，为了之后不重复的跑，把它存在表格里，以后直接读取就好。
#plot_genes_branched_heatmap(cds[row.names(subset(BEAM_res,qval < 1e-4)),],branch_point = 3 ，num_clusters = 4, use_gene_short_name = T,show_rownames = T)#有737个基因个gene，太多了，展示在图上不好看，所以下边选了top50的基因进行展示，这一步跑了比较久
BEAM_genes <- top_n(BEAM_res, n = 50, desc(qval)) %>% pull(gene_short_name) %>% as.character()
#在这里可以看不同的基因在不同的分化的时候的表达情况
p<-plot_genes_branched_heatmap(cds[BEAM_genes,],
                            branch_point = 3, #绘制的是哪个分支
                            num_clusters = 3, #分成几个cluster，根据需要调整
                            show_rownames = T,#热图上现实基因名
                            return_heatmap = T)

image.png

这时候还可以看一下基因的变化

keygenes <- head(BEAM_genes ,4)#在这里我们选择top4的gene看一下
cds_subset <- cds[keygenes,]
p1 <- plot_genes_in_pseudotime(cds_subset, color_by = "State")
p2 <- plot_genes_in_pseudotime(cds_subset, color_by = "cluster")
p3 <- plot_genes_in_pseudotime(cds_subset, color_by = "Pseudotime")
P<-p1|p2|p3
P
ggsave("gene_plot.pdf", plot = P, width = 16, height = 8) #width和height可以指定图片的长宽比，是很好用的参数。

image.png

关于monocle2的介绍就先到这里把，这个真的是太耗cpu了

Monocle3
软件安装,有些包如果安装了，就可不必重新安装，如果过程中遇到报错，缺啥补啥，总的来说，没有yum权限，装的痛苦无比。

BiocManager::install(c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
                       'limma', 'S4Vectors', 'SingleCellExperiment',
                       'SummarizedExperiment', 'batchelor', 'Matrix.utils'))
install.packages("devtools")
devtools::install_github('cole-trapnell-lab/leidenbase')
devtools::install_github('cole-trapnell-lab/monocle3')

1加载包

suppressMessages({
library(reshape2)
library(Seurat)
library(corrplot)
library(grid)
library(cowplot)
library(tidyverse)
library(monocle3)
library(ggsci)
library(ggpubr)
})

2 读入Seurat的数据，为了增加分析速度，还是单挑出来一个样本进行分析，monocle3的分析速度真的是比2感人很多.....

PRO<-readRDS('/test1/PRO1.rds')
Idents(PRO)<-'orig.ident'
PRO<-subset(PRO,idents = 'Bong_marrow1')

3 生成monocle格式的数据,并进行均一化，降维分析

data <- GetAssayData(PRO, assay = 'RNA', slot = 'counts')
cell_metadata <- PRO@meta.data
gene_annotation <- data.frame(gene_short_name = rownames(data))
rownames(gene_annotation) <- rownames(data)
cds <- new_cell_data_set(data,cell_metadata = PRO@meta.data,gene_metadata = gene_annotation)
cds <- preprocess_cds(cds, num_dim = 50) #preprocess_cds函数相当于seurat中NormalizeData+ScaleData+RunPCA
#plot_pc_variance_explained(cds)
cds <- reduce_dimension(cds, preprocess_method = "PCA")
cds <- cluster_cells(cds) 
##从seurat导入整合过的umap坐标
cds.embed <- cds@int_colData$reducedDims$UMAP
int.embed <- Embeddings(PRO, reduction = "umap")
int.embed <- int.embed[rownames(cds.embed),]
cds@int_colData$reducedDims$UMAP <- int.embed
#识别轨迹
cds <- learn_graph(cds)
plot_cells(cds, color_cells_by = "celltype")

image.png

4 细胞排序

#手动选择root
#cds <- order_cells(cds) 教程中使用的代码，但是由于juputer中不能出现交互界面，所以分析的时候不能用，没有使用jupyter的小伙伴可以试试是否可以出现交互界面，可以直接在交互界面上进行选择，否则，可能得想点其它的办法。
p<-plot_cells(cds, color_cells_by = "celltype")
p + geom_vline(xintercept = seq(8,9,0.25)) + geom_hline(yintercept = seq(-5,-4,0.25)) #使用坐标定位我们想要的root，定位完选择我们想要的格子，在这里我们选择交叉的右上角
embed <- data.frame(Embeddings(PRO, reduction = "umap"))
embed <- subset(embed, UMAP_1 > 8.75 & UMAP_1 < 9 & UMAP_2 > -4.25 & UMAP_2 < -4)
root.cell <- rownames(embed)
cds <- order_cells(cds, root_cells = root.cell)#可以选择多个root
p2<-plot_cells(cds, color_cells_by = "pseudotime", label_cell_groups = FALSE, label_leaves = FALSE,  label_branch_points = FALSE)
p<-p1|p2
p
ggsave("odercell.pdf", plot = p, width = 10, height = 5)

image.png

5 寻找轨迹分析差异基因

Track_genes <- graph_test(cds, neighbor_graph="principal_graph", cores=2) # 如果出现了Error: 'rBind' is defunct.的报错，运行trace('calculateLW', edit = T, where = asNamespace("monocle3"))然后将出现的脚本中的 Matrix::rBind 改为 rbind
#trace('calculateLW', edit = T, where = asNamespace("monocle3"))
write.csv(Track_genes, "Trackgenes1.csv", row.names = T)
Track_genes_sig <- Track_genes %>% top_n(n=10, morans_I) %>%pull(gene_short_name) %>% as.character()
p <- plot_genes_in_pseudotime(cds[Track_genes_sig,], color_cells_by="celltype", min_expr=0.5, ncol = 2)
p
ggsave("Genes_Jitterplot.pdf", plot = p, width = 8, height = 6)
# plot_cells(cds, genes=Track_genes_sig, show_trajectory_graph=FALSE,label_cell_groups=FALSE,  label_leaves=FALSE) #也可以查看这些基因在featureplot中的表达情况

image.png

6 共表达模块分析

genelist <- pull(Track_genes, gene_short_name) %>% as.character()
gene_module <- find_gene_modules(cds[genelist,], resolution=0.01, cores = 6)#resolution越大，在这里分出来的module越多
#write.csv(gene_module, "Genes_Module.csv", row.names = F)
cell_group <- tibble::tibble(cell=row.names(colData(cds)), cell_group=colData(cds)$celltype)
agg_mat <- aggregate_gene_expression(cds, gene_module, cell_group)
row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))
p <- pheatmap::pheatmap(agg_mat, scale="column", clustering_method="ward.D2")
ggsave("Genes_Module.pdf", plot = p, width = 8, height = 8)

共表达模块

3安装velocyto
直接在原来的singlecell_test环境下安装

pip install -i https://pypi.tuna.tsinghua.edu.cn/simple velocyto

1 首先我们先加载一下包

suppressMessages({
library(loomR)
library(velocyto.R)
library(Seurat)
library(SeuratWrappers)
library(SeuratDisk)
library(ggplot2)
})

2 接下来我们读入数据，还是使用同一套数据,为了降低运算量，我们就选择一个样本来做

PRO<-readRDS('./PRO1.rds')
Idents(PRO)<-'orig.ident'
PRO<-subset(PRO,idents = 'Bong_marrow1')
ldat <- ReadVelocity(file = "./looms.loom") #读入loom数据

3 整理数据格式

#统一seurat和velocity的名称
colnames(ldat$spliced)<-gsub("x","_1",colnames(ldat$spliced))
colnames(ldat$spliced)<-gsub("looms:","",colnames(ldat$spliced))
colnames(ldat$unspliced)<-colnames(ldat$spliced)
colnames(ldat$ambiguous)<-colnames(ldat$spliced)
## 由于Seurat是选用了筛选后的barcode，而velocyto并没有，在这里我们将细胞统一一下的
ldat$spliced<-ldat$spliced[,rownames(PRO@meta.data)]
ldat$unspliced<-ldat$unspliced[,rownames(PRO@meta.data)]
ldat$ambiguous<-ldat$ambiguous[,rownames(PRO@meta.data)]
emat <- ldat$spliced
#emat <- emat[,colSums(emat)>=1e3] 根据需要看数据需不需过滤
nmat <- ldat$unspliced # disregarding spanning reads, as there are too few of them

4 数据过滤

cluster.label <- PRO@active.ident
emat <- filter.genes.by.cluster.expression(emat,cluster.label,min.max.cluster.average = 0.2)#筛选基因
nmat <- filter.genes.by.cluster.expression(nmat,cluster.label,min.max.cluster.average = 0.05)
length(intersect(rownames(emat),rownames(nmat)))
PROumap <- PRO@reductions$umap@cell.embeddings 
cell.dist <- as.dist(1-armaCor(t(PRO@reductions$umap@cell.embeddings)))#计算细胞距离
fit.quantile <- 0.02
rvel.cd <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells=25,cell.dist=cell.dist,fit.quantile=fit.quantile) #计算RNA速率

5 绘图

Idents(PRO)<-'cluster'
gg<-UMAPPlot(PRO)
gg
ggplot_build(gg)$data
colors <- as.list(ggplot_build(gg)$data[[1]]$colour) 
p1 <- show.velocity.on.embedding.cor(PROumap,rvel.cd,n=1000,scale='sqrt',cell.colors=ac(colors,alpha=0.5),cex=0.8,arrow.scale=2,show.grid.flow=T,min.grid.cell.mass=1.0,grid.n=50,arrow.lwd=1,do.par=F,cell.border.alpha =0.1,USE_OPENMP=1,n.cores=24,main="Cell Velocity") #n可以好好设置一下，细胞数多的时候不宜设的太小
p1

p1

测试4：
使用R语言先输出需要的信息

suppressMessages({
library(loomR)
library(velocyto.R)
library(Seurat)
library(SeuratWrappers)
library(SeuratDisk)
library(ggplot2)
})

读入数据，输出我们需要的数据

PRO<-readRDS('./PRO1.rds')
Idents(PRO)<-'orig.ident'
PRO<-subset(PRO,idents = 'Bong_marrow1')
#提取坐标信息
write.csv(Embeddings(PRO, reduction = "umap"), file = "cell_embeddings.csv")
# 获取每个细胞的barcode
write.csv(Cells(PRO), file = "cellID_obs.csv", row.names = FALSE)
# 提取每个细胞的cluster信息
write.csv(PRO@meta.data[, 'RNA_snn_res.0.8', drop = FALSE], file = "cell_clusters.csv")
# 提取每个细胞的celltype信息
write.csv(PRO@meta.data[, 'cluster', drop = FALSE], file = "cell_celltype.csv")
library(scales)
# 获取celltype的颜色信息
hue_pal()(length(levels(PRO$cluster)))
# 获取cluster的颜色信息
hue_pal()(length(levels(PRO$cluster)))

下面使用python进行分析
1 导入包

import scanpy
import scvelo as scv
import pandas as pd
import numpy as np
import os

2 读入数据

loom_data = scv.read('./looms.loom', cache=False)
# barcode名字去重后缀x，与seurat导出的barcode名称一致
loom_data.obs = loom_data.obs.rename(index = lambda x: x.replace('looms:','').replace('x', '_1'))
# 读取seurat中的meta信息
sample_obs = pd.read_csv('./cellID_obs.csv')
cell_umap= pd.read_csv('./cell_embeddings.csv', header=0, names=["Cell ID", "UMAP_1", "UMAP_2"])
cell_clusters = pd.read_csv('./cell_clusters.csv', header=0, names=["Cell ID", "cluster"])
cell_celltype = pd.read_csv('./cell_celltype.csv', header=0, names=["Cell ID", "celltype"])
# 对细胞文件和RNA剪切速率文件取交集，保留关注的细胞类型
sample_one = loom_data[np.isin(loom_data.obs.index, sample_obs)]
# 构建umap坐标, cluster, celltype信息数据框
sample_one_index = pd.DataFrame(sample_one.obs.index)
sample_one_index = sample_one_index.rename(columns = {'CellID':'Cell ID'})
umap_ordered = sample_one_index.merge(cell_umap, on = 'Cell ID')
celltype_ordered = sample_one_index.merge(cell_celltype, on = "Cell ID")
celltype_ordered.head()
clusters_ordered = sample_one_index.merge(cell_clusters, on = "Cell ID")
clusters_ordered.head()
umap_ordered = umap_ordered.iloc[:,1:]
clusters_ordered = clusters_ordered.iloc[:,1:]
celltype_ordered = celltype_ordered.iloc[:,1:]
sample_one.obsm['X_umap'] = umap_ordered.values
sample_one.uns['clusters'] = clusters_ordered.values
sample_one.obs['celltype'] = celltype_ordered.values
adata = sample_one
# some gene labels are duplicated (Ensembl IDs are still unique!!)
adata.var_names_make_unique()
# save model to file
adata.write('Allcelltype_dynamicModel.h5ad', compression = 'gzip')
# read h5ad file
#adata= scv.read('Allcelltype_dynamicModel.h5ad')

3 进行数据整理

scv.pp.filter_and_normalize(adata,min_shared_counts=30, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)#在基本预处理（基因选择和归一化）之后，我们计算一阶和二阶矩（均值和非中心方差）以进行速度估计：
scv.tl.velocity(adata, mode = "stochastic")#核心步骤，估计转录速率
scv.tl.velocity_graph(adata) #[Getting Started — scVelo 0.2.5.dev6+g1805ab4.d20210829 documentation](https://scvelo.readthedocs.io/getting_started/)没读懂这一步的算法，感兴趣的小伙伴自己去看吧 
ident_colours = ['#F8766D','#E38900','#C49A00','#99A800','#53B400','#00BC56','#00C094','#00BFC4','#00B6EB','#06A4FF','#A58AFF','#DF70F8','#FB61D7','#FF66A8']#Seurat导出来的颜色，不用也没关系，会有自动生成的颜色
scv.pl.velocity_embedding_stream(adata, basis='X_umap',color = "celltype", palette = ident_colours, size = 20,alpha =0.8) #可以进行绘图了
scv.pl.velocity_embedding_grid(adata, basis='umap', color='celltype', save='embedding_grid.pdf', title='', scale=0.25)#另一种展示形式

image.png

image.png

scv.pl.proportions(adata, groupby='celltype') #未剪切的mRNA多可能说明这个cluster转录比较活跃

image.png

# 看一下我们关注的基因的表达情况
scv.pl.velocity(adata, var_names=['Camp'], color='celltype')

image.png

测试5
1 加载需要的包

import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib.pyplot as plt

2 读入数据并进行预处理

data = sc.read_h5ad("bonemarrow.h5ad")
results_file = './bomarrow.h5ad' #写入的文件
data.X = data.X.astype('float64')  
sc.tl.draw_graph(data)
sc.pl.draw_graph(data, color='clusters', legend_loc='on data',title = "")

image.png

3 去噪

sc.tl.diffmap(data)
sc.pp.neighbors(data, n_neighbors=10, use_rep='X_diffmap')
sc.tl.draw_graph(data)
sc.pl.draw_graph(data, color='clusters', legend_loc='on data')#legend_loc='right margin'

image.png

4 聚类于PAGA分析
看不同cluster的时候是可以指定颜色的，具体教程在Trajectory inference for hematopoiesis in mouse — Scanpy documentation (scanpy-tutorials.readthedocs.io)

sc.tl.paga(data,groups='clusters')
sc.tl.paga(data, groups='clusters')
sc.pl.paga(data, threshold=0.03, show=False)#查看不同cluster中的相互关系
sc.pl.paga(data,color=['Lyz2','Vcan'])#点的大小代表该群细胞数的大小，粗细代表2个群之间的联系，查看marker基因在不同cluster中的表达情况

image.png

image.png

#利用PAGA重新计算细胞之间的距离
sc.tl.draw_graph(data, init_pos='paga')
sc.pl.draw_graph(data, color=['clusters','Lyz2','Vcan','Fn1','Ccr2','Ly6c2'], legend_loc='on data')

image.png

sc.pl.paga_compare(data, threshold=0.03, title='', right_margin=0.2,size=10, edge_width_scale=0.5,legend_fontsize=12, fontsize=12, frameon=False, edges=True)

image.png

就画到这里吧，总的来说，这方面的教程比较少，理解起来和每步操作起来的逻辑还没有很顺，等有时间慢慢的再学习补充吧