Heatmap with star

(Carolina email, Dec 30th 2021- help with heatmap plotting)

code and data : /space/chen-syn01/1/data/cinliu/plots/heatmaps/heatmap_with_star

Goal plot:

Original instructions

I need to generate a heatmap with the genetic correlation data in the *.txt file attached. I am looking to generate something like the attached screenshot, where we plot our 24 brain phenotypes ("Region") along the y-axis, and the "Disorder" on the x-axis. For the color of the squares, we can plot the "rg" value, and if we could put a single * over the squares that have p<0.05.

For ordering of the variables, could we order the brain phenotypes as we have in previous figures for the paper - I have this in one of my codes:

reorder_pheno <-c("1.motorpremotor_area","2.dorsolateralprefrontal_area","3.dorsomedialfrontal_area","4.orbitofrontal_area",

"5.parsopercularis_area","6.superiortemporal_area","7.posterolateraltemporal_area","8.anteromedialtemporal_area",

"9.inferiorparietal_area","10.superiorparietal_area","11.precuneus_area","12.occipital_area","1.motor_premotor_SMA_thickness",

"2.superiorparietal_thickness","3.inferiorparietal_thickness","5.occipital_thickness","4.perisylvian_thickness","8.temporalpole_thickness", "9.medialtemporal_thickness" ,"10.middletemporal_thickness" , "6.ventromedialoccipital_thickness",

"11.dorsolateralprefrontal_thickness","7.ventralfrontal_thickness","12.medialprefrontal_thickness")

And for the disorders, if we could order them this way (the first eight are neurodevelopmental, the last 6 are neurodegenerative):

reorder_disorder<-c("ADHD","AN","ASD","BPD","MDD","OCD","PTSD","SZ","AD","ALS","EP","FTD","PD","PSP)

Basic steps

Setup

rm(list=ls())

library(graph4lg)

library(reshape2)

#dir <- "/space/chen-syn01/1/data/chen/gwas_gclust/genesis/ABCD/"

dir <- "/Users/nini/Desktop/2020lab/plots/CH/" #UPDATE

setwd(dir)

pheno_dir <- "/space/chen-syn01/1/data/cmakowski/ABCD/GCLUST/GWAS/Science_revisions/"

Import data

import the .txt file

dat <- read.delim("/Users/nini/Desktop/2020lab/carolina/heatmap/brainregion_disorder_ldsc.txt") #update

2. Manipulating data

Reordering the dataset.

####REORDERING####

reorder_pheno <-c("1.motorpremotor_area","2.dorsolateralprefrontal_area","3.dorsomedialfrontal_area","4.orbitofrontal_area", "5.parsopercularis_area","6.superiortemporal_area","7.posterolateraltemporal_area","8.anteromedialtemporal_area", "9.inferiorparietal_area","10.superiorparietal_area","11.precuneus_area","12.occipital_area","1.motor_premotor_SMA_thickness", "2.superiorparietal_thickness","3.inferiorparietal_thickness","5.occipital_thickness","4.perisylvian_thickness", "8.temporalpole_thickness", "9.medialtemporal_thickness" ,"10.middletemporal_thickness" , "6.ventromedialoccipital_thickness", "11.dorsolateralprefrontal_thickness","7.ventralfrontal_thickness","12.medialprefrontal_thickness") reorder_disorder <- c("ADHD","AN","ASD","BPD","MDD","OCD","PTSD","SZ","AD","ALS","EP","FTD","PD","PSP")

reorder_pheno <- sub(".*\\.","",reorder_pheno) #remove number and dot

dat$Region <- as.factor(dat$Region)

#change the differently spelled words to match the dataset spelling

setdiff(reorder_pheno, dat$Region) #looking for mismatch

setdiff(reorder_disorder, dat$Disorder)#looking for mismatch

n <- grep("motor_premotor_SMA_thickness", reorder_pheno) #change "motor_premotor_SMA_thickness" to "motor_thickness"

reorder_pheno[n] <- "motor_thickness"

n <- grep("SZ", reorder_disorder) #change "SZ" to "SCZ"

reorder_disorder[n] <- "SCZ"

dat$Disorder <- factor(dat$Disorder, levels =reorder_disorder)

dat$Region <- factor(dat$Region, levels = reorder_pheno)

3. Plotting

ggplot() is used for plotting (in particular geom_tile()).

geom_tile(aes(fill=rg),colour= "black",size=0.4)
geom_point(data = dat[dat$p < 0.05, ],shape = "*", size=4) #makes the *

plot1 <- ggplot(dat, aes(Disorder,Region)) +

geom_tile(aes(fill=rg),

colour= "black",

size=0.4) +

geom_point(data = dat[dat$p < 0.05, ],shape = "*", size=4) +

scale_x_discrete(position = "top")+

scale_y_discrete(limits=rev) +

theme(panel.background = element_blank(),

axis.title.x=element_blank(),

axis.title.y=element_blank(),

axis.ticks = element_blank(),

text = element_text(size=12),

legend.key.size = unit(1.9, "cm"),

legend.key.width = unit(0.4,"cm"),

legend.key = element_rect(colour="black"),

legend.margin = margin(t = -14, r = 0, b = 3, l = 0, unit = "pt"),

axis.text.x = element_text(angle = 90, vjust = 0, hjust=0)) +

scale_fill_gradient2(low = "#0288d1",mid = "white", high= "#e64a19")

plot1

ggsave(file="heatmap.png", plot1, width = 5.5, height = 5.4)

Full code

heatmap.R

#heatmap for carolina