reproducibility_report.Rmd

---
title: "Genetic variants associated with alcohol dependence co-ordinate regulation of *ADH* genes in gastrointestinal and adipose tissues"
author: "Rebecca Hibberd and Evgeniia Golovina"
date: "23/01/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message=FALSE)

# install.packages("pacman")
# load libraries
pacman::p_load(ggplot2, ggpubr, viridis, colormap, scales, reshape, dplyr, tidyr, reshape2,
               ComplexHeatmap, tibble, circlize)
```

This is a reproducibility report for "Genetic variants associated with alcohol dependence co-ordinate regulation of *ADH* genes in gastrointestinal and adipose tissues" study.  

Python (version 2.7.15), R (version 3.5.2) and RStudio (version 1.1.463) were used for data processing, analysis and visualisation.  

1. Cell type- and tissue-specific Hi-C data is available on [GEO](https://www.ncbi.nlm.nih.gov/geo/) database (accessions: GSE63525, GSE35156, GSE43070, GSE77565, GSE105194, GSE105513, GSE105544, GSE52457, GSE105914, GSE105957, GSE87112).  
2. RNA-seq and genotyping data (GTEx v7) are available via [dbGaP](https://www.ncbi.nlm.nih.gov/gap/) access (accession: phs000424.v7.p2).  
3. Human genome build hg19 release 75 (GRCh37) (“Homo_sapiens.GRCh37.75.dna.primary_assembly.fa.gz”) was downloaded from ftp://ftp.ensembl.org/pub/release-75/fasta/homo_sapiens/dna/.  
4. SNP genomic positions were obtained from ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/.  
5. Gene annotation for GENCODE v19 (“gencode.v19.transcripts.patched_contigs.gtf”) was downloaded from https://storage.googleapis.com/gtex_analysis_v7/reference/.  
6. SNPs associated with alcohol dependence were downloaded from the [GWAS Catalog](www.ebi.ac.uk/gwas/) on 07/12/2018.  
7. Promoter-enhancer histone marks and affected regulatory motifs were obtained from the [HaploReg](https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php) (version 4.1).  

### 1. Identification of significant spatial eQTL SNP-gene interactions using CoDeS3D.

The [CoDeS3D](https://github.com/Genome3d/codes3d-v1) pipeline was used to identify tissue-specific spatial interactions between regulatory regions (marked by SNPs associated with alcohol dependence) and their target genes.  

First, we run `python codes3d.py -i data/snps/73_gwas_alcohol_dependence_snps_2018-11-30_1E-06.txt -o results/codes3d_output/` to get significant tissue-specific spatial regulatory interactions.

### 2. Functional annotation of the eQTL SNPs.

Functional annotation of SNPs was performed using [wANNOVAR](http://wannovar.wglab.org/) tool.

Locations and functional annotations of eQTLs were plotted across the human genome (hg19):

i) extracted eQTL SNPs, position and chromosome from `results/codes3d_output/AD_significant_eqtls.txt`. Removed duplicates. 
ii) Combined 13 SNPs on Chromosome 4 (positions between 100214164 and 100301126) as a single entry due to close proximity. Position was the calculated average (mean) of all of their positions and assigned functional annotation as "multiple SNPS"
iii) Assigned numbers to different types of functional annotation as their "Type": 1=intergenic, 2= intronic, 3=exonic, 4=ncRNA intronic 5=multiple SNPs
This produced .txt file `results/functional_annotation/eQTL_SNP_grouped.txt`
iv) Obtained chromosome sizes for hg19 - `results/functional_annotation/chr_sizes_hg19.txt`

```{r functional_annotation}
chrdata <- read.table("results/functional_annotation/chr_sizes_hg19.txt", sep = "\t", header=TRUE)
eQTLpos<-read.table("results/functional_annotation/eQTL_SNP_grouped.txt", sep = "\t", header=TRUE)

chrdata$Chromosome <- factor(chrdata$Chromosome, levels = chrdata$Chromosome)

#tiff("eQTL_SNP_map.tiff", units="in", width=12, height=6, res=300)
ggplot(chrdata, aes(Chromosome,Position), legend(colors())) +
  geom_bar(stat="identity", fill = "grey70", width = 0.5) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.text.x = element_text(vjust=9, hjust=0.4),
        axis.ticks.x = element_blank(),
        legend.title = element_blank()) +
  scale_y_continuous(labels = scales::comma) +
  geom_segment(data=eQTLpos, aes(x=Chromosome-0.25, 
                                 xend=Chromosome+0.25, 
                                 y=Position-50000, 
                                 yend=Position+50000, 
                                 colour=factor(Type)), 
               size=0.75) +
  scale_colour_manual(labels=c("intergenic","intronic","exonic","ncRNA intronic","multiple SNPs"),
                      values = c(1,2,3,4,5,"red","red"))
#dev.off()
```

Then, we plotted positions and functional annotation of the 13 grouped eQTLs (from the previous plot) on chr4:

i) assigned functional annotations to numbers: 1= intergenic 2=intronic 3=UTR3 4=exonic 5=ncRNA 6=intronic 7=upstream - `chr4.txt`

```{r chr4}
chr4data <- read.table("results/functional_annotation/chr4_sizes_hg19.txt", sep = "\t", header=TRUE)
pos4 <- read.table("results/functional_annotation/chr4.txt", sep = "\t", header=TRUE)
chr4data$Chromosome <- factor(chr4data$Chromosome, levels = chr4data$Chromosome)

#tiff("chr_4_eQTL_SNP_map.tiff", units="in", width=3, height=6, res=300)
ggplot(pos4, aes(Chromosome, Position)) +
  geom_segment(data=pos4, aes(x=Chromosome, 
                              xend=Chromosome, 
                              y=Position-100, 
                              yend=Position+100, 
                              colour=factor(Type)), 
               size=7) +
  scale_colour_manual(labels=c("intergenic","intronic","UTR3","exonic","ncRNA intronic","upstream"), 
                      values = c(1,2,3,4,5,6,7)) +
  theme(panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.background = element_blank(),
      panel.border = element_blank(),
      axis.ticks.y = element_line(),
      legend.title = element_blank(),
      axis.title.x = element_blank(),
      axis.line.y = element_line(),
      axis.text.y = element_text(colour="black", margin = margin(t = 0, r = 5, b = 0, l = 0)),
      axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0))) +
  scale_y_continuous(labels = scales::comma, 
                     breaks = seq(100200000, 100320000, by = 10000),
                     limits = c(100210000,100310000)) +
  scale_x_discrete() +
  guides(colour = guide_legend(override.aes = list(size=1)))
#dev.off()
```

### 3. Percentage of eQTL vs non-eQTL SNPs.

```{r eQTL_vs_noneQTL}
counts <- read.table("results/functional_annotation/AD_eQTLs_vs_non.txt", sep = "\t", header=TRUE)

#tiff("eQTLvsnon.tiff", units="in", width=3, height=6, res=300)
barplot(as.matrix(counts), names.arg=c("eQTL", "non-eQTL"), 
        ylab = "Number of SNPs", ylim = c(0,40), col=c("grey70"))
#dev.off()
```

### 4. Correlation analysis between number of eQTL-eGene interactions and tissue sample size.

Next, we plotted the relationships between the number of significant eQTL-eGene interactions found and the sample size for each tissue: i) for all eQTLs and ii) for just ADH eQTLs

```{r eQTL_vs_sample_size}
#All eQTLs
data <- read.table("data/snps/ADH_eQTL_data.txt", sep = "\t", header=TRUE)

#plot no labels
#tiff("all eQTLs (no labels).tiff", units="in", width=6, height=6, res=300)
ggplot(data, aes(x= sample_size, y=number_of_eQTLs)) + 
  geom_point(size = 2.5, color="#701e7fff") + 
  geom_smooth(method=lm, se=T, linetype="solid", size=0.4, color="#701e7fff") +
  theme_classic() + 
  scale_x_continuous(name="Sample size") + # Number of samples in GTEx tissues
  scale_y_continuous(name="Number of eQTL-eGene interactions") +
  theme(plot.title = element_blank(),
        legend.text = element_text(size=17),
        axis.text=element_text(size=19, color = "black"),
        axis.title=element_text(size=20, color = "black")) +
  stat_cor(method = "pearson", size=8)
#dev.off()

#adding labels for some data points
#tiff("all eQTLs - labelled.tiff", units="in", width=6, height=6, res=300)
ggplot(data, aes(x= sample_size, y=number_of_eQTLs, label=tissue_name)) + 
  geom_point(size = 2.5, color="#701e7fff") + 
  geom_text(data=subset(data, (sample_size > 365 & number_of_eQTLs < 12) | number_of_eQTLs > 16|(number_of_ADH_eQTLs==0 & sample_size==153)),nudge_y=-0.5, size=3.5) +
  geom_smooth(method=lm, se=T, linetype="solid", size=0.4, color="#701e7fff") +
  theme_classic() + 
  scale_x_continuous(name="Sample size", limits=as.numeric(c(NA,530))) +
  scale_y_continuous(name="Number of eQTL-eGene interactions") +
  theme(plot.title = element_blank(),
        legend.text = element_text(size=17),
        axis.text=element_text(size=19, color = "black"),
        axis.title=element_text(size=19, color = "black")) +
  stat_cor(method = "pearson", size=6, label.y = 23)
#dev.off()

#ADH only plot 
#no labels
#tiff("ADH (no labels).tiff", units="in", width=6, height=6, res=300)
ggplot(data, aes(x= sample_size, y=number_of_ADH_eQTLs)) + 
  geom_point(size = 2.5, color="#701e7fff") + 
  geom_smooth(method=lm, se=T, linetype="solid", size=0.4, color="#701e7fff") +
  theme_classic() + 
  scale_x_continuous(name="Sample size") + # Number of samples in GTEx tissues
  scale_y_continuous(name="Number of ADH eQTL-eGene interactions",
                     expand = expand_scale(mult = c(0, .1))) +
  theme(plot.title = element_blank(),
        legend.text = element_text(size=17),
        axis.text=element_text(size=19, color = "black"),
        axis.title=element_text(size=20, color = "black")) +
  stat_cor(method = "pearson", size=8, label.y = 15)
#dev.off()

#with specfic labels
#tiff("ADH - labelled.tiff", units="in", width=6, height=6, res=300)
ggplot(data, aes(x= sample_size, y=number_of_ADH_eQTLs, label=tissue_name)) + 
  geom_point(size = 2.5, color="#701e7fff") + 
  geom_text(data=subset(data, (sample_size > 250 & number_of_ADH_eQTLs < 3) | number_of_ADH_eQTLs > 10 |(number_of_ADH_eQTLs==0 & sample_size==153)),
            nudge_y=-0.4, size=3.5) +
  geom_smooth(method=lm, se=T, linetype="solid", size=0.4, color="#701e7fff") +
  theme_classic() + 
  scale_x_continuous(name="Sample size", limits=as.numeric(c(NA,530)))+
  scale_y_continuous(name="Number of ADH eQTL-eGene interactions") +
  theme(plot.title = element_blank(),
        legend.text = element_text(size=17),
        axis.text=element_text(size=19, color = "black"),
        axis.title=element_text(size=19, color = "black")) +
  stat_cor(method = "pearson", size=6, label.x = 80, label.y = 13)
#dev.off()
```

### 5. Gene Ontology and Pathway Analyses

Gene enrichment and pathway analyses were performed using the g:GOSt module of [g:Profiler](https://biit.cs.ut.ee/gprofiler/). 

### 6. LD analysis

Analysis of linkage disequilibrium (LD) was performed using teh LDmatrix module of the [LDlink](https://ldlink.nci.nih.gov/?tab=home) (version 3.7) tool.

```{r LD_analysis}
get_lower_tri<-function(ld_dt){
  ld_dt[upper.tri(ld_dt)] <- NA
  return(ld_dt)
}

get_upper_tri<-function(ld_dt){
  ld_dt[lower.tri(ld_dt, diag=T)] <- NA
  return(ld_dt)
}

## SNPs on chr3
snps <- c('rs904092', 'rs1789882', 'rs1693457', 'rs9307239', 'rs1789891', 'rs2173201', 'rs1614972',
          'rs2241894', 'rs12639833', 'rs1789924', 'rs1826907')

# R2
# No LD info for rs2066702 --> removed from the analysis
ld_data_r <- read.delim("results/LD/r2_12388_new.txt", header = T)
lower_tri_r <- get_lower_tri(ld_data_r)
snps_sorted <- lower_tri_r$RS_number
lower_ld_melted_r <- melt(lower_tri_r, na.rm = F)

#jpeg("AD_LD_R2.jpeg", units="in", width=8, height=8, res=300)
ggplot(data = lower_ld_melted_r, 
       aes(y=factor(RS_number, levels=lower_tri_r$RS_number), x=variable,
           fill=value)) +
  geom_tile(color='white') +
  scale_fill_gradient2(low='white', high = 'red',space='Lab',
                       name= bquote('LD score '~(R^2)), na.value='transparent') +
  scale_x_discrete("") +
  scale_y_discrete("") + 
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "none") +
  coord_fixed()
#dev.off()

# D'
# No LD info for rs2066702 --> removed from the analysis
ld_data_d <- read.delim("results/LD/d_prime_12388_new.txt", header = T)
upper_tri_d <- get_upper_tri(ld_data_d)
upper_tri_d$RS_number <- snps_sorted
upper_ld_melted_d <- melt(upper_tri_d, na.rm = TRUE)

#jpeg("AD_LD_D.jpeg", units="in", width=8, height=8, res=300)
ggplot(data = upper_ld_melted_d,
       aes(y=factor(RS_number, levels=upper_tri_d$RS_number), x=variable,
           fill=value)) + 
  geom_tile(color='white') + 
  scale_fill_gradient2(low='white', high = '#606060',space='Lab',
                       name="LD score (D')", na.value='transparent') +
  scale_x_discrete("") +
  scale_y_discrete("") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle=90),
        legend.position = "none") +
  coord_fixed()
#dev.off()

ld_data_r <- read.delim("results/LD/d_prime_12388_new.txt", header = T)
lower_tri_r <- get_lower_tri(ld_data_r)
snps_sorted <- lower_tri_r$RS_number
lower_ld_melted_r <- melt(lower_tri_r, na.rm = F)

#jpeg("AD_LD_D.jpeg", units="in", width=8, height=8, res=300)
ggplot(data = lower_ld_melted_r, 
       aes(y=factor(RS_number, levels=lower_tri_r$RS_number), x=variable,
           fill=value)) +
  geom_tile(color='white') +
  scale_fill_gradient2(low='white', high = '#606060',space='Lab',
                       name="LD score (D')", na.value='transparent') +
  scale_x_discrete("") +
  scale_y_discrete("") + 
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "none") +
  coord_fixed()
#dev.off()
```

### 7. Tissue-specific *ADH* gene expression and the effect sizes of SNPs impacting *ADH* genes.

```{r ADH1A, fig.width=10, fig.height=8}
# ADH gene regulation
snps <- c('rs2066702','rs904092', 'rs1789882', 'rs1693457', 'rs9307239', 'rs1789891',
          'rs2173201', 'rs1614972', 'rs2241894', 'rs12639833', 'rs1789924', 'rs1826907')
tissues <- c("Adipose - Subcutaneous", "Adipose - Visceral (Omentum)", "Adrenal Gland",
             "Artery - Aorta", "Artery - Tibial", "Brain - Caudate (basal ganglia)",
             "Breast - Mammary Tissue", "Cells - Transformed fibroblasts", 
             "Colon - Sigmoid", "Colon - Transverse",
             "Esophagus - Gastroesophageal Junction", "Esophagus - Mucosa",
             "Esophagus - Muscularis", "Heart - Atrial Appendage",
             "Heart - Left Ventricle", "Lung", "Muscle - Skeletal", "Nerve - Tibial",
             "Pancreas", "Skin - Not Sun Exposed (Suprapubic)", "Skin - Sun Exposed (Lower leg)",
             "Small Intestine - Terminal Ileum", "Spleen", "Thyroid")

adh.df <- read.delim('data/eGenes/ADH_significant_eQTLs.txt', sep='\t', header=TRUE)
adh.sub <- subset(adh.df, select=c(SNP, Gene_Name, Tissue, Effect_Size))
adh1a <- subset(adh.sub, Gene_Name=='ADH1A')
adh1b <- subset(adh.sub, Gene_Name=='ADH1B')
adh1c <- subset(adh.sub, Gene_Name=='ADH1C')
adh4 <- subset(adh.sub, Gene_Name=='ADH4')

# Plotting ADH1A eQTLs
adh1a_melt <- melt(adh1a, id=c("SNP", "Tissue"), measure.vars = c("Effect_Size"))
#unique(adh1a_melt$SNP)
adh1a_dcast <- dcast(adh1a_melt, Tissue ~ SNP)

i <- c()
for(s in 1:length(snps)){
  if(!(snps[s] %in% unique(adh1a_melt$SNP))){
    i <- c(i , snps[s])
    #adh1a_dcast$i <- NA
  }
}

adh1a_dcast$rs9307239 <- NA; adh1a_dcast$rs2173201 <- NA; adh1a_dcast$rs1614972 <- NA;
adh1a_dcast$rs2241894 <- NA; adh1a_dcast$rs12639833 <- NA; adh1a_dcast$rs1789924 <- NA
adh1a_dcast <- remove_rownames(adh1a_dcast) %>% column_to_rownames(var="Tissue")
adh1a_mat <- as.matrix(adh1a_dcast)

colors = colorRamp2(c(-1,0,1), c('#d7191c','#ffffbf','#1a9641'), space = "RGB")

#pdf("figures/ADH1A_gene_regulation.pdf", width = 16, height = 16)
adh1a_ht <- Heatmap(adh1a_mat, name = "Effect size",
                    col = colors, na_col = "white",
                    cluster_rows = FALSE, cluster_columns = FALSE, row_names_side = "left",
                    column_names_side = "top", column_order = snps,
                    width = unit(10, "cm"), height = unit(15, "cm"),
                    border = TRUE, rect_gp = gpar(col = "grey"),
                    heatmap_legend_param = list(legend_direction = "horizontal",
                                                legend_width = unit(5, "cm")))
draw(adh1a_ht, heatmap_legend_side = "bottom")
#dev.off()
```

```{r ADH1B, fig.width=8, fig.height=4}
# Plotting ADH1B eQTLs
adh1b_melt <- melt(adh1b, id=c("SNP", "Tissue"), measure.vars = c("Effect_Size"))
adh1b_dcast <- dcast(adh1b_melt, Tissue ~ SNP)
i <- c()
for(s in 1:length(snps)){
  if(!(snps[s] %in% unique(adh1b_melt$SNP))){
    i <- c(i , snps[s])
    #adh1a_dcast$i <- NA
  }
}

adh1b_dcast$rs2066702 <- NA; adh1b_dcast$rs904092 <- NA; adh1b_dcast$rs1789882 <- NA;
adh1b_dcast$rs1693457 <- NA; adh1b_dcast$rs9307239 <- NA; adh1b_dcast$rs1789891 <- NA
adh1b_dcast$rs2173201 <- NA; adh1b_dcast$rs1614972 <- NA; adh1b_dcast$rs2241894 <- NA
adh1b_dcast$rs12639833 <- NA; adh1b_dcast$rs1789924 <- NA
adh1b_dcast <- remove_rownames(adh1b_dcast) %>% column_to_rownames(var="Tissue")
adh1b_mat <- as.matrix(adh1b_dcast)

#pdf("figures/ADH1B_gene_regulation.pdf", width = 16, height = 16)
adh1b_ht <- Heatmap(adh1b_mat, name = "Effect size",
                    col = colors, na_col = "white",
                    cluster_rows = FALSE, cluster_columns = FALSE, row_names_side = "left",
                    column_names_side = "top", column_order = snps,
                    width = unit(10, "cm"), height = unit(6, "cm"),
                    border = TRUE, rect_gp = gpar(col = "grey"),
                    heatmap_legend_param = list(legend_direction = "horizontal",
                                                legend_width = unit(5, "cm")))
draw(adh1b_ht, heatmap_legend_side = "bottom")
#dev.off()
```

```{r ADH1C, fig.width=10, fig.height=7}
# Plotting ADH1C eQTLs
adh1c_melt <- melt(adh1c, id=c("SNP", "Tissue"), measure.vars = c("Effect_Size"))
adh1c_dcast <- dcast(adh1c_melt, Tissue ~ SNP)
i <- c()
for(s in 1:length(snps)){
  if(!(snps[s] %in% unique(adh1c_melt$SNP))){
    i <- c(i , snps[s])
    #adh1a_dcast$i <- NA
  }
}

adh1c_dcast$rs2066702 <- NA; adh1c_dcast$rs904092 <- NA; adh1c_dcast$rs1789882 <- NA;
adh1c_dcast$rs1693457 <- NA
adh1c_dcast <- remove_rownames(adh1c_dcast) %>% column_to_rownames(var="Tissue")
adh1c_mat <- as.matrix(adh1c_dcast)

#pdf("figures/ADH1C_gene_regulation.pdf", width = 16, height = 16)
adh1c_ht <- Heatmap(adh1c_mat, name = "Effect size",
                    col = colors, na_col = "white",
                    cluster_rows = FALSE, cluster_columns = FALSE, row_names_side = "left",
                    column_names_side = "top", column_order = snps,
                    width = unit(10, "cm"), height = unit(10, "cm"),
                    border = TRUE, rect_gp = gpar(col = "grey"),
                    heatmap_legend_param = list(legend_direction = "horizontal",
                                                legend_width = unit(5, "cm")))
draw(adh1c_ht, heatmap_legend_side = "bottom")
#dev.off()
```

```{r ADH4, fig.width=8, fig.height=4}
# Plotting ADH4 eQTLs
adh4_melt <- melt(adh4, id=c("SNP", "Tissue"), measure.vars = c("Effect_Size"))
adh4_dcast <- dcast(adh4_melt, Tissue ~ SNP)
i <- c()
for(s in 1:length(snps)){
  if(!(snps[s] %in% unique(adh4_melt$SNP))){
    i <- c(i , snps[s])
    #adh1a_dcast$i <- NA
  }
}

adh4_dcast$rs2066702 <- NA; adh4_dcast$rs904092 <- NA; adh4_dcast$rs1789882 <- NA;
adh4_dcast$rs1693457 <- NA; adh4_dcast$rs1789891 <- NA; adh4_dcast$rs1826907 <- NA
adh4_dcast <- remove_rownames(adh4_dcast) %>% column_to_rownames(var="Tissue")
adh4_mat <- as.matrix(adh4_dcast)

#pdf("figures/ADH4_gene_regulation.pdf", width = 16, height = 16)
adh4_ht <- Heatmap(adh4_mat, name = "Effect size",
                   col = colors, na_col = "white",
                   cluster_rows = FALSE, cluster_columns = FALSE, row_names_side = "left",
                   column_names_side = "top", column_order = snps,
                   width = unit(10, "cm"), height = unit(2, "cm"),
                   border = TRUE, rect_gp = gpar(col = "grey"),
                   heatmap_legend_param = list(legend_direction = "horizontal",
                                               legend_width = unit(5, "cm")))
draw(adh4_ht, heatmap_legend_side = "bottom")
#dev.off()
```

```{r ADH_expression, fig.width=10, fig.height=7}
# ADH gene expression
expr.df <- read.delim('data/eGenes/ADH_gene_expression.txt', sep='\t', header=TRUE)
expr.df <- na.omit(expr.df)

expr_melt <- melt(expr.df, id=c("Gene_Name", "Tissue"), measure.vars = c("log.TPM.1."))
expr_dcast <- dcast(expr_melt, Tissue ~ Gene_Name)
expr_dcast <- remove_rownames(expr_dcast) %>% column_to_rownames(var="Tissue")
expr_mat <- as.matrix(expr_dcast)

colors = colorRamp2(c(0, max(expr_mat)), c("white","black"), space = "RGB")

#pdf("figures/adh_aldh_gene_expression_heatmap.pdf", width = 16, height = 16)
ht_ge <- Heatmap(expr_mat, name = "log10(TPM+1)",
                 col = colors, na_col = "white",
                 cluster_rows = FALSE, cluster_columns = FALSE, row_names_side = "left",
                 column_names_side = "top", column_names_gp = gpar(fontface = "italic"),
                 width = unit(6, "cm"), height = unit(13, "cm"),
                 border = TRUE, rect_gp = gpar(col = "grey"),
                 heatmap_legend_param = list(legend_direction = "horizontal",
                                             legend_width = unit(5, "cm")))
draw(ht_ge, heatmap_legend_side = "bottom")
#dev.off()
```

### 8. Regulatory annotation of *ADH* SNPs

Promoter-enhancer histone marks and affected regulatory motifs were obtained from the [HaploReg](https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php) website (version 4.1).