covid19_hg.Rmd

---
title: "Reviewing summary statistics from the COVID19-HGI"
author: "[Sander W. van der Laan, PhD](https://swvanderlaan.github.io) | @swvanderlaan | s.w.vanderlaan@gmail.com"
date: "`r Sys.Date()`"
output:
  html_notebook: 
    cache: yes
    code_folding: hide
    collapse: yes
    df_print: paged
    fig.align: center
    fig_caption: yes
    fig_width: 12
    fig_height: 8
    fig_retina: 2
    highlight: default #tango kate haddock
    theme: spacelab #cosmo lumen paper yeti
    toc: yes
    toc_float:
      collapsed: no
      smooth_scroll: yes
mainfont: Arial
subtitle: ""
editor_options:
  chunk_output_type: inline
---

```{r global_options, include = FALSE}
# further define some knitr-options.
knitr::opts_chunk$set(fig.width = 12, fig.height = 8, fig.path = 'Figures/', 
                      warning = TRUE, # show warnings during codebook generation
                      message = TRUE, # show messages during codebook generation
                      error = TRUE, # do not interrupt codebook generation in case of errors; usually better for debugging
                      echo = TRUE,  # show R code
                      eval = TRUE)
ggplot2::theme_set(ggplot2::theme_minimal())
pander::panderOptions("table.split.table", Inf)

```

# Setup
We will clean the environment, setup the locations, define colors, and create a datestamp.

_Clean the environment._
```{r echo = FALSE}
rm(list = ls())
```

_Set locations and working directories..._
```{r LocalSystem, echo = FALSE}
### Operating System Version
### Mac Pro
# ROOT_loc = "/Volumes/EliteProQx2Media"
# GENOMIC_loc = "/Users/svanderlaan/iCloud/Genomics"
# AEDB_loc = paste0(GENOMIC_loc, "/AE-AAA_GS_DBs")
# LAB_loc = paste0(GENOMIC_loc, "/LabBusiness")
# PLINK_loc=paste0(ROOT_loc,"/PLINK")
# AEGSQC_loc =  paste0(PLINK_loc, "/_AE_ORIGINALS/AEGS_COMBINED_QC2018")
# PROJECT_loc = paste0(ROOT_loc, "/PLINK/analyses/baseline")

### MacBook
ROOT_loc = "/Users/swvanderlaan"
# GENOMIC_loc = paste0(ROOT_loc, "/iCloud/Genomics")
# LAB_loc = paste0(GENOMIC_loc, "/LabBusiness")
PLINK_loc = paste0(ROOT_loc,"/PLINK")
COVID19_loc = paste0(PLINK_loc,"/analyses/covid19_hgi")
PROJECT_loc = paste0(COVID19_loc,"/covid19_hg")
SCRIPTS_loc = paste0(PROJECT_loc,"/SCRIPTS")
# RAWDATA_loc = paste0(COVID19_loc,"/RESULTS/COVID19_HGI_ANA5_20200429") # first release
RAWDATA_loc = paste0(COVID19_loc,"/RESULTS/COVID19_HGI_ANA5_20200508") # second release

### SOME VARIABLES WE NEED DOWN THE LINE
# TRAIT_OF_INTEREST = "UCORBIO" # Phenotype
PROJECTNAME = "COVID19" # primary goal of this analysis
SUBPROJECTNAME = "PRELIMINARYRELEASE2" # e.g. study name

### FILTERING
# we set a threshold on the allele frequency
AF_THRESHOLD = 0.005 

### DEBUGGING MODE
# we set the debugging mode, or real mode, options: YES/NO
DEBUGGING_FLAG="YES"

cat("\nCreate a new analysis directory...\n")
ifelse(!dir.exists(file.path(PROJECT_loc, "/",paste0(PROJECTNAME,".",SUBPROJECTNAME))), 
       dir.create(file.path(PROJECT_loc, "/",paste0(PROJECTNAME,".",SUBPROJECTNAME))), 
       FALSE)
ANALYSIS_loc = paste0(PROJECT_loc,"/",paste0(PROJECTNAME,".",SUBPROJECTNAME))

ifelse(!dir.exists(file.path(ANALYSIS_loc, "/PLOTS")), 
       dir.create(file.path(ANALYSIS_loc, "/PLOTS")), 
       FALSE)
PLOT_loc = paste0(ANALYSIS_loc,"/PLOTS")

ifelse(!dir.exists(file.path(PLOT_loc, "/QC")), 
       dir.create(file.path(PLOT_loc, "/QC")), 
       FALSE)
QC_loc = paste0(PLOT_loc,"/QC")

ifelse(!dir.exists(file.path(ANALYSIS_loc, "/OUTPUT")), 
       dir.create(file.path(ANALYSIS_loc, "/OUTPUT")), 
       FALSE)
OUT_loc = paste0(ANALYSIS_loc, "/OUTPUT")

ifelse(!dir.exists(file.path(ANALYSIS_loc, "/BASELINE")), 
       dir.create(file.path(ANALYSIS_loc, "/BASELINE")), 
       FALSE)
BASELINE_loc = paste0(ANALYSIS_loc, "/BASELINE")

setwd(paste0(PROJECT_loc))
getwd()
list.files()

```

_... a package-installation function ..._
```{r Function: installations}
install.packages.auto <- function(x) { 
  x <- as.character(substitute(x)) 
  if (isTRUE(x %in% .packages(all.available = TRUE))) { 
    eval(parse(text = sprintf("require(\"%s\")", x)))
  } else { 
    # Update installed packages - this may mean a full upgrade of R, which in turn
    # may not be warrented. 
    #update.packages(ask = FALSE) 
    eval(parse(text = sprintf("install.packages(\"%s\", dependencies = TRUE, repos = \"http://cran-mirror.cs.uu.nl/\")", x)))
  }
  if (isTRUE(x %in% .packages(all.available = TRUE))) { 
    eval(parse(text = sprintf("require(\"%s\")", x)))
  } else {
    source("http://bioconductor.org/biocLite.R")
    # Update installed packages - this may mean a full upgrade of R, which in turn
    # may not be warrented.
    #biocLite(character(), ask = FALSE) 
    eval(parse(text = sprintf("biocLite(\"%s\")", x)))
    eval(parse(text = sprintf("require(\"%s\")", x)))
  }
}
```

_... and load those packages._
```{r loading_packages}
install.packages.auto("readr")
install.packages.auto("optparse")
install.packages.auto("tools")
install.packages.auto("dplyr")
install.packages.auto("tidyr")
install.packages.auto("naniar")

# To get 'data.table' with 'fwrite' to be able to directly write gzipped-files
# Ref: https://stackoverflow.com/questions/42788401/is-possible-to-use-fwrite-from-data-table-with-gzfile
# install.packages("data.table", repos = "https://Rdatatable.gitlab.io/data.table")
library(data.table)

install.packages.auto("tidyverse")
install.packages.auto("knitr")
install.packages.auto("DT")
install.packages.auto("eeptools")

install.packages.auto("haven")
install.packages.auto("tableone")

install.packages.auto("ggpubr")

install.packages.auto("BlandAltmanLeh")

install.packages.auto("ggplot2")

install.packages.auto("CMplot")

# For regional association plots
devtools::install_github("jrs95/gassocplot")
```


_We will create a datestamp and define the Utrecht Science Park Colour Scheme_.
```{r Setting: Colors}

Today = format(as.Date(as.POSIXlt(Sys.time())), "%Y%m%d")
Today.Report = format(as.Date(as.POSIXlt(Sys.time())), "%A, %B %d, %Y")

### UtrechtScienceParkColoursScheme
###
### WebsitetoconvertHEXtoRGB:http://hex.colorrrs.com.
### Forsomefunctionsyoushoulddividethesenumbersby255.
###
###	No.	Color			      HEX	(RGB)						              CHR		  MAF/INFO
###---------------------------------------------------------------------------------------
###	1	  yellow			    #FBB820 (251,184,32)				      =>	1		or 1.0>INFO
###	2	  gold			      #F59D10 (245,157,16)				      =>	2		
###	3	  salmon			    #E55738 (229,87,56)				      =>	3		or 0.05<MAF<0.2 or 0.4<INFO<0.6
###	4	  darkpink		    #DB003F ((219,0,63)				      =>	4		
###	5	  lightpink		    #E35493 (227,84,147)				      =>	5		or 0.8<INFO<1.0
###	6	  pink			      #D5267B (213,38,123)				      =>	6		
###	7	  hardpink		    #CC0071 (204,0,113)				      =>	7		
###	8	  lightpurple	    #A8448A (168,68,138)				      =>	8		
###	9	  purple			    #9A3480 (154,52,128)				      =>	9		
###	10	lavendel		    #8D5B9A (141,91,154)				      =>	10		
###	11	bluepurple		  #705296 (112,82,150)				      =>	11		
###	12	purpleblue		  #686AA9 (104,106,169)			      =>	12		
###	13	lightpurpleblue	#6173AD (97,115,173/101,120,180)	=>	13		
###	14	seablue			    #4C81BF (76,129,191)				      =>	14		
###	15	skyblue			    #2F8BC9 (47,139,201)				      =>	15		
###	16	azurblue		    #1290D9 (18,144,217)				      =>	16		or 0.01<MAF<0.05 or 0.2<INFO<0.4
###	17	lightazurblue	  #1396D8 (19,150,216)				      =>	17		
###	18	greenblue		    #15A6C1 (21,166,193)				      =>	18		
###	19	seaweedgreen	  #5EB17F (94,177,127)				      =>	19		
###	20	yellowgreen		  #86B833 (134,184,51)				      =>	20		
###	21	lightmossgreen	#C5D220 (197,210,32)				      =>	21		
###	22	mossgreen		    #9FC228 (159,194,40)				      =>	22		or MAF>0.20 or 0.6<INFO<0.8
###	23	lightgreen	  	#78B113 (120,177,19)				      =>	23/X
###	24	green			      #49A01D (73,160,29)				      =>	24/Y
###	25	grey			      #595A5C (89,90,92)				        =>	25/XY	or MAF<0.01 or 0.0<INFO<0.2
###	26	lightgrey		    #A2A3A4	(162,163,164)			      =>	26/MT
###
###	ADDITIONAL COLORS
###	27	midgrey			#D7D8D7
###	28	verylightgrey	#ECECEC"
###	29	white			#FFFFFF
###	30	black			#000000
###----------------------------------------------------------------------------------------------

uithof_color = c("#FBB820","#F59D10","#E55738","#DB003F","#E35493","#D5267B",
                 "#CC0071","#A8448A","#9A3480","#8D5B9A","#705296","#686AA9",
                 "#6173AD","#4C81BF","#2F8BC9","#1290D9","#1396D8","#15A6C1",
                 "#5EB17F","#86B833","#C5D220","#9FC228","#78B113","#49A01D",
                 "#595A5C","#A2A3A4", "#D7D8D7", "#ECECEC", "#FFFFFF", "#000000")

uithof_color_legend = c("#FBB820", "#F59D10", "#E55738", "#DB003F", "#E35493",
                        "#D5267B", "#CC0071", "#A8448A", "#9A3480", "#8D5B9A",
                        "#705296", "#686AA9", "#6173AD", "#4C81BF", "#2F8BC9",
                        "#1290D9", "#1396D8", "#15A6C1", "#5EB17F", "#86B833",
                        "#C5D220", "#9FC228", "#78B113", "#49A01D", "#595A5C",
                        "#A2A3A4", "#D7D8D7", "#ECECEC", "#FFFFFF", "#000000")
### ----------------------------------------------------------------------------
```


# COVID19 Host Genetic Initiative

Some friends and colleagues asked about the preliminary results from the genetics studies on COVID19. One large-scale effort is the [COVID-19 Host Genetics Initiative (COVID19-HGI)](https://www.covid19hg.org){target="_blank"}[^1] which aims to bring together the human genetics community to generate, share and analyze data to learn the genetic determinants of COVID-19 susceptibility, severity and outcomes. Ultimately, these discoveries could help to generate hypotheses for drug repurposing, identify individuals at unusually high or low risk, and contribute to global knowledge of the biology of SARS-CoV-2 infection and disease.

In this repository I just wanted to share a quick look at the (preliminary) results from the rounds of meta-analysis of genome-wide association studies (GWAS) included. 

> First Release Notes: An important caveat here is that these data really are just that: preliminary. The phenotype definition used here does not accurately capture disease susceptibility and severity. Read more about this [here](https://www.covid19hg.org/blog/2020-04-30-first-round-of-meta-analysis-results/){target="_blank"}.

> Second Release Notes: Meta-analysis was done with inverse variance weighting. Positions are on GRCh37. An allele frequency (AF) filter of 0.0001 and an imputation quality (INFO) filter of 0.6 was applied to each study prior to meta-analysis. Read more [here](https://www.covid19hg.org/results/){target="_blank"}.


## Studies

The participating studies:

- *UK Biobank (UKBB)* (analyzed by Danny Wilson, Konrad Karczewski, Tomoko Nakanishi and Manuel Rivas)
- *FinnGen (FinnGen)* (analyzed by Juha Karjalainen)
- *Genes and Health (GNH)* (analyzed by David Van Heel)
- *BioMe biobank (BioMe_RGN)* (analyzed by Arden Moscati)
- *Helix Exome+ COVID-19 Phenotypes (Helix)* 
- *LifeLines Global Screening Array (LifelinesCyto)* (analyzed by Patrick Deelen)
- *LifeLines CytoSNP (LifelinesGsa)* (analyzed by Patrick Deelen)
- *Netherlands Twin Register (NTR)*
- *Partners Healthcare Biobank (PHBB)*


## Data description

Data column descriptions:

- #CHR: chromosome, 
- POS: chromosome position in build 37, 
- REF: non-effect allele, 
- ALT: effect allele (beta is for this allele), 
- SNP: #CHR:POS:REF:ALT, 
- {STUDY}_AF_Allele2: allele frequency in {STUDY}, 
- all_meta_N: number of studies that had the variant after AF and INFO filtering and as such were used for the meta, 
- all_inv_var_meta_beta: effect size on log(OR) scale, 
- all_inv_var_meta_sebeta: standard error of effect size, 
- all_inv_var_meta_p: p-value, 
- all_inv_var_het_p: p-value from Cochran's Q heterogeneity test.


## Analyses

There are several analyses carried out for COVID-19 susceptibility, severity and outcomes. 

- ANA2, Hospitalized vs. non-hospitalized
- ANA5, susceptibility (affected vs. population)
- ANA7, COVID-19 predicted by flu-like symptoms


### Analysis 2: 20200508-results-ANA2_inv_var_meta

*Phenotype*: ANA2, Hospitalized vs. non-hospitalized. *Population*: All.

*Total cases*: 716. *Total controls*: 616.

*Contributing studies (n_cases, n_controls)*: 

- FinnGen_FIN (30, 179), 
- GNH_SAS (34, 13), 
- LifelinesGsa_EUR (6, 183), 
- PHBB_EUR (17, 57), 
- UKBB_EUR (629, 184).

Downloads:

- COVID19_HGI_ANA2_20200513.txt.gz
- COVID19_HGI_ANA2_20200513.txt.gz.tbi


### Analysis 5: 20200508-results-ANA5_ALL_inv_var_meta

*Phenotype*: ANA5, susceptibility (affected vs. population). *Population*: All.

*Total cases*: 1,678. *Total controls*: 674,635.

*Contributing studies (n_cases, n_controls)*: 

- BioMe_RGN_EUR (20, 10169), 
- FinnGen_FIN (209, 203431), 
- GNH_SAS (64, 27351), 
- LifelinesCyto_EUR (62, 5750), 
- LifelinesGsa_EUR (189, 15975), 
- NTR_CEU (163, 3160), 
- PHBB_AFR (34, 2169), 
- PHBB_EUR (74, 28950), 
- UKBB_AFR (51, 6585), 
- UKBB_EUR (812, 371095).

Downloads: 

- COVID19_HGI_ANA5_20200513.txt.gz
- COVID19_HGI_ANA5_20200513.txt.gz.tbi


### Analysis 7: 20200508-results-ANA7_inv_var_meta

*Phenotype*: ANA7, COVID-19 predicted by flu-like symptoms. *Population*: European.

*Total cases*: 1,294. *Total controls*: 26,969.

*Contributing studies (n_cases, n_controls)*: 

- LifelinesGsa_EUR (443, 15406), 
- LifelinesCyto_EUR (194, 5451), 
- Helix_EUR (657, 6112).

Downloads:

- COVID19_HGI_ANA7_20200513.txt.gz
- COVID19_HGI_ANA7_20200513.txt.gz.tbi


[^1]: [The COVID-19 Host Genetics Initiative _The COVID-19 Host Genetics Initiative, a global initiative to elucidate the role of host genetic factors in susceptibility and severity of the SARS-CoV-2 virus pandemic_ **EJHG 2020**.](https://www.nature.com/articles/s41431-020-0636-6).

# Loading data

We have downloaded the data from the [site](https://www.covid19hg.org/){target="_blank"} and stored this locally.

```{r LoadData}
cat("Getting results ...\n")

cat("\n* analysis 2 ...\n")
prelim.gwas.raw <- fread(paste0(RAWDATA_loc,"/COVID19_HGI_ANA2_20200513.txt.gz"),
                         header = TRUE,
                         verbose = FALSE,
                         showProgress = TRUE)

cat("\n* analysis 5 ...\n")
prelim.gwas.raw5 <- fread(paste0(RAWDATA_loc,"/COVID19_HGI_ANA5_20200513.txt.gz"),
                         header = TRUE,
                         verbose = FALSE,
                         showProgress = TRUE)

cat("\n* analysis 7 ...\n")
prelim.gwas.raw7 <- fread(paste0(RAWDATA_loc,"/COVID19_HGI_ANA7_20200513.txt.gz"),
                         header = TRUE,
                         verbose = FALSE,
                         showProgress = TRUE)

```


Here we take a random sample of the data for debugging this script if `DEBUGGING_FLAG` equals `YES`.
```{r RandomSample}
# for debug
# random sampling

if (DEBUGGING_FLAG == "YES") {
  cat("We are in debugging mode.")
  
  sample.size = 500000
  
  prelim.gwas.raw <- prelim.gwas.raw %>%
    sample_n(., sample.size)
  
  prelim.gwas.raw5 <- prelim.gwas.raw5 %>%
    sample_n(., sample.size)
  
  prelim.gwas.raw7 <- prelim.gwas.raw7 %>%
    sample_n(., sample.size)


} else if (DEBUGGING_FLAG == "NO") {
  cat("Real-mode: using all data.")
  
  } else {
    cat("This message must be attributable to human error. Please check the debugging flag.")
    
    }

```

Summarizing the data, and get a `head` of each dataset.
```{r SummaryData ANA2}
cat("\n* analysis 2 ...\n")
str(prelim.gwas.raw)
dim(prelim.gwas.raw)

# summary(prelim.gwas.raw)
head(prelim.gwas.raw)

```

```{r SummaryData ANA5}
cat("\n* analysis 5 ...\n")
str(prelim.gwas.raw5)
dim(prelim.gwas.raw5)

# summary(prelim.gwas.raw5)
head(prelim.gwas.raw5)

```

```{r SummaryData ANA7}
cat("\n* analysis 7 ...\n")
str(prelim.gwas.raw7)
dim(prelim.gwas.raw7)

# summary(prelim.gwas.raw7)
head(prelim.gwas.raw7)

```

# Calculate mean AF

Here we will calculate the mean allele frequency across cohorts to plot these and to create stratified QQ-plots.
```{r Prep histogram AF}
prelim.gwas.raw <- prelim.gwas.raw %>%
  mutate(meanAF = rowMeans(select(., ends_with("_AF_Allele2")),
                          na.rm = TRUE))

prelim.gwas.raw5 <- prelim.gwas.raw5 %>%
  mutate(meanAF = rowMeans(select(., ends_with("_AF_Allele2")),
                          na.rm = TRUE))

prelim.gwas.raw7 <- prelim.gwas.raw7 %>%
  mutate(meanAF = rowMeans(select(., ends_with("_AF_Allele2")),
                          na.rm = TRUE))
```

# Visualisation 

Here we plot the raw results, as is, so no filtering on allele frequency etc. First we make a regular QQ plot, than we stratify it on allele frequency and imputation quality score.

## N studies
  
```{r NGWAS}

p1 <- ggpubr::gghistogram(prelim.gwas.raw,
                    x = "all_meta_N",
                    y = "..count..",
                    title = "COVID-19 hospitalisation",
                    xlab = "n studies",
                    ylab = "n variants",
                    color = uithof_color[3],
                    fill = uithof_color[3],
                    add = "mean", add.params = list(color = uithof_color[3], linetype = 2))

p2 <- ggpubr::gghistogram(prelim.gwas.raw,
                    x = "all_meta_N",
                    y = "..count..",
                    title = "COVID-19 susceptibility",
                    xlab = "n studies",
                    ylab = "n variants",
                    color = uithof_color[16],
                    fill = uithof_color[16],
                    add = "mean", add.params = list(color = uithof_color[3], linetype = 2))

p3 <- ggpubr::gghistogram(prelim.gwas.raw,
                    x = "all_meta_N",
                    y = "..count..",
                    title = "COVID-19 predicted by flu-like symptoms",
                    xlab = "n studies",
                    ylab = "n variants",
                    color = uithof_color[23],
                    fill = uithof_color[23],
                    add = "mean", add.params = list(color = uithof_color[3], linetype = 2))

require(patchwork)

p1 / p2 / p3

rm(p1, p2, p3)

```


## QQ-plots

A simple QQ-plot can nicely show deviation from the expected distribution; this could indidate a lack of power, subtle effects originating from confounders (technical artefacts from the experiment, or population substructure), or true polygenecity[^2].

[^2]: McCarthy, M. I. _et al._ Genome-wide association studies for complex traits: consensus, uncertainty and challenges. **Nat. Rev. Genet.** 9, 356–369 (2008).

_Prepping the data for plotting._
```{r PrepQQMH}
Pmap <- prelim.gwas.raw %>%
  select(., SNP, `#CHR`, POS, all_inv_var_meta_p)

Pmap5 <- prelim.gwas.raw5 %>%
  select(., SNP, `#CHR`, POS, all_inv_var_meta_p)

Pmap7 <- prelim.gwas.raw7 %>%
  select(., SNP, `#CHR`, POS, all_inv_var_meta_p)

```

Plotting.
```{r QQPlotRaw}

CMplot(Pmap, plot.type = "q",
       col = uithof_color[3], cex.lab = 1,25,
       ylim = 10, 
       file.output = FALSE,
       main = "COVID-19 hospitalisation",
       verbose = TRUE)

CMplot(Pmap5, plot.type = "q",
       col = uithof_color[16], cex.lab = 1,25,
       ylim = 10, 
       file.output = FALSE,
       main = "COVID-19 susceptibility",
       verbose = TRUE)

CMplot(Pmap7, plot.type = "q",
       col = uithof_color[23], cex.lab = 1,25,
       ylim = 10, 
       file.output = FALSE,
       main = "COVID-19 predicted by flu-like symptoms",
       verbose = TRUE)

```


## Manhattan plot

The Manhattan visualizes all the results.

```{r MHPlotRaw}

CMplot(Pmap, plot.type = "m",
       col = uithof_color, cex.lab = 1.25,
       # col = c(uithof_color[3], uithof_color[26]), cex.lab = 1.25,
       file.output = FALSE,
       main = "COVID-19 hospitalisation",
       ylim = 10,
       threshold = 5e-8, threshold.col = uithof_color[26], threshold.lty = 2, threshold.lwd = 1,
       amplify = TRUE,
       chr.labels = c(1:22,"X"),
       verbose = TRUE)

CMplot(Pmap5, plot.type = "m",
       col = uithof_color, cex.lab = 1.25,
       # col = c(uithof_color[16], uithof_color[26]), cex.lab = 1.25,
       file.output = FALSE,
       main = "COVID-19 susceptibility",
       ylim = 10,
       threshold = 5e-8, threshold.col = uithof_color[26], threshold.lty = 2, threshold.lwd = 1,
       amplify = TRUE,
       chr.labels = c(1:22,"X"),
       verbose = TRUE)

CMplot(Pmap7, plot.type = "m",
       col = uithof_color, cex.lab = 1.25,
       # col = c(uithof_color[23], uithof_color[26]), cex.lab = 1.25,
       file.output = FALSE,
       main = "COVID-19 predicted by flu-like symptoms",
       ylim = 10,
       threshold = 5e-8, threshold.col = uithof_color[26], threshold.lty = 2, threshold.lwd = 1,
       amplify = TRUE,
       chr.labels = c(1:22), #,"X", "Y"),
       verbose = TRUE)

```


## Allele frequencies

We will plot the distribution of the allele frequencies. 

```{r histogram AF}
p1 <- ggpubr::gghistogram(prelim.gwas.raw,
                    x = "meanAF", y = "..count..",
                    title = "COVID-19 hospitalisation",
                    xlab = "",
                    ylab = "",
                    color = uithof_color[3],
                    fill = uithof_color[3],
                    add = "mean", add.params = list(color = uithof_color[26], linetype = 2, size = 1.25)) #+
  # annotate("text", label = "mean AF", x = 0.30, y = 5e5, colour = uithof_color[26]) #+
  # geom_vline(xintercept = 0.01, colour = uithof_color[3], lty = 1) + #annotate("text", label = "AF = 1%", x = 0.1, y = 5e5, colour = uithof_color[3]) +
  # geom_vline(xintercept = 0.02, colour = uithof_color[13], lty = 1) + #annotate("text", label = "AF = 2%", x = 0.1, y = 4.5e5, colour = uithof_color[13]) +
  # geom_vline(xintercept = 0.03, colour = uithof_color[16], lty = 1) #+ annotate("text", label = "AF = 3%", x = 0.1, y = 4e5, colour = uithof_color[16])

p2 <- ggpubr::gghistogram(prelim.gwas.raw5,
                    x = "meanAF", y = "..count..",
                    title = "COVID-19 susceptibility",
                    xlab = "",
                    ylab = "n variants",
                    color = uithof_color[16],
                    fill = uithof_color[16],
                    add = "mean", add.params = list(color = uithof_color[26], linetype = 2, size = 1.25)) #+
  # annotate("text", label = "mean AF", x = 0.30, y = 5e5, colour = uithof_color[26]) 

p3 <- ggpubr::gghistogram(prelim.gwas.raw7,
                    x = "meanAF", y = "..count..",
                    title = "COVID-19 predicted by flu-like symptoms",
                    xlab = "allele frequency\n(mean across studies)",
                    ylab = "",
                    color = uithof_color[23],
                    fill = uithof_color[23],
                    add = "mean", add.params = list(color = uithof_color[26], linetype = 2, size = 1.25)) #+
  # annotate("text", label = "mean AF", x = 0.30, y = 5e5, colour = uithof_color[26]) 


require(patchwork)

p1 / p2 / p3

rm(p1, p2, p3)

```

## Stratified QQ-plots

Here we will add in QQ-plots stratified on the mean allele frequency, this will show the lambda and aids in interpreting the inflation per AF-stratum.


# Filtering data

Given that sample sizes (of the cases) remain small, we will apply a filter at mean AF = `r AF_THRESHOLD`.
```{r filtering data}

prelim.gwas.raw.qc <- prelim.gwas.raw %>%
  filter(meanAF > AF_THRESHOLD)

prelim.gwas.raw5.qc <- prelim.gwas.raw5 %>%
  filter(meanAF > AF_THRESHOLD)

prelim.gwas.raw7.qc <- prelim.gwas.raw7 %>%
  filter(meanAF > AF_THRESHOLD)

rm(prelim.gwas.raw, prelim.gwas.raw5, prelim.gwas.raw7)

```

## QQ-plot after filter

```{r prepping QQPlotQC}
PmapQC <- prelim.gwas.raw.qc %>%
  select(., SNP, `#CHR`, POS, all_inv_var_meta_p)

PmapQC5 <- prelim.gwas.raw5.qc %>%
  select(., SNP, `#CHR`, POS, all_inv_var_meta_p)

PmapQC7 <- prelim.gwas.raw7.qc %>%
  select(., SNP, `#CHR`, POS, all_inv_var_meta_p)


```

```{r QQPlotQC}
CMplot(PmapQC, plot.type = "q",
       col = uithof_color[3], cex.lab = 1,25,
       ylim = 10, 
       file.output = FALSE,
       main = "COVID-19 hospitalisation",
       verbose = TRUE)

CMplot(PmapQC5, plot.type = "q",
       col = uithof_color[16], cex.lab = 1,25,
       ylim = 10, 
       file.output = FALSE,
       main = "COVID-19 susceptibility",
       verbose = TRUE)

CMplot(PmapQC7, plot.type = "q",
       col = uithof_color[23], cex.lab = 1,25,
       ylim = 10, 
       file.output = FALSE,
       main = "COVID-19 predicted by flu-like symptoms",
       verbose = TRUE)

```


## Manhattan plot after filter

```{r MHPlotQC}

CMplot(PmapQC, plot.type = "m",
       col = uithof_color, cex.lab = 1.25,
       # col = c(uithof_color[3], uithof_color[26]), cex.lab = 1.25,
       file.output = FALSE,
       main = "COVID-19 hospitalisation",
       ylim = 10,
       threshold = 5e-8, threshold.col = uithof_color[26], threshold.lty = 2, threshold.lwd = 1,
       amplify = TRUE,
       chr.labels = c(1:22,"X"),
       verbose = TRUE)

CMplot(PmapQC5, plot.type = "m",
       col = uithof_color, cex.lab = 1.25,
       # col = c(uithof_color[16], uithof_color[26]), cex.lab = 1.25,
       file.output = FALSE,
       main = "COVID-19 susceptibility",
       ylim = 10,
       threshold = 5e-8, threshold.col = uithof_color[26], threshold.lty = 2, threshold.lwd = 1,
       amplify = TRUE,
       chr.labels = c(1:22,"X"),
       verbose = TRUE)

CMplot(PmapQC7, plot.type = "m",
       col = uithof_color, cex.lab = 1.25,
       # col = c(uithof_color[23], uithof_color[26]), cex.lab = 1.25,
       file.output = FALSE,
       main = "COVID-19 predicted by flu-like symptoms",
       ylim = 10,
       threshold = 5e-8, threshold.col = uithof_color[26], threshold.lty = 2, threshold.lwd = 1,
       amplify = TRUE,
       chr.labels = c(1:22), #,"X", "Y"),
       verbose = TRUE)

```

# Version & change log

------

    Version:      v1.0.1
    Last update:  2020-05-16
    Written by:   Sander W. van der Laan (s.w.vanderlaan-2[at]umcutrecht.nl).
    Description:  Script to review the (preliminary) meta-GWAS results from the COVID19-HGI.
    Minimum requirements: R version 3.4.3 (2017-06-30) -- 'Single Candle', Mac OS X El Capitan
    
    
    Changes log
    * v1.0.1 Second round meta-analysis, preliminary results.
    * v1.0.0 Initial version. First round meta-analysis, preliminary results. 

------


# Session information

```{r eval = TRUE}
sessionInfo()
```

# Saving environment
Actually: we don't want to save the data, this makes the `.RData`-file too large...
```{r Saving}

rm(prelim.gwas.raw, prelim.gwas.raw5, prelim.gwas.raw7,
   prelim.gwas.raw.qc, prelim.gwas.raw5.qc, prelim.gwas.raw7.qc,
   Pmap, Pmap5, Pmap7,
   PmapQC, PmapQC5, PmapQC7)
save.image(paste0(PROJECT_loc, "/",Today,".",PROJECTNAME,".",SUBPROJECTNAME,".RData"))

```

------
<sup>&copy; 1979-2020 Sander W. van der Laan | s.w.vanderlaan-2[at]umcutrecht.nl | [swvanderlaan.github.io](https://swvanderlaan.github.io){target="_blank"}.</sup>
------


*Footnotes*