写在前面

之前曾抱怨过complexheatmap太复杂了，现在我想说是我太菜了。折腾了一圈complexheatmap后，发现complexheatmap真的是极致全面的一个热图包了。

TidyTuesday

如果你想入门R或者python，但是没有数据给你折腾而发愁。那么我推荐你去看看TidyTueday。点这里

TT

每周二一个新的数据包，覆盖各种类型。比如今天我们折腾的数据就是来自TidyTuesday。想要获取数据非常简单，只需要安装TidyTuesday的包。

install.packages("tidytuesdayR")

首先读入需要的包。

#Loading necessary packages
library(tidytuesdayR)
library(tidyverse)
library(DataExplorer) # easy way for data exploring

然后

# Import data from tidytuesday
tuesdata <- tidytuesdayR::tt_load('2020-09-01')
# or
tuesdata <- tidytuesdayR::tt_load(2020, week = 36)
--- Compiling #TidyTuesday Information for 2020-09-01 ----
--- There are 5 files available ---
--- Starting Download ---

    Downloading file 1 of 5: `arable_land_pin.csv`
    Downloading file 2 of 5: `cereal_crop_yield_vs_fertilizer_application.csv`
    Downloading file 3 of 5: `cereal_yields_vs_tractor_inputs_in_agriculture.csv`
    Downloading file 4 of 5: `key_crop_yields.csv`
    Downloading file 5 of 5: `land_use_vs_yield_change_in_cereal_production.csv`

--- Download complete ---

当然，具体请考虑自己的网速。因为数据保存在github上……

数据可视化

key_crop_yields <- tuesdata$key_crop_yields

key_crop_yields %>%
  filter(Code != ""& Entity != "World") %>% 
  select(-(2)) %>% 
  create_report(
    config = configure_report(
      add_plot_bar = FALSE
    )
  )

这次使用DataExplorer自动生成报告，当然你可以通过函数查看更具体的内容。

key_crop_yields %>%
  filter(
    Entity %in% c(
      "Japan",
      "China",
      "United States",
      "Brazil",
      "Indonesia",
      "India",
      "Pakistan",
      "Nigeria",
      "Bangladesh",
      "Russia",
      "Mexico"
    )
  ) %>%
  select(-(2)) %>%
  mutate(Year = as.factor(Year),
         Entity = as.factor(Entity)) %>% plot_boxplot(by = c("Year"))
  

可可的产量很有意思啊。
或者

key_crop_yields %>%
  drop_na() %>% 
  filter(Code != "" & Entity != "World" & Year == 2018) %>%
  select(-(1:3)) %>% 
  plot_correlation(type = "c")

当然，也能用ggplot2

key_crop_yields %>%
  filter(
    Entity %in% c(
      "Japan",
      "China",
      "United States",
      "Brazil",
      "Indonesia",
      "India",
      "Pakistan",
      "Nigeria",
      "Bangladesh",
      "Russia",
      "Mexico"
    )
  ) %>%
  ggplot(., aes(x = Year, y = `Wheat (tonnes per hectare)`, color = Entity)) +
  geom_point(size = .1) +
  geom_smooth(method = 'lm') +
  theme_bw()

尼日利亚的小麦产量居然是下降的，听说他们主粮是一种豆子。中国小麦公顷产增长超级快，那么中国小麦公顷产达到世界前列了么？一张热图告诉你。
注：这里选择的是人口前11的国家。

世界小麦公顷产量热图

hmmm 看来不是…… 什么你问热图怎么做的？complexheatmap 啊！
代码具体参考kevinblighe/E-MTAB-6141: Data from Lewis, Barnes, Blighe et al., Cell Rep. 2019 Aug 27; 28(9): 2455–2470.e5. (github.com)

#载入包
require(RColorBrewer)
require(ComplexHeatmap)
require(circlize)
require(cluster)
require(countrycode)

然后整理一下数据

set.seed(2021)
yield_heatmap <- key_crop_yields %>%
  filter(!is.na(Code)) %>%
  select(Entity, Year, `Wheat (tonnes per hectare)`) %>%
  spread(Year, `Wheat (tonnes per hectare)`) %>% 
  filter(Entity != "World") %>% 
  column_to_rownames(var = "Entity") %>% 
  drop_na()

key_crop_yields %>% 
  filter(str_detect(Entity,"Sudan")) %>% 
  select(Entity, Year, `Wheat (tonnes per hectare)`) %>% 
  drop_na() %>% 
  mutate(Entity = "Sudan") %>% 
  arrange(Year) %>% 
  spread(Year, `Wheat (tonnes per hectare)`) %>% 
  column_to_rownames(var = "Entity") %>% 
  bind_rows(yield_heatmap) -> heat

上面做了两件事，一个是把数据里各种大洲和世界数据去掉了，另一个就是把苏丹数据跟前苏丹的数据合并了。
注：这里指的是Entity不是Country，防杠还是去掉了一些有争议的地区。

然后joint一下地理信息，后来发现起始直接通过Code就能获得准确信息了……以后说吧。

rownames(heat) %>% data.frame(country = .) -> df
df$Continent <-  countrycode(sourcevar = df[, "country"],
                             origin = "country.name",
                             destination = "continent")

接下来是定义五大洲板块的色块注释。

ann <- data.frame(Continent = df$Continent,
                  row.names = df$country,
                  stringsAsFactors =  FALSE)
colours <- list(
  Continent = c('Africa' = '#2171B5', 'Asia' = '#6A51A3', 'Europe' = 'red2', 'Americas' = 'green3', 'Oceania' = 'orange')
)

colAnn <- HeatmapAnnotation(
  df = ann,
  which = 'row', # 'col' (samples) or 'row' (gene) annotation?
  col = colours,
  annotation_height = 0.6,
  annotation_width = unit(1, 'cm'),
  gap = unit(1, 'mm'),
  annotation_legend_param = list(
    Continent = list(
      nrow = 5,
      title = 'Continent',
      title_position = 'topcenter',
      legend_direction = 'vertical',
      title_gp = gpar(fontsize = 12, fontface = 'bold'),
      labels_gp = gpar(fontsize = 12, fontface = 'bold'))
  )
)

顶部的箱线图注释

pick.col <- brewer.pal(5, 'Greens')
boxplotCol <- HeatmapAnnotation(
    Average = anno_boxplot(
      heat,
      border = FALSE,
      gp = gpar(fill = colorRampPalette(pick.col)(length(1:58))),
      pch = '.',
      size = unit(2, 'mm'),
      axis = TRUE,
      axis_param = list(
        gp = gpar(fontsize = 12),
        side = 'left')),
      annotation_width = unit(c(2.0), 'cm'),
      which = 'col')

然后是顶部线图的注释

heat_anno <- HeatmapAnnotation(
    Average = anno_lines(
      apply(heat, 2, mean),
      border = TRUE,
      gp = gpar(col = "red2"),
      size = unit(2, 'mm'),
      axis = TRUE,
      axis_param = list(
        gp = gpar(fontsize = 10),
        side = 'left')),
    Distribution = anno_boxplot(
      heat,
      border = FALSE,
      gp = gpar(fill = colorRampPalette(pick.col)(length(1:58))),
      pch = '.',
      size = unit(2, 'mm'),
      axis = TRUE,
      axis_param = list(
        gp = gpar(fontsize = 12),
        side = 'left')),
      annotation_width = unit(c(2.0), 'cm'),
      which = 'col')

当然你也可以换成柱状图

barplotCol <- HeatmapAnnotation(
    Average = anno_barplot(
      apply(heat, 2, mean),
      border = TRUE,
      gp = gpar(fill = colorRampPalette(pick.col)(length(1:58))),
      size = unit(2, 'mm'),
      axis = TRUE,
      axis_param = list(
        gp = gpar(fontsize = 12),
        side = 'left')),
      annotation_width = unit(c(2.0), 'cm'),
      which = 'col')

然后开始画图热图吧
说个题外话，如果你对魔改complexheatmap有兴趣，强烈建议看看这儿
ComplexHeatmap Complete Reference (jokergoo.github.io)
虽然看的人想喊师傅别念了别念了……

myCol <- colorRampPalette(c('dodgerblue', 'black', 'yellow'))(100) #三色
myBreaks <- seq(1, 8, length.out = 100)  #设置色彩范围
pamClusters <- cluster::pam(heat, k = 5) #设置聚类数

Heatmap(
  heat,
  name = 'Wheat
  (tonnes per hectare)', 
  col = colorRamp2(myBreaks, myCol), #这个颜色设定值得学习，做出来很高级
  row_split = factor(pamClusters$clustering, levels = c(5,4,3,1,2)),  #聚类排序，
  cluster_row_slices = FALSE, #不解释，自己看manual
  heatmap_legend_param = list( #设置标签
    color_bar = 'continuous',
    legend_direction = 'vertical',
    legend_width = unit(8, 'cm'),
    legend_height = unit(5.0, 'cm'),
    title_position = 'topcenter',
    title_gp = gpar(fontsize = 12, fontface = 'bold'),
    labels_gp = gpar(fontsize = 12, fontface = 'bold')
  ),
  # row parameters
  cluster_rows = TRUE, #聚类
  show_row_dend = TRUE, #显示聚类
  row_title = "Cluster",
  row_title_gp = gpar(fontsize = 10,  fontface = 'bold'), #字号，粗体
  row_title_rot = 90, 
  show_row_names = TRUE,
  row_names_gp = gpar(fontsize = 10, fontface = 'bold'),
  row_names_side = 'right',
  row_dend_width = unit(25, 'mm'),
  
  # column parameters
  cluster_columns = FALSE,
  column_title = 'Year',
  column_title_side = 'bottom',
  column_title_gp = gpar(fontsize = 10, fontface = 'bold'),
  column_title_rot = 0,
  show_column_names = TRUE,
  column_names_gp = gpar(fontsize = 10, fontface = 'bold'),
  column_names_max_height = unit(10, 'cm'),
  column_dend_height = unit(25, 'mm'),
  # cluster methods for rows and columns
  clustering_distance_rows = function(m) dist(m),
  row_dend_reorder = T,
  clustering_method_rows = 'ward.D2',
  top_annotation = heat_anno,
  right_annotation = colAnn
) -> heat_map
heat_map

是不是很简单呢！做完我都不记得那些参数是干嘛的了……
那么有没有简单的办法呢！有个低配版，效果嘛……

yield_mean <- apply(heat, 2, mean) %>% data.frame(average = .) %>% mutate(year = rownames(.))

heat %>%
  mutate(country = rownames(.)) %>%
  left_join(., df, by = "country") %>%
  gather(year, yield, -country, -Continent ) %>%
  left_join(., yield_mean, by = "year") %>%
  as.tibble() %>%
  tidyHeatmap::heatmap(country, year, yield,
                       cluster_columns = FALSE,
                       .scale = "none",
                       transform = NULL) %>%
  tidyHeatmap::add_tile(Continent) %>%
  tidyHeatmap::add_bar(average)

怎么说呢，又不是不能用！