Creating Heatmaps in R

Introduction

This is a brief tutorial on retrieving and preparing data for generating a heatmap in R using genomic data. We’ll be using publicly available data from “Comparative transcriptomic and proteomic analyses provide insights into the key genes involved in high-altitude adaptation in the Tibetan pig” [Scientific Reports 2017]. This is a comparative analysis of the transcriptomic and proteomic profiles of heart tissues obtained from Tibetan and Yorkshire pigs raised at high (TH and YH) and low (TL and YL) altitudes. The study aimed to identify key genes and molecular mechanisms involved in the high-altitude adaptations of the Tibetan pig.

The full text of the study is available through PubMed Central. The raw and processed data is deposited in GEO under series GSE92981.

Load packages

Load the relevant packages.

library(tidyverse)
library(magrittr)
library(pheatmap)
library(RColorBrewer)
library(rio)

If you get an error there is no package, you need to install some packages first. If you install some, others should install automatically as dependencies. Uncomment (remove #) the relevant commands to install.

# install.packages("tidyverse")
# install.packages("pheatmap")
# install.packages("rio")

Retrieve data

Download values

FPKMs (Fragments Per Kilobase of transcript per Million mapped reads) are a common RNA-seq expression unit. First, download the FPKM table. In the GEO submission, this is the supplementary file “GSE92981_expression_all_FPKM.txt.gz”.

fpkm_full = read_tsv("ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE92nnn/GSE92981/suppl/GSE92981_expression_all_FPKM.txt.gz")

Parsed with column specification:
cols(
  .default = col_double(),
  gene_id = col_character(),
  gene_name = col_character(),
  GeneID = col_character(),
  description = col_character(),
  locus = col_character(),
  `TH/YH` = col_character(),
  `TH/TL` = col_character(),
  `TL/YL` = col_character(),
  `YH/YL` = col_character()
)
See spec(...) for full column specifications.

Let’s check the dimensions of the downloaded FPKM table.

dim(fpkm_full)

[1] 18975    29

We should also check the contents of the table. It’s a good idea to use head() when dealing with potentially large tables to prevent excessive output.

head(fpkm_full)

Download stats

We have the FPKMs, but the table contains the full gene list. We probably want to use just an interesting subset. The study also provides a table of differentially expressed genes. It’s not in GEO, but it’s in the supplemental data (Supplementary Table S4). Unfortunately, it’s an Excel file, which makes it slightly more difficult to import. We’ll use the “Table S4-1” sheet (DEGs of TH vs YH).

stats_full = import("https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5473931/bin/41598_2017_3976_MOESM2_ESM.xls", sheet = "Table S4-1", skip = 1)

Check the dimensions of the stats table.

dim(stats_full)

[1] 299   6

Check the contents of the stats table.

head(stats_full)

Prepare data

Clean up values

The FPKM table has a lot of columns. Let’s check just the FPKM table column names.

colnames(fpkm_full)

 [1] "gene_id"     "gene_name"   "GeneID"      "description"
 [5] "locus"       "TH1_FPKM"    "TH2_FPKM"    "TL1_FPKM"   
 [9] "TL2_FPKM"    "YH1_FPKM"    "YH2_FPKM"    "YL1_FPKM"   
[13] "YL2_FPKM"    "TH1_count"   "TH2_count"   "TL1_count"  
[17] "TL2_count"   "YH1_count"   "YH2_count"   "YL1_count"  
[21] "YL2_count"   "TH"          "TL"          "YH"         
[25] "YL"          "TH/YH"       "TH/TL"       "TL/YL"      
[29] "YH/YL"      

The magrittr package adds the ability to use %>% to “pipe” the expression forward, allowing the conversion of nested function calls into a simple pipeline of operations. This makes the code easier to read and work with if you want to perform multiple operations. For example, we can repeat the previous colnames() operation with a pipe.

fpkm_full %>% colnames()

 [1] "gene_id"     "gene_name"   "GeneID"      "description"
 [5] "locus"       "TH1_FPKM"    "TH2_FPKM"    "TL1_FPKM"   
 [9] "TL2_FPKM"    "YH1_FPKM"    "YH2_FPKM"    "YL1_FPKM"   
[13] "YL2_FPKM"    "TH1_count"   "TH2_count"   "TL1_count"  
[17] "TL2_count"   "YH1_count"   "YH2_count"   "YL1_count"  
[21] "YL2_count"   "TH"          "TL"          "YH"         
[25] "YL"          "TH/YH"       "TH/TL"       "TL/YL"      
[29] "YH/YL"      

Let’s clean up the table. The dplyr package has many helper functions to make this easier. We’ll create a new column with both gene names and gene IDs using mutate(), keep only relevant columns with select(), convert to a standard data frame (from tibble), and set gene column as row names.

fpkm_clean = fpkm_full %>%
  mutate(gene = paste0(gene_id, ":", gene_name)) %>%
  select(gene, ends_with("FPKM")) %>%
  as.data.frame() %>%
  column_to_rownames("gene")
head(fpkm_clean)

We can get rid of “FPKM” in column names to keep them simple.

colnames(fpkm_clean) = gsub("_FPKM", "", colnames(fpkm_clean))
head(fpkm_clean)

In addition to the shape of the table, it’s important to check the contents. With thousands of rows, that is hard to do manually. A boxplot can quickly show the distribution.

# "las = 2" makes the sample labels perpendicular
boxplot(fpkm_clean, las = 2)

The range of values is so high, the boxplots are not even visible. We can log-transform the values to make variation more similar across orders of magnitude.

fpkm_log2 = log2(fpkm_clean + 1)
boxplot(fpkm_log2, las = 2)

Clean up stats

As we did with FPKMs, we’ll create a new column with both gene names and gene IDs so that we have that in common.

stats_full = stats_full %>% mutate(gene = paste0(`Gene ID`, ":", `Gene Name`))
head(stats_full)

We can sort by significance.

stats_full %>% arrange(`P-value`) %>% head

Get the list of most significant genes.

top_genes = stats_full %>% arrange(`P-value`) %>% head(50) %$% gene
head(top_genes)

[1] "ENSSSCG00000025268:CKM"    "ENSSSCG00000018092:MT-ND6"
[3] "ENSSSCG00000003431:NPPB"   "ENSSSCG00000003430:NPPA"  
[5] "ENSSSCG00000018084:MT-ND3" "ENSSSCG00000009668:CLU"   

Generate heatmaps

Basic heatmaps

Select 100 random genes to use. We’ll only use these genes for the next few examples to save time.

random_genes = sample(rownames(fpkm_clean), 100)
head(random_genes)

[1] "ENSSSCG00000003168:ALDH16A1" "ENSSSCG00000023857:ARHGEF40"
[3] "ENSSSCG00000002996:RAB4B"    "ENSSSCG00000024394:IKZF4"   
[5] "ENSSSCG00000029664:-"        "ENSSSCG00000003138:HSD17B14"

Generate a simple heatmap using the original values for 100 random genes.

pheatmap(fpkm_clean[random_genes, ])

Generate a simple heatmap using the log-transformed values for 100 random genes.

pheatmap(fpkm_log2[random_genes, ])

Center and scale in the row direction.

pheatmap(fpkm_log2[random_genes, ], scale = "row")

Color options

Manually create color range.

my_colors = c("green", "black", "red")
my_colors

[1] "green" "black" "red"  

Expand the color range from 3 to 50 values.

my_colors = colorRampPalette(my_colors)(50)
my_colors

 [1] "#00FF00" "#00F400" "#00EA00" "#00DF00" "#00D500" "#00CA00"
 [7] "#00C000" "#00B600" "#00AB00" "#00A100" "#009600" "#008C00"
[13] "#008200" "#007700" "#006D00" "#006200" "#005800" "#004E00"
[19] "#004300" "#003900" "#002E00" "#002400" "#001A00" "#000F00"
[25] "#000500" "#050000" "#0F0000" "#1A0000" "#240000" "#2E0000"
[31] "#390000" "#430000" "#4E0000" "#580000" "#620000" "#6D0000"
[37] "#770000" "#820000" "#8C0000" "#960000" "#A10000" "#AB0000"
[43] "#B60000" "#C00000" "#CA0000" "#D50000" "#DF0000" "#EA0000"
[49] "#F40000" "#FF0000"

Generate a heatmap with a green/red color scheme.

pheatmap(fpkm_log2[random_genes, ], scale = "row", color = my_colors)

RColorBrewer provides a collection of pre-made color palettes.

display.brewer.all()

Use an RColorBrewer color palette.

my_colors = brewer.pal(n = 11, name = "RdBu")
my_colors = colorRampPalette(my_colors)(50)
my_colors = rev(my_colors)
my_colors

 [1] "#053061" "#0A3B70" "#10467F" "#16518E" "#1B5C9E" "#2166AC"
 [7] "#2870B1" "#2F79B5" "#3682BA" "#3D8BBF" "#4695C4" "#569FC9"
[13] "#66A9CF" "#76B3D4" "#86BDDA" "#95C6DF" "#A2CDE2" "#AFD4E6"
[19] "#BCDAEA" "#C9E1ED" "#D4E6F0" "#DBEAF2" "#E3EDF3" "#EBF1F4"
[25] "#F3F5F6" "#F7F4F2" "#F8EEE8" "#FAE8DE" "#FBE3D4" "#FCDDCA"
[31] "#FBD4BE" "#FAC9B0" "#F8BEA2" "#F6B394" "#F4A886" "#EF9B7A"
[37] "#E98D6F" "#E37E64" "#DD7059" "#D7624F" "#D05447" "#C84540"
[43] "#C13639" "#BA2832" "#B2192B" "#A41328" "#940E26" "#850923"
[49] "#760421" "#67001F"

Generate a heatmap with a blue/red color scheme.

pheatmap(fpkm_log2[random_genes, ], scale = "row", color = my_colors)

Significant genes

Visualize the most significant genes.

pheatmap(fpkm_log2[top_genes, ], scale = "row", color = my_colors)

Remove the border around each cell and reduce the row label font size.

pheatmap(fpkm_log2[top_genes, ],
         scale = "row",
         color = my_colors,
         border_color = NA,
         fontsize_row = 4)

Additional notes

This is just one approach to generating a heatmap. There are many other ways. You can even do it in Excel (sort of). Check this listing of graphical (no coding needed) and R-based alternatives: http://bit.ly/heatmapoptions.