Seurat Dimplot函数学习总结

今天为了画这个cluster中怎么显示标签的图,研究了一个Seurat中怎么画这个图的,下面是学习过程中做的总结

运行例子

rm(list=ls())
library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
pbmc3k.final <- LoadData("pbmc3k", type = "pbmc3k.final")
#print(DimPlot(pbmc3k.final,label = TRUE, repel = TRUE))

# remove the legend
p=DimPlot(pbmc3k.final,label = TRUE, repel = TRUE)+theme(legend.position = "none")
print(p)

结果如下
在这里插入图片描述这个就是我想要达到的效果,把类别的标签放在类的中心显示,但是这个是Seurat包自己实现的,我要做的事就是把Seurat包的这个Dimplot函数给抽出来,变成自己的ggplot函数

调试

在这里插入图片描述其中核心部分是下面两个函数
在这里插入图片描述一个是SingleDimPlot函数, 另一个是LabelClusters函数,

首先发现Dimplot函数会进入函数SingleDimPlot函数,其函数实现如下
在这里插入图片描述虽然这个函数看起来很复杂,其实大部分都是在判断条件,所以对于自己的需求而言,很多代码可以全删了,经过SingleDimPlot函数后,画出来的效果如图
在这里插入图片描述可以看到,此时这个结果只是普普通通的ggplot画图,所以调整图例到细胞的类别中是在LabelClusters函数里, 其中LabelClusters函数实现如下
在这里插入图片描述
这个函数也是很多条件的,其实核心代码就那么几行, 大致看完这个代码之后,可以着手删代码来

自己实现

将Seurat对象转化到画图的csv文件

在这里插入图片描述

write.csv(data,file="./pbmc_UMAP.csv")

在这里插入图片描述

这里其实有一个问题,就是当把Seurat对象存成csv文件,然后再读入时,有些变量的格式是改变的,最常见的就是因子factor类型和chacracter类型,这个一定要注意,否则很容易出错, 比如在Seurat中,data的ident列是因子类型的
在这里插入图片描述当我把pbmc_UMAP.csv重新读取进来时,该数据框各列的类型如下
在这里插入图片描述可以看到,我重新读取进来的文件的ident列变成了character类型,这一点一定要多加注意,下面的代码也说明了这一点,需要额外进行类型转换

简化代码1

rm(list=ls())
setwd("/Users/yxk/Desktop/test/Seurat_label测试/")
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(cowplot))
suppressPackageStartupMessages(library(ggrepel))

data = read.csv("./pbmc_UMAP.csv",row.names = 1) # 相当于python里面的index
plot <- ggplot(data = data)

col.by ="ident"

plot = plot + geom_point(mapping = aes_string(x = "UMAP_1", 
      y = "UMAP_2", color = paste0("`", col.by, "`"), shape = NULL, 
      alpha = NULL), size = 0.6000758)

plot <- plot + guides(color = guide_legend(override.aes = list(size = 3))) + 
labs(color = NULL, title = col.by) 

plot <- plot + theme_cowplot()
print(plot)


# clusters这个是有问题的
LabelClusters = function(plot, id, clusters = NULL, labels = NULL, split.by = NULL, 
  repel = TRUE, box = FALSE, geom = "GeomPoint", position = "median", 
  ...) 
{
#   xynames <- unlist(x = GetXYAesthetics(plot = plot, geom = geom), 
#     use.names = TRUE)
  xynames=c("UMAP_1","UMAP_2")
  names(xynames)=c("x","y")
  if (!id %in% colnames(x = plot$data)) {
    stop("Cannot find variable ", id, " in plotting data")
  }
  if (!is.null(x = split.by) && !split.by %in% colnames(x = plot$data)) {
    warning("Cannot find splitting variable ", id, " in plotting data")
    split.by <- NULL
  }
  data <- plot$data[, c(xynames, id, split.by)]
  possible.clusters <- as.character(x = na.omit(object = unique(x = data[, 
    id])))
  groups <- possible.clusters
    
  if (any(!groups %in% possible.clusters)) {
    stop("The following clusters were not found: ", paste(groups[!groups %in% 
      possible.clusters], collapse = ","))
  }
  pb <- ggplot_build(plot = plot)
  if (geom == "GeomSpatial") {
    xrange.save <- layer_scales(plot = plot)$x$range$range
    yrange.save <- layer_scales(plot = plot)$y$range$range
    data[, xynames["y"]] = max(data[, xynames["y"]]) - data[, 
      xynames["y"]] + min(data[, xynames["y"]])
    if (!pb$plot$plot_env$crop) {
      y.transform <- c(0, nrow(x = pb$plot$plot_env$image)) - 
        pb$layout$panel_params[[1]]$y.range
      data[, xynames["y"]] <- data[, xynames["y"]] + sum(y.transform)
    }
  }
  data <- cbind(data, color = pb$data[[1]][[1]])
  labels.loc <- lapply(X = groups, FUN = function(group) {
    data.use <- data[data[, id] == group, , drop = FALSE]
    data.medians <- if (!is.null(x = split.by)) {
      do.call(what = "rbind", args = lapply(X = unique(x = data.use[, 
        split.by]), FUN = function(split) {
        medians <- apply(X = data.use[data.use[, split.by] == 
          split, xynames, drop = FALSE], MARGIN = 2, 
          FUN = median, na.rm = TRUE)
        medians <- as.data.frame(x = t(x = medians))
        medians[, split.by] <- split
        return(medians)
      }))
    }
    else {
      as.data.frame(x = t(x = apply(X = data.use[, xynames, 
        drop = FALSE], MARGIN = 2, FUN = median, na.rm = TRUE)))
    }
    data.medians[, id] <- group
    data.medians$color <- data.use$color[1]
    return(data.medians)
  })
  if (position == "nearest") {
    labels.loc <- lapply(X = labels.loc, FUN = function(x) {
      group.data <- data[as.character(x = data[, id]) == 
        as.character(x[3]), ]
      nearest.point <- nn2(data = group.data[, 1:2], query = as.matrix(x = x[c(1, 
        2)]), k = 1)$nn.idx
      x[1:2] <- group.data[nearest.point, 1:2]
      return(x)
    })
  }
  labels.loc <- do.call(what = "rbind", args = labels.loc)
#   labels.loc[, id] <- factor(x = labels.loc[, id], levels = levels(data[, 
#     id]))
  labels.loc$'ident'= as.factor(labels.loc$'ident') 
    
  labels <-  groups
  if (length(x = unique(x = labels.loc[, id])) != length(x = labels)) {
    stop("Length of labels (", length(x = labels), ") must be equal to the number of clusters being labeled (", 
      length(x = labels.loc), ").")
  }
  names(x = labels) <- groups
  for (group in groups) {
    labels.loc[labels.loc[, id] == group, id] <- labels[group]
  }
  #print(labels.loc)
  if (box) {
    geom.use <- ifelse(test = repel, yes = geom_label_repel, 
      no = geom_label)
    plot <- plot + geom.use(data = labels.loc, mapping = aes_string(x = xynames["x"], 
      y = xynames["y"], label = id, fill = id), show.legend = FALSE, 
      ...) + scale_fill_manual(values = labels.loc$color[order(labels.loc[, 
      id])])
  }
  else {
    geom.use <- ifelse(test = repel, yes = geom_text_repel, 
      no = geom_text)
    plot <- plot + geom.use(data = labels.loc, mapping = aes_string(x = xynames["x"], 
      y = xynames["y"], label = id), show.legend = FALSE, 
      ...)
  }
  if (geom == "GeomSpatial") {
    plot <- suppressMessages(expr = plot + coord_fixed(xlim = xrange.save, 
      ylim = yrange.save))
  }
  #print(plot)
  return(plot)
}

plot <- LabelClusters(plot = plot, id = "ident", repel = TRUE, 
  size = 4, split.by = NULL, box = FALSE, 
  color = "black")

print(plot)

结果如下
在这里插入图片描述
在这里插入图片描述

简化代码2

rm(list=ls())
setwd("/Users/yxk/Desktop/test/Seurat_label测试/")
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(cowplot))
suppressPackageStartupMessages(library(ggrepel))

data = read.csv("./pbmc_UMAP.csv",row.names = 1) # 相当于python里面的index

plot <- ggplot(data = data)

col.by ="ident"
plot = plot + geom_point(mapping = aes_string(x = "UMAP_1", 
                                              y = "UMAP_2", color = paste0("`", col.by, "`"), shape = NULL, 
                                              alpha = NULL), size = 0.6000758)

plot <- plot + guides(color = guide_legend(override.aes = list(size = 3))) + 
  labs(color = NULL, title = col.by) 

plot <- plot + theme_cowplot()
print(plot)


LabelClusters = function(plot, id, clusters = NULL, labels = NULL, split.by = NULL, 
                         repel = TRUE, box = FALSE, geom = "GeomPoint", position = "median", 
                         ...) 
{
  #   xynames <- unlist(x = GetXYAesthetics(plot = plot, geom = geom), 
  #     use.names = TRUE)
  xynames=c("UMAP_1","UMAP_2")
  names(xynames)=c("x","y")
  data <- plot$data[, c(xynames, id, split.by)]
  possible.clusters <- as.character(x = na.omit(object = unique(x = data[, id])))
  groups <- possible.clusters
  pb <- ggplot_build(plot = plot)
  data <- cbind(data, color = pb$data[[1]][[1]])
  labels.loc <- lapply(X = groups, FUN = function(group) {
    data.use <- data[data[, id] == group, , drop = FALSE]
    data.medians <- as.data.frame(x = t(x = apply(X = data.use[, xynames, drop = FALSE], MARGIN = 2, FUN = median, na.rm = TRUE)))
    
    data.medians[, id] <- group
    data.medians$color <- data.use$color[1]
    return(data.medians)
  })
  labels.loc <- do.call(what = "rbind", args = labels.loc)
  #labels.loc[, id] <- factor(x = labels.loc[, id], levels = levels(data[,  id]))
  labels.loc$'ident'= as.factor(labels.loc$'ident')
  
  labels <-  groups

  names(x = labels) <- groups
  for (group in groups) {
    labels.loc[labels.loc[, id] == group, id] <- labels[group]
  }
  #print(labels.loc)

  geom.use <- ifelse(test = repel, yes = geom_text_repel, 
                     no = geom_text)
  plot <- plot + geom.use(data = labels.loc, mapping = aes_string(x = xynames["x"], 
                                                                  y = xynames["y"], label = id), show.legend = FALSE,  ...)
  return(plot)
}

plot <- LabelClusters(plot = plot, id = "ident", repel = TRUE, 
                      size = 4, split.by = NULL, box = FALSE, 
                      color = "black")

print(plot)

结果如下
在这里插入图片描述

简化代码3

 rm(list=ls())
setwd("/Users/yxk/Desktop/test/Seurat_label测试/")
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(cowplot))
suppressPackageStartupMessages(library(ggrepel))

data = read.csv("./pbmc_UMAP.csv",row.names = 1) # 相当于python里面的index
plot <- ggplot(data = data)
col.by ="ident"
plot = plot + geom_point(mapping = aes_string(x = "UMAP_1", 
                                              y = "UMAP_2", color = paste0("`", col.by, "`"), shape = NULL, 
                                              alpha = NULL), size = 0.6000758)

plot <- plot + guides(color = guide_legend(override.aes = list(size = 3))) + 
  labs(color = NULL, title = col.by) 

plot <- plot + theme_cowplot()
print(plot)

###########

id = "ident"
size = 4
color = "black"
xynames=c("UMAP_1","UMAP_2")
names(xynames)=c("x","y")
data <- plot$data[, c(xynames, id)]
clusters<- as.character(x = na.omit(object = unique(x = data[, id])))
groups <- clusters
pb <- ggplot_build(plot = plot)
data <- cbind(data, color = pb$data[[1]][[1]])
labels.loc <- lapply(X = groups, FUN = function(group) {
  data.use <- data[data[, id] == group, , drop = FALSE]
  data.medians <- as.data.frame(x = t(x = apply(X = data.use[, xynames, drop = FALSE], MARGIN = 2, FUN = median, na.rm = TRUE)))
  data.medians[, id] <- group
  data.medians$color <- data.use$color[1]
  return(data.medians)
})
labels.loc <- do.call(what = "rbind", args = labels.loc)
#labels.loc[, id] <- factor(x = labels.loc[, id], levels = levels(data[,  id]))
labels.loc$'ident'= as.factor(labels.loc$'ident')

labels <-  groups

names(x = labels) <- groups
for (group in groups) {
  labels.loc[labels.loc[, id] == group, id] <- labels[group]
}
#print(labels.loc)

geom.use <- ifelse(test = TRUE, yes = geom_text_repel, 
                   no = geom_text)
plot <- plot + geom.use(data = labels.loc, mapping = aes_string(x = xynames["x"], 
                                                                y = xynames["y"], label = id), show.legend = FALSE)
print(plot)


结果如下
在这里插入图片描述

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/648941.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

学浪的缓存怎么导出来

学浪的缓存导出问题困扰着许多用户&#xff0c;备份和管理数据变得至关重要。在数字化时代&#xff0c;保护和利用数据是企业和个人不可或缺的需求。在这篇文章中&#xff0c;我们将深入探讨学浪缓存导出的方法&#xff0c;为您解决疑惑&#xff0c;让您轻松掌握数据的安全与便…

自定义类型详解(结构体,位段,枚举,联合体)

目录 结构体 结构体的自引用 匿名结构体 结构体内存对齐 结构体内存对齐的意义 修改默认对齐数 位段 位段的内存分配 位段的跨平台问题 位段的应用 枚举 枚举类型的定义 枚举和#define的区别 联合体&#xff08;共用体&#xff09; 联合体大小的计算 结构体 C语…

斐讯N1刷OpenWRT并安装内网穿透服务实现远程管理旁路由

文章目录 前言1. 制作刷机固件U盘1.1 制作刷机U盘需要准备以下软件&#xff1a;1.2 制作步骤 2. N1盒子降级与U盘启动2.1 N1盒子降级2.2 N1盒子U盘启动设置2.3 使用U盘刷入OpenWRT2.4 OpenWRT后台IP地址修改2.5 设置旁路由&无线上网 3. 安装cpolar内网穿透3.1 下载公钥3.2 …

08Django项目--用户管理系统--查(前后端)

对应视频链接点击直达 TOC 一些朋友加我Q反馈&#xff0c;希望有每个阶段的完整项目代码&#xff0c;那从今天开始&#xff0c;我会上传完整的项目代码。 用户管理&#xff0c;简而言之就是用户的增删改查。 08项目点击下载&#xff0c;可直接运行&#xff08;含数据库&…

3个完美恢复方案!苹果手机恢复视频就是这么简单

在日常生活中&#xff0c;我们使用苹果手机记录了许多珍贵的视频&#xff0c;包括家庭聚会、旅行经历等。但是我们不小心删除了&#xff0c;在“最近删除”文件夹中也找不到。 别担心&#xff0c;这并不困难&#xff01;无论您是意外删除了重要的视频还是遇到了其他数据丢失问…

PawSQL: 企业级SQL审核工具的新玩家

随着数据库应用在企业中的广泛使用&#xff0c;确保SQL代码质量的重要性日益凸显。现有的SQL审核工具很多&#xff0c;包括Yearning、goInception、Bytebase、爱可生的SQLE、云和恩墨的SQM等等&#xff0c;但是它们或者规则覆盖度、或者是在正确率等方面存在明显不足&#xff1…

【stm32】江科协听课笔记

[3-1] GPIO输出_哔哩哔哩_bilibili 5.GPIO输出 这里&#xff0c;寄存器就是一段特殊的存储器&#xff0c;内核可以通过APB2总线队寄存器进行读写&#xff0c;这样就可以完成输出/读取电平的功能。寄存器的每一位对应一个引脚&#xff0c;stm32是32位的&#xff0c;这里的寄存器…

K8S有了Service,为什么还要Ingress?

1、有了Service为什么还要Ingress? NodePort对外暴露端口存在的不足&#xff1a; 一个端口只能一个服务使用, 端口需要提前规划。 随着业务扩展, 端口的管理将是一个头疼的问题 只支持4层的负载均衡 LoadBalancer存在的不足&#xff1a; 贵、贵、贵。 要上云(俗话说上云…

前端Vue自定义顶部搜索框:实现热门搜索与历史搜索功能

前端Vue自定义顶部搜索框&#xff1a;实现热门搜索与历史搜索功能 摘要&#xff1a; 随着前端开发复杂性的增加&#xff0c;组件化开发成为了提高效率和降低维护成本的有效手段。本文介绍了一个基于Vue的前端自定义顶部搜索框组件&#xff0c;该组件不仅具备基本的搜索功能&am…

Android端 可使用Yolov5训练 路标识别

相信大家对于路标识别&#xff0c;红绿灯识别&#xff0c;图形识别opencv中也是一件烦人的事情&#xff0c;其参数是及其受到现实环境因素的影响的&#xff0c;那么今天我就给大家推荐一种方式&#xff0c;缺点是周期长&#xff0c;但其优点是如果训练效果好久对于环境的各种变…

安卓开发日志采集和分析面面谈

日志面面谈 为什么需要日志 复现问题&#xff0c;回溯到问题产生时候的系统状态&#xff0c;有利于定位和分析问题。 安卓日志有哪些&#xff1f; cpu 关注的纬度&#xff1a; 单个应用使用系统cpu分配温度 有什么用&#xff1a; App卡顿、ANRApp异常退出 怎么用&…

记一次 .NET某企业数字化平台 崩溃分析

一&#xff1a;背景 1. 讲故事 前些天群里有一个朋友说他们软件会偶发崩溃&#xff0c;想分析看看是怎么回事&#xff0c;所幸的是自己会抓dump文件&#xff0c;有了dump就比较好分析了&#xff0c;接下来我们开始吧。 二&#xff1a;WinDbg 分析 1. 程序为什么会崩溃 win…

从0开始回顾ElasticSearch

1 elasticsearch概述 1.1 elasticsearch简介 官网: https://www.elastic.co/ ElasticSearch是一个基于Lucene的搜索服务器。它提供了一个分布式多用户能力的全文搜索引擎&#xff0c;基于RESTful web接口。Elasticsearch是用Java开发的&#xff0c;并作为Apache许可条款下的…

芯课堂 | 芯片抗干扰测试方案

MCU芯片对所在环境中存在的电磁干扰须具有一定程度的抗扰度&#xff0c;确保使用该芯片的设备能正常运行。国际电工委员会&#xff08;IEC&#xff09;制定了多项国际标准&#xff0c;其中与MCU芯片相关的有IEC61000-4-2 &#xff08;静电&#xff09;&#xff0c; IEC61000-4-…

RK3568笔记二十六:音频应用

若该文为原创文章&#xff0c;转载请注明原文出处。 一、介绍 音频是我们最常用到的功能&#xff0c;音频也是 linux 和安卓的重点应用场合。 测试使用的是ATK-DLR3568板子&#xff0c;板载外挂RK809 CODEC芯片&#xff0c;RK官方驱动是写好的&#xff0c;不用在自己重新写。…

家居的3D交互展示用什么工具比较专业?

家居的3D交互展示可以使用多种专业工具来实现&#xff0c;这些工具不仅能够在手机和电脑上查看&#xff0c;还能在手机上进行交互操作&#xff0c;如放缩、旋转等&#xff0c;并且支持高清流畅的画面展示。以下是一些推荐的3D交互展示工具&#xff1a; 1、在线3D展示软件&…

牛客热题:寻找第K大

&#x1f4df;作者主页&#xff1a;慢热的陕西人 &#x1f334;专栏链接&#xff1a;力扣刷题日记 &#x1f4e3;欢迎各位大佬&#x1f44d;点赞&#x1f525;关注&#x1f693;收藏&#xff0c;&#x1f349;留言 文章目录 牛客热题&#xff1a;寻找第K大题目链接方法一&#…

Docker基础篇之Docker入门介绍

文章目录 1. 为什么要有Docker&#xff1f;2. Docker简介3. 容器和虚拟机的区别4. Docker下载 1. 为什么要有Docker&#xff1f; 假设我们现在正在开发一个项目&#xff0c;使用的是一台笔记本电脑而且开发环境具有特定的配置&#xff0c;其他开发人员身处的环境配置也各不相同…

ZeroTier+Nomachine远程

目录 前述&#xff1a;一、Zero二、Nomachine 前述&#xff1a; 需要远程控制时&#xff0c;服务端与客户端都必须下载这两个软件&#xff01;远程主机&#xff08;被控制的主机&#xff09;和远程客户端&#xff08;控制主机的用户&#xff09;都必须具有网络连接&#xff0c;…

地铁判官:啥时候B端系统界面,也出个“判官”,讲好不准打脸。

小编所在的城市——山东青岛&#xff0c;出了个地铁判官&#xff0c;我看了视频&#xff0c;哈哈哈&#xff0c;俗世的判断标准就是那么简单直接&#xff0c;而放到B端系统那就难说啦。 如何判断B端系统的优劣&#xff0c;各位看官&#xff0c;各抒己见吧。 判断B端系统界面的…