关于南丁格尔图的“绘后感”

原创

HenryLei

修改于 2024-05-10 14:28:59

2490

修改于 2024-05-10 14:28:59

文章被收录于专栏：生信笔记生信笔记

前言

先说点跑题的内容。

通过绘图，体会到，数据分析的核心在于数据整理。数据整理好了，画图是几行代码的事儿；数据没有整理好，再厉害的码农也画不出想要的图。
关于数据整理，原则是根据你的呈现目标整理&根据R语言函数对数据的要求整理。即既要满足想要呈现的内容又要满足代码对输入数据的要求。因此，要用计算机语言的思考方式，根据自己的目标整理数据。
不同的数据整理的方式会有不同。即使作相同的图，也没法完全照套相同图形的代码。即“一图一码”。

再说点其他跑题的内容。

不久前，我同学委托我帮助其画图，于是给了我如下的样图，让我照着画。

但我画出来是这样的：

虽然有点差距，但作为新手，我自己已经很满意了。所以今天就总结一下，温故知新。

数据的准备

收到的数据是一个Excel表，通常大家会按照下面的形式进行分类整理数据。

但是，准确的说，上面这种数据排布形式只是方便填写和阅读，并不能用于作为R语言的输入数据的排布形式。因此，我们需要按照计算机语言能够理解的思维方式重新整理数据。

导入R前的数据整理

一、数据整理的原则

我自己总结的原则是，如果你画的是二维图，即只有X和Y轴的图，那么你的数据需要整理成核心只有两列的数据表。假设你画的是三维图（当然，我没有画过，暂时这样预设），即有X、Y、Z轴，那么我想，你需要将数据整理成核心有三列的数据表。

二、什么叫“核心只有两列的数据表”

因为最终画出来的图只有x和y轴，无论你将数据分了多少组，将样本分了多少组，即你要做多少种标记（颜色、形状等等）或者你重复测了多少次，有多少平行数据等等，图像要表现的关系核心，只有x和y轴体现出来的数据。

只要你的数据需要画在一张图上，那么除了核心数据，其他都需要看成分组。

这样，我们需要将x轴的数据整理成1列，将y轴的数据整理成1列，将各种分组的方式，按照需要整理的若干列，与x和y列的数据对应起来即可。

三、针对上述表格的具体操作

有了上面叙述的原则，我们尝试将原始获得的表格进行整理。

在上面的表格中，我们需要表现的是微生物种名和两种方式的值之间的关系。即，微生物种名和值分别是x和y轴表现的数据，两种方式测得的值是“值”的分类，真菌、病毒、细菌是“微生物种名”的分类。

因此我们将tNGS和mNGS合并成1列，增加1列“值的分类”，对应数据的单元格内标上对应的tNGS和mNGS。另外增加1列“名称的分类”，与物种名称对应填上真菌、病毒和细菌。

由于数据比较少，也比较简单，上面这些前期整理的步骤我就在Excel表里直接用鼠标拖动几下就完成了。当然你也可以导入R里面，用函数进行处理也是可以的。整理完成后，另存为成.CSV格式的文件，便于R读取。

导入R后的预处理

一、导入数据

data <- read.csv("raw_data.csv", header = T)
dim(data)
[1] 52  4
data.clean <- data[,1:4]

这里第一次导入的时候还有一个小插曲，我用dim查看的时候，发现有5列，于是点进表格查看，发现多了一列空列x，可能是由于在保存csv文件的时候，Excel表的一个空列被认为做过修改，所以也作为空列导入了，于是强迫症的我还把空列删除了一下。

上面这里导入是正常的4列，是因为后面在处理数据的时候，发现这个数据里面还有一个坑，我用代码调整了半天，发现还不如直接在Excel表里面整理来得快，于是上面的csv文件是后续修改过的，这个后面再细说。

二、基于绘图要求的数据整理思路

由于要画成“南丁格尔图”，我查了一下，普遍的画法是将柱状图再加一层极坐标的图层就可以实现旋转。但这里的问题重点在文字标注。

如果柱状图带着X轴的刻度标签添加极坐标图层，X轴的标签是不旋转的。即原来是水平方向放在X轴的下方，添加极坐标后，标签依然水平围绕着极坐标。如下图（这里是调整成了90度方向）
如果在柱状图上方标注好文字信息，再添加极坐标旋转，同样文字不会旋转。

因此，如果需要画成像文章开始那样的文字围绕图形旋转的样式，只能图形和文字分别在2个图层中，各自按照角度旋转，再匹配上。

三、基于以上思路需要准备的数据

图形数据，即是表格中已有的数据（Counts值）
分组数据和物种分类数据，表格中也有（Target和Classification）
文字标签的数据，包括
1. 各物种的名称
2. 名称的旋转角度
3. 物种类别名称
4. 类别名称的旋转角度

需要说明的是，物种名称我们可以用Species列的数据，但是你会发现每个名称有2个重复，如果用这个数据，那么标签文字就会有重复。同样，物种类别也可以用Classification列，但同样会有重复。因此，我们需要单独准备去重后的Species和Classification。

旋转角度，我们可以设想将360度按照去重后Species的数量进行平均分配，那么每个分配到的度数就是每个标签旋转的角度。

Classification标签文字的旋转相对简单一些，因为只有3种，我们可以根据Species标签画好后的具体位置进行手动设置。

（一）去重Species与旋转角度的数据准备

由于最终需要按照物种所属的3个类别集中在一起呈现，因此最终x轴的物种顺序应该与上图表格中的顺序一致（或者Fungus，Virus与Bacterium任意的前后顺序），但是在这种情况下，对Species列去重后，由于每一类的重复数量不同，对应生成的新列会稍微复杂一点（也可以生成）。因此，我先将Species列按照字母排序后，再进行去重。

1.去重Species数据准备

> table(dt.cl.resorted$Classification)#可见Classification重复数量不同

Bacterium    Fungus     Virus 
       30        10        12 
library(tibble)#rownames_to_column()函数在tibble包中
data.clean <- rownames_to_column(data.clean,var = "ID")#由于没有指定行名，所以行名是序号。将行名转换成列，便于后续再将表格排序还原。
> data.clean %>% head()
  ID Classification                Species Counts Target
1  1         Fungus       Candida albicans      6   tNGS
2  2         Fungus  Aspergillus fumigatus      2   tNGS
3  3         Fungus Pneumocystis jirovecii      3   tNGS
4  4         Fungus   Candida dubliniensis      0   tNGS
5  5         Fungus      Candida nigricans      0   tNGS
6  6         Fungus       Candida albicans      5   mNGS

#按种名排序
data.clean.sorted <- data.clean[order(data.clean$Species),]#order函数排序，返回的是排序后的行号；sort函数排序，返回的是排序后的内容。因此这里用排序后的行号来调整整个表格的顺序 
ID.sorted <- seq(1,nrow(data.clean.sorted))
data.clean.sorted$ID.sorted <- ID.sorted#给排序后的表格也添加一个序号
> data.clean.sorted %>% head()
   ID Classification                 Species Counts Target ID.sorted            
24 24      Bacterium Acinetobacter baumannii      6   tNGS         1
39 39      Bacterium Acinetobacter baumannii      6   mNGS         2
2   2         Fungus   Aspergillus fumigatus      2   tNGS         3
7   7         Fungus   Aspergillus fumigatus      1   mNGS         4
1   1         Fungus        Candida albicans      6   tNGS         5
6   6         Fungus        Candida albicans      5   mNGS         6

#种名去重
data.clean.sorted$uniq.species <- NA #这里首先新增1列，赋为缺失值（NA）
data.clean.sorted$uniq.species[seq(1,nrow(data.clean.sorted),2)] <- unique(data.clean.sorted$Species)#再将去重后的种名，间隔写入新增列中

去重后的种名需要编号，以便后续用于设置旋转角度，但是我在这里踩了坑，直接在这里编号了。

#给去重后的种名编号
data.clean.sorted$uniq.ID <- NA #同样先增加1列缺失值
data.clean.sorted$uniq.ID[seq(1,nrow(data.clean.sorted),2)] <- seq(1,26) #间隔写入编号。理论上这里26最好用函数nrow(data.clean.sorted)/2代替。这里总行数是偶数，若是奇数，应该怎么写函数呢？
> data.clean.sorted %>% head()
   ID Classification                 Species Counts Target ID.sorted            uniq.species uniq.ID
24 24      Bacterium Acinetobacter baumannii      6   tNGS         1 Acinetobacter baumannii       1
39 39      Bacterium Acinetobacter baumannii      6   mNGS         2                    <NA>      NA
2   2         Fungus   Aspergillus fumigatus      2   tNGS         3   Aspergillus fumigatus       2
7   7         Fungus   Aspergillus fumigatus      1   mNGS         4                    <NA>      NA
1   1         Fungus        Candida albicans      6   tNGS         5        Candida albicans       3
6   6         Fungus        Candida albicans      5   mNGS         6                    <NA>      NA

然后生成原顺序表，发现不对。。。

#顺序表
dt.cl.resorted <- data.clean.sorted[order(as.numeric(data.clean.sorted$ID)),]
> data.clean.sorted[order(as.numeric(data.clean.sorted$ID)),]
   ID Classification                    Species Counts Target ID.sorted               uniq.species uniq.ID
1   1         Fungus           Candida albicans      6   tNGS         5           Candida albicans       3
2   2         Fungus      Aspergillus fumigatus      2   tNGS         3      Aspergillus fumigatus       2
3   3         Fungus     Pneumocystis jirovecii      3   tNGS        41     Pneumocystis jirovecii      21
4   4         Fungus       Candida dubliniensis      0   tNGS         7       Candida dubliniensis       4
5   5         Fungus          Candida nigricans      0   tNGS         9          Candida nigricans       5
6   6         Fungus           Candida albicans      5   mNGS         6                       <NA>      NA
7   7         Fungus      Aspergillus fumigatus      1   mNGS         4                       <NA>      NA
8   8         Fungus     Pneumocystis jirovecii      1   mNGS        42                       <NA>      NA
9   9         Fungus       Candida dubliniensis      1   mNGS         8                       <NA>      NA
10 10         Fungus          Candida nigricans      2   mNGS        10                       <NA>      NA

uniq.ID应该在恢复原顺序之后，按照从1到26的顺序编号，如今这样就乱套了。。。

正确应该是，上表中，uniq.ID为NA，然后根据uniq.species列对应的非NA行填入顺序编号1到26，于是我重新编号。

#修改uniq.id
uniq.id <- 1
for (i in 1:52){ #理论上，这里52也应该写成nrow(dt.cl.resorted),我这里仗着数据表内容比较少，就偷懒。。。
  if (is.na(dt.cl.resorted$uniq.species[i])){
    next
  }else{
    dt.cl.resorted$uniq.ID[i] <- uniq.id
    uniq.id=uniq.id+1 #本来变量自增，我这里写的是 uniq.id += 1,发现不行。似乎R里面不能这样写
  }
}

> head(dt.cl.resorted, 10)
   ID Classification                Species Counts Target ID.sorted           uniq.species uniq.ID
1   1         Fungus       Candida albicans      6   tNGS         5       Candida albicans       1
2   2         Fungus  Aspergillus fumigatus      2   tNGS         3  Aspergillus fumigatus       2
3   3         Fungus Pneumocystis jirovecii      3   tNGS        41 Pneumocystis jirovecii       3
4   4         Fungus   Candida dubliniensis      0   tNGS         7   Candida dubliniensis       4
5   5         Fungus      Candida nigricans      0   tNGS         9      Candida nigricans       5
6   6         Fungus       Candida albicans      5   mNGS         6                   <NA>      NA
7   7         Fungus  Aspergillus fumigatus      1   mNGS         4                   <NA>      NA
8   8         Fungus Pneumocystis jirovecii      1   mNGS        42                   <NA>      NA
9   9         Fungus   Candida dubliniensis      1   mNGS         8                   <NA>      NA
10 10         Fungus      Candida nigricans      2   mNGS        10                   <NA>      NA

2.旋转角度的数据准备

#角度
angel <- 90-360*(dt.cl.resorted$uniq.ID-0.5)/(nrow(dt.cl.resorted)/2) #将360度分成26份，按照uniq.ID的份数分配转角度数。但由于我们的图形是双柱状图，为了让文字在两个柱体中间，所以每个份额让出0.5。由于从90度开始，因此前面有90-
dt.cl.resorted$angel <- ifelse(angel < -90, angel + 180, angel) #angle指定的度数是以对象的中心点作旋转的。因此按照angel的数值旋转后，左半圆的字体会倒置。为了避免这个现象，将左半圆的度数+180（相当于以自己本身中心点旋转180，或者说理解成左右半圆的度数对称相等）
dt.cl.resorted$angel2 <- angel #保留一列非对称旋转的数据，后面会用到，再解释

默认情况下，极坐标函数coord_polar作图的原点是以时钟的0点开始的。但在ggplot2中的各图层函数的angle参数（设置旋转度数）的值是以直角坐标系为参照，以角度为单位。因此时钟的0点对应直角坐标系的90°

（二）外圈文字准备

即Classification的文字

#外圈文字准备
dt.cl.resorted$uniq.class <- NA
dt.cl.resorted$uniq.class[c(3,13,31)] <- c("Fungus","Virus","Bacterium") #3个文字放置的位置，我是根据上面的图画好后，对比着与哪个Species的名称对应好看，即放置在哪个位置

#文字旋转度数
dt.cl.resorted$class.angel <- NA
dt.cl.resorted$class.angel[c(3,13,31)] <- c(-27.6931,-135.002,90)#参照上图中Species的旋转角度，手动设置。当然也可以按照每类中Species数量的多少，按照比例瓜分360度来设置，类似与上面的angel

四、基于函数要求的数据处理

#分类变量映射的因子化
Groups <- factor(dt.cl.resorted$Target)  #映射中，分类变量一定要用因子，这是函数本身的要求
#自定义x轴顺序的因子化处理
dt.cl.resorted$Species <- factor(dt.cl.resorted$Species,levels=dt.cl.resorted$uniq.species) #因子水平不能有重复。必须与变量中的值对应，因子水平中没有的变量会被设置成缺失值(NA)

关于x轴的顺序。由于本次数据x轴本身也是分类变量，理论上也要先因子化，才能进行映射画图。但是画柱状图的时候，默认会将x轴的分类变量自动因子化然后作图。自动因子化的时候，因子水平按照字母顺序排列，因此作图后x轴的顺序是字母顺序。因此需要手动指定因子水平的顺序。

关于因子
因子相当于是给分类变量设置顺序。即因子水平中指定的顺序即为分类变量的顺序。这与分类变量本身在向量中的排列顺序无关。

作图

柱状图

先作柱状图

library(ggplot2)
library(ggprism)
library(ggthemes)

p <- ggplot(
  dt.cl.resorted, #由于刚开始我不停调整其他图层的映射变量，所以主函数这里只指定了data，后面看来所有图层的x和y轴都是相同映射
)+
  geom_bar(
    aes(x = Species,
        y = Counts,
        fill = Groups),#用因子化的Target列分类变量作为填充柱状图的颜色分类
    position="dodge", #并排分类变量Groups（Target），默认是stack，即上下堆叠。jitter是左右重叠堆放
    stat = "identity", #指定y轴数据是统计个数还是直接用数值。默认count，即系将Y轴数据作为分类变量统计个数
    show.legend = T, #是否展示图例
  )+
  theme_minimal()+
  theme(
    axis.text = element_blank(),  #坐标轴标签设置为空白，因为后面自己设置    
    axis.title = element_blank(), #坐标轴标题，这里设置为空白
    #panel.grid = element_blank(), 可以取消图型背景的格子线
    plot.margin = unit(rep(1,4),"cm") #边界
  )+
  scale_fill_prism(palette = "candy_bright")+
  scale_fill_manual(values = c('#fec79e','#8ec4cb')) #填充颜色

p

南丁格尔图

再作极坐标图（南丁格尔图），并标注文字

pp <- p+
  coord_polar()+  #转换成极坐标
  ylim(c(-5,18))+ #用y轴的极值范围来设置圆心的留白大小和外围圆圈的范围
  geom_text(
    aes(
      x=Species,
      y=Counts,
      label=uniq.species, #制定文字内容，由于映射全部的数据，所以标签要用包含NA值的数据，从而一一对应
      color=Classification #字体颜色用物种类别区分
    ),
    size = 3,
    angle=dt.cl.resorted$angel, #旋转角度
    hjust=ifelse(dt.cl.resorted$angel2 <  -90, 1, 0),#调整文字与柱状图的相对位置，默认写到直方图上面。0到1的值。1表示右适应，0表示左适应。这是为什么前面保留angel原始角度的原因：在这里用于判断。    
  )+
guides(color=F) #单独移除名称分类图例（即本层图层的图例）
  
pp

外圈文字添加

p3 <- pp+
  geom_text(
    aes(
      x=Species,
      y=Counts,
      label=uniq.class,
      color=uniq.class
    ),
    angle=dt.cl.resorted$class.angel,
    hjust=0.5,
    vjust=-12,
    frontface="bold",
    size=5,
    )

p3

保存

dev.off()
png("pix_draft.png", 
    width = 2800, 
    height = 1800, 
    res = 300)
pp
dev.off()

png("draft.png",
    width = 2800,
    height = 1800,
    res = 300)
p3
dev.off()