\documentclass[9pt,a4paper,]{extarticle}

\usepackage{f1000_styles}

\usepackage[pdfborder={0 0 0}]{hyperref}

\usepackage[super]{natbib}
\bibliographystyle{unsrtnat}


%% maxwidth is the original width if it is less than linewidth
%% otherwise use linewidth (to make sure the graphics do not exceed the margin)
\makeatletter
\def\maxwidth{ %
  \ifdim\Gin@nat@width>\linewidth
    \linewidth
  \else
    \Gin@nat@width
  \fi
}
\makeatother

\usepackage{color}
\usepackage{fancyvrb}
\newcommand{\VerbBar}{|}
\newcommand{\VERB}{\Verb[commandchars=\\\{\}]}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
% Add ',fontsize=\small' for more characters per line
\usepackage{framed}
\definecolor{shadecolor}{RGB}{248,248,248}
\newenvironment{Shaded}{\begin{snugshade}}{\end{snugshade}}
\newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{#1}}
\newcommand{\DecValTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\FloatTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\CharTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\StringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\ImportTok}[1]{#1}
\newcommand{\CommentTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\OtherTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{#1}}
\newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\VariableTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.81,0.36,0.00}{\textbf{#1}}}
\newcommand{\BuiltInTok}[1]{#1}
\newcommand{\ExtensionTok}[1]{#1}
\newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.77,0.63,0.00}{#1}}
\newcommand{\RegionMarkerTok}[1]{#1}
\newcommand{\InformationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\WarningTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AlertTok}[1]{\textcolor[rgb]{0.94,0.16,0.16}{#1}}
\newcommand{\ErrorTok}[1]{\textcolor[rgb]{0.64,0.00,0.00}{\textbf{#1}}}
\newcommand{\NormalTok}[1]{#1}

% disable code chunks background
%\renewenvironment{Shaded}{}{}

% disable section numbers
\setcounter{secnumdepth}{0}

\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}



\usepackage{amsthm}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\theoremstyle{definition}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}
\theoremstyle{definition}
\newtheorem{example}{Example}
\theoremstyle{definition}
\newtheorem{exercise}{Exercise}
\theoremstyle{remark}
\newtheorem*{remark}{Remark}
\newtheorem*{solution}{Solution}
\begin{document}
\pagestyle{front}

\title{Using regulatory genomics data to interpret the function of disease variants and prioritise genes from expression studies}

\author{Enrico Ferrero}

\maketitle
\thispagestyle{front}

\begin{abstract}
The identification of therapeutic targets is a critical step in the research and developement of new drugs, with several drug discovery programmes failing because of a weak linkage between target and disease. Genome-wide association studies and large-scale gene expression experiments are providing insights into the biology of several common diseases, but the complexity of transcriptional regulation mechanisms often limits our understanding of how genetic variation can influence changes in gene expression. Several initiatives in the field of regulatory genomics are aiming to close this gap by systematically identifying and cataloguing regulatory elements such as promoters and enhacers across different tissues and cell types. In this Bioconductor workflow, we will explore how different types of regulatory genomic data can be used for the functional interpretation of disease-associated variants and for the prioritisation of gene lists from gene expression experiments.
\end{abstract}

\section*{Keywords}
bioconductor; r; rstats; regulatory genomics; functional genomics; genetics; gwas; transcriptomics; integration; multiomics


\clearpage
\pagestyle{main}

\section{Introduction}\label{introduction}

Discovering and bringing new drugs to the market is a long, expensive and inefficient process \citep{Waring2015, DiMasi2016}.
The majority of drug discovery programmes fail for efficacy reasons \citep{Harrison2016}, with up to 40\% of these failures due to lack of a clear link between the target and the disease under investigation \citep{Cook2014}.
Target selection, the first step in drug discovery programmes, is thus a critical decision point.
It has previously been shown that therapeutic targets with a genetic link to the disease under investigation are more likely to progress through the drug discovery pipeline, suggesting that genetics can be used as a tool to prioritise and validate drug targets in early discovery \citep{Plenge2013, Nelson2015}.

One of the biggest challenges in translating findings from genome-wide association studies (GWASs) to therapies is that the great majority of single nucleotide polymorphisms (SNPs) associated with disease are found in non-coding regions of the genome, and therefore cannot be easily linked to a target gene \citep{Maurano2012}.
Many of these SNPs could be regulatory variants, affecting the expression of nearby or distal genes by interfering with the transcriptional process \citep{Ward2012}.

The most established way to map disease-associated regulatory variants to target genes is to use expression quantitative trait loci (eQTLs) \citep{Albert2015}, variants that affect the expression of specific genes.
The GTEx consortium profiled eQTLs across 44 human tissues by performing a large-scale mapping of genome-wide correlations between genetic variants and gene expression \citep{GTEx2017}.
However, depending on the power of the study, it might not be possible to detect all existing regulatory variants as eQTLs.
An alternative is to use information on the location of promoters and distal enhancers across the genome and link these regulatory elements to their target genes.
Large, multi-centre initiatives such as ENCODE \citep{ENCODE2012}, Roadmap Epigenomics \citep{Roadmap2015} and BLUEPRINT \citep{Adams2012, Stunnenberg2016} mapped regulatory elements in the genome by profiling a number of chromatin features, including DNase hypersensitive sites (DHSs), several types of histone marks and binding of chromatin-associated proteins in a large number of cells and tissues.
Similarly, the FANTOM consortium used cap analysis of gene expression (CAGE) to identify promoters and enhancers across hundreds of cells and tissues \citep{Fantom2014}.

Knowing that a certain stretch of DNA is an enhancer is however not informative of the target gene(s).
One way to infer links between enhancers and promoters \emph{in silico} is to identify significant correlations across a large panel of cell types, an approach that was used for distal and promoter DHSs \citep{Thurman2012} as well as for CAGE-defined promoters and enhancers \citep{Andersson2014}.
Experimental methods to assay interactions between regulatory elements also exist.
Chromatin interaction analysis by paired-end tag sequencing (ChIA-PET) \citep{Fullwood2009, Zhang2013} couples chromatin immunoprecipitation with DNA ligation to identify DNA regions interacting thanks to the binding of a specific protein.
Promoter capture Hi-C \citep{Mifsud2015, Javierre2016} extends chromatin conformation capture by using ``baits'' to enrich for promoter interactions and increase resolution.

Overall, linking genetic variants to their candidate target genes is not straightforward, not only because of the complexity of the human genome and transcriptional regulation, but also because of the variety of data types and approaches that can be used.
To address this problem, we developed STOPGAP, a database of disease variants mapped to their most likely target gene(s) using several different types of regulatory genomic data \citep{Shen2017}.
The database is currently undergoing a major overhaul and will eventually be superseded by \href{https://github.com/Ensembl/postgap}{POSTGAP}.
A valid and recent alternative is INFERNO \citep{Amlie-Wolf2017}, though it does only rely on eQTL data for target gene assignment.
These resources implement some or all of the approaches that will be reviewed in the workflow and constitute good entry points for identifying the most likely target gene(s) of regulatory SNPs. However, as they tend to hide much of the complexity involved in the process, we will not use them and rely on the original datasets instead.

In this workflow, we will explore how regulatory genomic data can be used to connect the genetic and transcriptional layers by providing a framework for the discovery of novel therapeutic targets.
We will use eQTL data from GTEx \citep{GTEx2017}, FANTOM5 correlations between promoters and enhancers \citep{Andersson2014} and promoter capture Hi-C data \citep{Javierre2016} to annotate significant GWAS variants to putative target genes and to prioritise genes obtained from a differential expression analysis (Figure 1).

\hypertarget{htmlwidget-6c56f25cf019d83b2621}{}

\label{fig:figure1}Diagram showing a schematic representation of the workflow and the steps involved.

\section{Workflow}\label{workflow}

\subsection{Install required packages}\label{install-required-packages}

R version 3.4.2 and Bioconductor version 3.6 were used for the analysis.
The code below will install all required packages and dependencies from Bioconductor and CRAN:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{source}\NormalTok{(}\StringTok{"https://bioconductor.org/biocLite.R"}\NormalTok{)}
\CommentTok{# uncomment the following line to install packages}
\CommentTok{#biocLite(c("clusterProfiler", "DESeq2", "GenomicFeatures", "GenomicInteractions", "GenomicRanges", "ggplot2", "Gviz", "gwascat", "InteractionSet", "recount", "pheatmap", "RColorBrewer", "rtracklayer", "R.utils", "splitstackshape", "VariantAnnotation"))}
\end{Highlighting}
\end{Shaded}

\subsection{Gene expression data and differential gene expression analysis}\label{gene-expression-data-and-differential-gene-expression-analysis}

We start with a common scenario: we ran a RNA-seq experiment comparing patients with a disease and healthy individuals, and would like to discover key disease genes and potential therapeutic targets by integrating genetic information in our analysis.

The RNA-seq data we will be using comes from blood of patients with systemic lupus erythematosus (SLE) and healthy controls \citep{Hung2015}.
SLE is a chronic autoimmune disorder that can affect several organs with a significant unmet medical need \citep{Kaul2016}.
It is a complex and remarkably heterogeneous disease, in terms of both genetics and clinical manifestations \citep{Marion2014}.
Early diagnosis and classification of SLE remain extremely challenging \citep{Amezcua-Guerra2015}.

In the original study \citep{Hung2015}, the authors explore transcripts bound by Ro60, an RNA-binding protein against which some SLE patients produce autoantibodies.
They identify Alu retroelements among these transcripts and use RNA-seq data to check their expression levels, observing that Alu elements are significantly more expressed in SLE patients, and particularly in those patients with anti-Ro antibodies and with a higher interferon signature metric (ISM).

We are going to use \texttt{recount} \citep{Collado-Torres2017} to obtain gene-level counts:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(recount)}
\CommentTok{# uncomment the following line to download dataset}
\CommentTok{#download_study("SRP062966")}
\KeywordTok{load}\NormalTok{(}\KeywordTok{file.path}\NormalTok{(}\StringTok{"SRP062966"}\NormalTok{, }\StringTok{"rse_gene.Rdata"}\NormalTok{))}
\NormalTok{rse <-}\StringTok{ }\KeywordTok{scale_counts}\NormalTok{(rse_gene)}
\end{Highlighting}
\end{Shaded}

Other Bioconductor packages that can be used to access data from gene expression experiments directly in R are \texttt{GEOquery} \citep{Davis2007} and \texttt{ArrayExpress} \citep{Kauffmann2009}.

We have 117 samples overall.
This is what the matrix of counts looks like:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{assay}\NormalTok{(rse)[}\DecValTok{1}\OperatorTok{:}\DecValTok{3}\NormalTok{, }\DecValTok{1}\OperatorTok{:}\DecValTok{3}\NormalTok{]}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                    SRR2443263 SRR2443262 SRR2443261
## ENSG00000000003.14         19          6         10
## ENSG00000000005.5           0          0          0
## ENSG00000000419.12        489        238        224
\end{verbatim}

Each gene is a row and each sample is a column.
We note that genes are annotated using the GENCODE \citep{Harrow2012} v25 annotation, which will be useful later on.

To check how we can split samples between cases and controls, we can have a look at the metadata contained in the \texttt{characteristics} column, which is a \texttt{CharacterList} object:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{head}\NormalTok{(rse}\OperatorTok{$}\NormalTok{characteristics, }\DecValTok{3}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## CharacterList of length 3
## [[1]] disease status: healthy tissue: whole blood anti-ro: control ism: control
## [[2]] disease status: healthy tissue: whole blood anti-ro: control ism: control
## [[3]] disease status: healthy tissue: whole blood anti-ro: control ism: control
\end{verbatim}

We have information about the disease status of the sample, the tissue of origin, the presence and level of anti-ro autoantibodies and the value of the ISM.
However, we note that basic information such as age or gender is missing.

We can create some new columns with the available information so that they can be used for downstream analyses.
We will also make sure that they are encoded as factors and that the correct reference layer is used:

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# disease status}
\NormalTok{rse}\OperatorTok{$}\NormalTok{disease_status <-}\StringTok{ }\KeywordTok{sapply}\NormalTok{(rse}\OperatorTok{$}\NormalTok{characteristics, }\StringTok{"["}\NormalTok{, }\DecValTok{1}\NormalTok{)}
\NormalTok{rse}\OperatorTok{$}\NormalTok{disease_status <-}\StringTok{ }\KeywordTok{sub}\NormalTok{(}\StringTok{"disease status: "}\NormalTok{, }\StringTok{""}\NormalTok{, rse}\OperatorTok{$}\NormalTok{disease_status)}
\NormalTok{rse}\OperatorTok{$}\NormalTok{disease_status <-}\StringTok{ }\KeywordTok{sub}\NormalTok{(}\StringTok{"systemic lupus erythematosus }\CharTok{\textbackslash{}\textbackslash{}}\StringTok{(SLE}\CharTok{\textbackslash{}\textbackslash{}}\StringTok{)"}\NormalTok{, }\StringTok{"SLE"}\NormalTok{, rse}\OperatorTok{$}\NormalTok{disease_status)}
\NormalTok{rse}\OperatorTok{$}\NormalTok{disease_status <-}\StringTok{ }\KeywordTok{factor}\NormalTok{(rse}\OperatorTok{$}\NormalTok{disease_status, }\DataTypeTok{levels =} \KeywordTok{c}\NormalTok{(}\StringTok{"healthy"}\NormalTok{, }\StringTok{"SLE"}\NormalTok{))}
\CommentTok{# tissue}
\NormalTok{rse}\OperatorTok{$}\NormalTok{tissue <-}\StringTok{ }\KeywordTok{sapply}\NormalTok{(rse}\OperatorTok{$}\NormalTok{characteristics, }\StringTok{"["}\NormalTok{, }\DecValTok{2}\NormalTok{)}
\NormalTok{rse}\OperatorTok{$}\NormalTok{tissue <-}\StringTok{ }\KeywordTok{sub}\NormalTok{(}\StringTok{"tissue: "}\NormalTok{, }\StringTok{""}\NormalTok{, rse}\OperatorTok{$}\NormalTok{tissue)}
\NormalTok{rse}\OperatorTok{$}\NormalTok{tissue <-}\StringTok{ }\KeywordTok{factor}\NormalTok{(rse}\OperatorTok{$}\NormalTok{tissue)}
\CommentTok{# anti-ro}
\NormalTok{rse}\OperatorTok{$}\NormalTok{anti_ro <-}\StringTok{ }\KeywordTok{sapply}\NormalTok{(rse}\OperatorTok{$}\NormalTok{characteristics, }\StringTok{"["}\NormalTok{, }\DecValTok{3}\NormalTok{)}
\NormalTok{rse}\OperatorTok{$}\NormalTok{anti_ro <-}\StringTok{ }\KeywordTok{sub}\NormalTok{(}\StringTok{"anti-ro: "}\NormalTok{, }\StringTok{""}\NormalTok{, rse}\OperatorTok{$}\NormalTok{anti_ro)}
\NormalTok{rse}\OperatorTok{$}\NormalTok{anti_ro <-}\StringTok{ }\KeywordTok{factor}\NormalTok{(rse}\OperatorTok{$}\NormalTok{anti_ro)}
\CommentTok{# ism}
\NormalTok{rse}\OperatorTok{$}\NormalTok{ism <-}\StringTok{ }\KeywordTok{sapply}\NormalTok{(rse}\OperatorTok{$}\NormalTok{characteristics, }\StringTok{"["}\NormalTok{, }\DecValTok{4}\NormalTok{)}
\NormalTok{rse}\OperatorTok{$}\NormalTok{ism <-}\KeywordTok{sub}\NormalTok{(}\StringTok{"ism: "}\NormalTok{, }\StringTok{""}\NormalTok{, rse}\OperatorTok{$}\NormalTok{ism)}
\NormalTok{rse}\OperatorTok{$}\NormalTok{ism <-}\StringTok{ }\KeywordTok{factor}\NormalTok{(rse}\OperatorTok{$}\NormalTok{ism)}
\end{Highlighting}
\end{Shaded}

We can check how many samples we have in each group (note that we ignore \texttt{tissue} as it's always \texttt{whole\ blood}):

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{metadata <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{disease_status =}\NormalTok{ rse}\OperatorTok{$}\NormalTok{disease_status, }\DataTypeTok{anti_ro.ism =} \KeywordTok{paste}\NormalTok{(rse}\OperatorTok{$}\NormalTok{anti_ro, rse}\OperatorTok{$}\NormalTok{ism, }\DataTypeTok{sep =} \StringTok{"."}\NormalTok{))}
\KeywordTok{table}\NormalTok{(metadata)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##               anti_ro.ism
## disease_status control.control high.ISM_high high.ISM_low med.ISM_high
##        healthy              18             0            0            0
##        SLE                   0            23            1           21
##               anti_ro.ism
## disease_status med.ISM_low none.ISM_high none.ISM_low
##        healthy           0             0            0
##        SLE               2            31           21
\end{verbatim}

Now we are ready to perform a simple differential gene expression analysis with \texttt{DESeq2} \citep{Love2014}.
Note that we remove genes with a low number of counts (less than 50 across all 117 samples) to speed up execution and reduce the memory footprint:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(DESeq2)}
\NormalTok{dds <-}\StringTok{ }\KeywordTok{DESeqDataSet}\NormalTok{(rse, }\OperatorTok{~}\StringTok{ }\NormalTok{disease_status)}
\NormalTok{dds <-}\StringTok{ }\KeywordTok{DESeq}\NormalTok{(dds)}
\NormalTok{dds <-}\StringTok{ }\NormalTok{dds[}\KeywordTok{rowSums}\NormalTok{(}\KeywordTok{counts}\NormalTok{(dds)) }\OperatorTok{>=}\StringTok{ }\DecValTok{50}\NormalTok{, ]}
\end{Highlighting}
\end{Shaded}

We used an extremely simple model; in the real world we should be accounting for co-variables, potential confounders and interactions between them.
For example, age and gender are usually included in this type of analysis, but we don't have access to this information for this dataset.
Similarly, the value of the ISM and the presence of anti-Ro autoantibodies can't be included in the analysis due to the fact that these variables are collinear with the disease status variable (i.e.: the value of both \texttt{anti\_ro} and \texttt{ism} is \texttt{control} for all samples with \texttt{disease\_status} equal to \texttt{healthy}.)
Like \texttt{DESeq2}, \texttt{edgeR} \citep{Robinson2010} and \texttt{limma} \citep{Ritchie2015} can also deal with multiple cofactors and different experimental designs, and constitute good alternatives for performing differential expression analyses.

We can now look at the data in more detail to assess if we can observe a separation between the SLE and healthy samples and whether any batch effect is visible.\\
We use the variance stabilising transformation (VST) \citep{Anders2010} for visualisation purposes:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{vsd <-}\StringTok{ }\KeywordTok{vst}\NormalTok{(dds, }\DataTypeTok{blind =} \OtherTok{FALSE}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

We will use the \href{http://cran.r-project.org/package=pheatmap}{\texttt{pheatmap}} and \href{https://cran.r-project.org/package=RColorBrewer}{\texttt{RColorBrewer}} packages to perform hierarchical clustering of the samples (Figure 2):

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(pheatmap)}
\KeywordTok{library}\NormalTok{(RColorBrewer)}
\NormalTok{sampleDists <-}\StringTok{ }\KeywordTok{dist}\NormalTok{(}\KeywordTok{t}\NormalTok{(}\KeywordTok{assay}\NormalTok{(vsd)))}
\NormalTok{sampleDistMatrix <-}\StringTok{ }\KeywordTok{as.matrix}\NormalTok{(sampleDists)}
\NormalTok{annotation =}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\KeywordTok{colData}\NormalTok{(vsd)[}\KeywordTok{c}\NormalTok{(}\StringTok{"anti_ro"}\NormalTok{, }\StringTok{"ism"}\NormalTok{, }\StringTok{"disease_status"}\NormalTok{)], }\DataTypeTok{row.names =} \KeywordTok{rownames}\NormalTok{(sampleDistMatrix))}
\NormalTok{colors <-}\StringTok{ }\KeywordTok{colorRampPalette}\NormalTok{(}\KeywordTok{rev}\NormalTok{(}\KeywordTok{brewer.pal}\NormalTok{(}\DecValTok{9}\NormalTok{, }\StringTok{"Blues"}\NormalTok{)))(}\DecValTok{255}\NormalTok{)}
\KeywordTok{pheatmap}\NormalTok{(sampleDistMatrix, }\DataTypeTok{clustering_distance_rows =}\NormalTok{ sampleDists, }\DataTypeTok{clustering_distance_cols =}\NormalTok{ sampleDists, }\DataTypeTok{clustering_method =} \StringTok{"complete"}\NormalTok{, }\DataTypeTok{annotation_col =}\NormalTok{ annotation, }\DataTypeTok{col =}\NormalTok{ colors, }\DataTypeTok{show_rownames =} \OtherTok{FALSE}\NormalTok{, }\DataTypeTok{show_colnames =} \OtherTok{FALSE}\NormalTok{, }\DataTypeTok{cellwidth =} \DecValTok{2}\NormalTok{, }\DataTypeTok{cellheight =} \DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{figure}

{\centering \includegraphics{biocondutor-regulatory-genomics-workflow_files/figure-latex/figure2-1} 

}

\caption{Heatmap showing Euclidean distances between samples clustered using complete linkage. Disease status and other experimental factors are visualised as column annotations.}\label{fig:figure2}
\end{figure}

While there isn't an unambiguous split between healthy and disease samples, the most distinct clusters (bottom right and top left) are entirely composed of SLE samples, with the central cluster containing all healthy samples and a number of SLE ones.
The clusters don't appear to be due to the ISM or the presence of anti-Ro autoantibodies.

Similarly, we can perform a principal component analysis (PCA) on the most variable 500 genes (Figure 3).
Note that we load \texttt{ggplot2} \citep{Wickham2009} to modify the look of the plot:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(ggplot2)}
\KeywordTok{plotPCA}\NormalTok{(vsd, }\DataTypeTok{intgroup =} \StringTok{"disease_status"}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{coord_fixed}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\begin{figure}

{\centering \includegraphics{biocondutor-regulatory-genomics-workflow_files/figure-latex/figure3-1} 

}

\caption{Scatter plot showing results of a PCA with samples coloured according to their disease status.}\label{fig:figure3}
\end{figure}

We can see some separation of healthy and SLE samples along both PC1 and PC2, though some SLE samples appear very similar to the healthy ones.
No obvious batch effects are visible from this plot.

Next, we select genes that are differentially expressed below a 0.05 adjusted \emph{p}-value threshold:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{res <-}\StringTok{ }\KeywordTok{results}\NormalTok{(dds, }\DataTypeTok{alpha =} \FloatTok{0.05}\NormalTok{)}
\KeywordTok{summary}\NormalTok{(res)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## out of 32820 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)     : 4829, 15% 
## LFC < 0 (down)   : 2709, 8.3% 
## outliers [1]     : 0, 0% 
## low counts [2]   : 2548, 7.8% 
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
\end{verbatim}

We can visualise the shrunken log2 fold changes using an MA plot (Figure 4):

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{res_lfc <-}\StringTok{ }\KeywordTok{lfcShrink}\NormalTok{(dds, }\DataTypeTok{coef =} \DecValTok{2}\NormalTok{)}
\KeywordTok{plotMA}\NormalTok{(res_lfc, }\DataTypeTok{ylim =} \KeywordTok{c}\NormalTok{(}\OperatorTok{-}\DecValTok{3}\NormalTok{, }\DecValTok{3}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\begin{figure}

{\centering \includegraphics{biocondutor-regulatory-genomics-workflow_files/figure-latex/figure4-1} 

}

\caption{MA plot showing genes differentially expressed (red dots) in SLE patients compared to healthy patients.}\label{fig:figure4}
\end{figure}

We observe large numbers of genes differentially expressed in both directions and across a range of fold changes, though the majority of significant genes appear to be upregulated in disease.

For convenience, we will save our differentially expressed genes (DEGs) in another object and map the GENCODE gene IDs to gene symbols using the annotation in the original \texttt{RangedSummarizedExperiment} object

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{degs <-}\StringTok{ }\KeywordTok{subset}\NormalTok{(res, padj }\OperatorTok{<}\StringTok{ }\FloatTok{0.05}\NormalTok{)}
\NormalTok{degs <-}\StringTok{ }\KeywordTok{merge}\NormalTok{(}\KeywordTok{rowData}\NormalTok{(rse), }\KeywordTok{as.data.frame}\NormalTok{(degs), }\DataTypeTok{by.x =} \StringTok{"gene_id"}\NormalTok{, }\DataTypeTok{by.y =} \StringTok{"row.names"}\NormalTok{, }\DataTypeTok{all =} \OtherTok{FALSE}\NormalTok{)}
\KeywordTok{head}\NormalTok{(degs, }\DecValTok{3}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## DataFrame with 3 rows and 9 columns
##                            gene_id bp_length symbol   baseMean
##                        <character> <integer> <list>  <numeric>
## ENSG00000000003 ENSG00000000003.14      4535 TSPAN6   8.739822
## ENSG00000000419 ENSG00000000419.12      1207   DPM1 431.485085
## ENSG00000000457 ENSG00000000457.13      6883  SCYL3 686.579323
##                 log2FoldChange      lfcSE      stat       pvalue
##                      <numeric>  <numeric> <numeric>    <numeric>
## ENSG00000000003     -0.4750382 0.18374822 -2.585267 9.730366e-03
## ENSG00000000419      0.5559772 0.10117967  5.494950 3.908216e-08
## ENSG00000000457      0.1927081 0.05928191  3.250707 1.151185e-03
##                         padj
##                    <numeric>
## ENSG00000000003 4.161281e-02
## ENSG00000000419 2.182977e-06
## ENSG00000000457 7.922475e-03
\end{verbatim}

\subsection{Accessing GWAS data}\label{accessing-gwas-data}

The differential expression analysis resulted in several thousands of DEGs.
Since we know that genes with high levels of differential expression are more likely to harbour disease-associated variants \citep{Chen2008} and that therapeutic targets with genetic evidence are more likely to progress through the drug discovery pipeline \citep{Nelson2015}, one way to prioritise them is to check which of these can be genetically linked to SLE.
To get hold of relevant GWAS data, we will be using the \texttt{gwascat} Bioconductor package \citep{Carey2017a}, which provides an interface to the GWAS catalog \citep{MacArthur2017}.
An alternative is to use the GRASP \citep{Eicher2015} database with the \texttt{grasp2db} \citep{Carey2017b} package.

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(gwascat)}
\CommentTok{# uncomment the following line to download file and build the gwasloc object all in one step}
\CommentTok{#snps <- makeCurrentGwascat()}
\CommentTok{# uncomment the following line to download file}
\CommentTok{#download.file("http://www.ebi.ac.uk/gwas/api/search/downloads/alternative", destfile = "gwas_catalog_v1.0.1-associations_e90_r2017-12-04.tsv")}
\NormalTok{snps <-}\StringTok{ }\KeywordTok{read.delim}\NormalTok{(}\StringTok{"gwas_catalog_v1.0.1-associations_e90_r2017-12-04.tsv"}\NormalTok{, }\DataTypeTok{check.names =} \OtherTok{FALSE}\NormalTok{, }\DataTypeTok{stringsAsFactors =} \OtherTok{FALSE}\NormalTok{)}
\NormalTok{snps <-}\StringTok{ }\NormalTok{gwascat}\OperatorTok{:::}\KeywordTok{gwdf2GRanges}\NormalTok{(snps, }\DataTypeTok{extractDate =} \StringTok{"2017-12-04"}\NormalTok{)}
\KeywordTok{genome}\NormalTok{(snps) <-}\StringTok{ "GRCh38"}
\KeywordTok{head}\NormalTok{(snps, }\DecValTok{3}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## gwasloc instance with 3 records and 37 attributes per record.
## Extracted:  2017-12-04 
## Genome:  GRCh38 
## Excerpt:
## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames                 ranges strand | DISEASE/TRAIT        SNPS
##          <Rle>              <IRanges>  <Rle> |   <character> <character>
##   [1]     chr1 [203186754, 203186754]      * | YKL-40 levels   rs4950928
##   [2]    chr13 [ 39776775,  39776775]      * |     Psoriasis   rs7993214
##   [3]    chr15 [ 78513681,  78513681]      * |   Lung cancer   rs8034191
##         P-VALUE
##       <numeric>
##   [1]     1e-13
##   [2]     2e-06
##   [3]     3e-18
##   -------
##   seqinfo: 23 sequences from GRCh38 genome; no seqlengths
\end{verbatim}

\texttt{snps} is a \texttt{gwasloc} object which is simply a wrapper around a \texttt{GRanges} object, the standard way to represent genomic ranges in Bioconductor.

We note here that the GWAS catalog uses GRCh38 coordinates, the same assembly used in the GENCODE v25 annotation.
When integrating genomic datasets from different sources it is essential to ensure that the same genome assembly is used, especially because many datasets in the public domain are still using GRCh37 coordinates.
As we will see below, it is possible and relatively straightforward to convert genomic coordinates between genome assemblies.

We can select only SNPs that are associated with SLE:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{snps <-}\StringTok{ }\KeywordTok{subsetByTraits}\NormalTok{(snps, }\DataTypeTok{tr =} \StringTok{"Systemic lupus erythematosus"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

We can visualise these as a Manhattan plot to look at the distribution of GWAS \emph{p}-values over chromosomes on a negative log\textsubscript{10} scale (Figure 5):
Note that \emph{p}-values lower than 1 x 10\textsuperscript{-25} are truncated in the figure:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{traitsManh}\NormalTok{(}\DataTypeTok{gwr =}\NormalTok{ snps, }\DataTypeTok{sel =}\NormalTok{ snps, }\DataTypeTok{traits =} \StringTok{"Systemic lupus erythematosus"}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{xlab}\NormalTok{(}\StringTok{"SLE GWAS SNPs"}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{ylab}\NormalTok{(}\StringTok{"-log10(p-value)"}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{theme}\NormalTok{(}\DataTypeTok{legend.position =} \StringTok{"none"}\NormalTok{,}
        \DataTypeTok{strip.text.x =} \KeywordTok{element_text}\NormalTok{(}\DataTypeTok{size =} \DecValTok{6}\NormalTok{),}
        \DataTypeTok{axis.text.x =} \KeywordTok{element_blank}\NormalTok{(),}
        \DataTypeTok{axis.ticks.x =} \KeywordTok{element_blank}\NormalTok{())}
\end{Highlighting}
\end{Shaded}

\begin{figure}

{\centering \includegraphics{biocondutor-regulatory-genomics-workflow_files/figure-latex/figure5-1} 

}

\caption{Manhattan plot showing GWAS variants significantly associated with SLE.}\label{fig:figure5}
\end{figure}

We observe several hits across most chromosomes, with many of them below a genome-wide significant threshold (\emph{p}-value \textless{} 1 x 10\textsuperscript{-8}), suggesting that genetics plays an important role in the pathogenesis of SLE.

We note here that genotyping arrays typically include a very small fraction of all possible SNPs in the human genome, and there is no guarantee that the \emph{tag} SNPs on the array are the true casual SNPs \citep{Bush2012}.
The alleles of other SNPs can be imputed from tag SNPs thanks to the structure of linkage disequilibrium (LD) blocks present in chromosomes.
Thus, when linking variants to target genes in a real-world setting, it is important to take into consideration neighbouring SNPs that are in high LD (e.g.: r\textsuperscript{2} \textgreater{} 0.8) and inherited with the tag SNPs.
Unfortunately, at the time of writing there is no straightforward way to perform this LD expansion step using R or Bioconductor packages, possibly because of the large amount of reference data required.
The \texttt{ldblock} package \citep{Carey2017c} used to provide this functionality by downloading the HapMap data from the NCBI website, but the dataset was retired in 2016.
At present, the best option to do this programmatically is probably to query the Ensembl REST API \citep{Yates2015}.

\subsection{Annotation of coding and proximal SNPs to target genes}\label{annotation-of-coding-and-proximal-snps-to-target-genes}

In order to annotate these variants, we need a a \texttt{TxDb} object, a reference of where transcripts are located on the genome.
We can build this using the \texttt{GenomicFeatures} \citep{Lawrence2013} package and the GENCODE v25 gene annotation:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(GenomicFeatures)}
\CommentTok{# uncomment the following line to download file}
\CommentTok{#download.file("ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.annotation.gff3.gz", destfile = "gencode.v25.annotation.gff3.gz")}
\NormalTok{txdb <-}\StringTok{ }\KeywordTok{makeTxDbFromGFF}\NormalTok{(}\StringTok{"gencode.v25.annotation.gff3.gz"}\NormalTok{)}
\NormalTok{txdb <-}\StringTok{ }\KeywordTok{keepStandardChromosomes}\NormalTok{(txdb)}
\end{Highlighting}
\end{Shaded}

We also have to convert the \texttt{gwasloc} object into a standard \texttt{GRanges} object:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{snps <-}\StringTok{ }\KeywordTok{GRanges}\NormalTok{(snps)}
\end{Highlighting}
\end{Shaded}

Let's check if the \texttt{gwasloc} and \texttt{TxDb} object use the same notation for chromosomes:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{seqlevelsStyle}\NormalTok{(snps)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] "UCSC"
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{seqlevelsStyle}\NormalTok{(txdb)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] "UCSC"
\end{verbatim}

OK, they do.
Now we can annotate our SNPs to genes using the \texttt{VariantAnnotation} \citep{Obenchain2014} package:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(VariantAnnotation)}
\NormalTok{snps_anno <-}\StringTok{ }\KeywordTok{locateVariants}\NormalTok{(snps, txdb, }\KeywordTok{AllVariants}\NormalTok{())}
\NormalTok{snps_anno <-}\StringTok{ }\KeywordTok{unique}\NormalTok{(snps_anno)}
\end{Highlighting}
\end{Shaded}

We use the \texttt{QUERYID} column in \texttt{snps\_anno} to recover metadata such as SNP IDs and GWAS \emph{p}-values from the original \texttt{snps} object:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{snps_metadata <-}\StringTok{ }\NormalTok{snps[snps_anno}\OperatorTok{$}\NormalTok{QUERYID]}
\KeywordTok{mcols}\NormalTok{(snps_anno) <-}\StringTok{ }\KeywordTok{cbind}\NormalTok{(}\KeywordTok{mcols}\NormalTok{(snps_metadata)[}\KeywordTok{c}\NormalTok{(}\StringTok{"SNPS"}\NormalTok{, }\StringTok{"P-VALUE"}\NormalTok{)], }\KeywordTok{mcols}\NormalTok{(snps_anno))}
\end{Highlighting}
\end{Shaded}

We can visualise where these SNPs are located (Figure 6):

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{loc <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\KeywordTok{table}\NormalTok{(snps_anno}\OperatorTok{$}\NormalTok{LOCATION))}
\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data =}\NormalTok{ loc, }\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =} \KeywordTok{reorder}\NormalTok{(Var1, }\OperatorTok{-}\NormalTok{Freq), }\DataTypeTok{y =}\NormalTok{ Freq)) }\OperatorTok{+}
\StringTok{  }\KeywordTok{geom_bar}\NormalTok{(}\DataTypeTok{stat =} \StringTok{"identity"}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{xlab}\NormalTok{(}\StringTok{"Genomic location of SNPs"}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{ylab}\NormalTok{(}\StringTok{"Number of SNPs"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{figure}

{\centering \includegraphics{biocondutor-regulatory-genomics-workflow_files/figure-latex/figure6-1} 

}

\caption{Bar plot showing genomic locations associated with SLE variants.}\label{fig:figure6}
\end{figure}

As expected \citep{Maurano2012}, the great majority of SNPs are located within introns and in intergenic regions.
For the moment, we will focus on SNPs that are either coding or in promoter and UTR regions, as these can be assigned to target genes rather unambiguously:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{snps_easy <-}\StringTok{ }\KeywordTok{subset}\NormalTok{(snps_anno, LOCATION }\OperatorTok{==}\StringTok{ "coding"} \OperatorTok{|}\StringTok{ }\NormalTok{LOCATION }\OperatorTok{==}\StringTok{ "promoter"} \OperatorTok{|}\StringTok{ }\NormalTok{LOCATION }\OperatorTok{==}\StringTok{ "threeUTR"} \OperatorTok{|}\StringTok{ }\NormalTok{LOCATION }\OperatorTok{==}\StringTok{ "fiveUTR"}\NormalTok{)}
\NormalTok{snps_easy <-}\StringTok{ }\KeywordTok{as.data.frame}\NormalTok{(snps_easy)}
\end{Highlighting}
\end{Shaded}

Now we can check if any of the genes we found to be differentially expressed in SLE is also genetically associated with the disease:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{snps_easy_in_degs <-}\StringTok{ }\KeywordTok{merge}\NormalTok{(degs, snps_easy, }\DataTypeTok{by.x =} \StringTok{"gene_id"}\NormalTok{, }\DataTypeTok{by.y =} \StringTok{"GENEID"}\NormalTok{, }\DataTypeTok{all =} \OtherTok{FALSE}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

We have 14 genes showing differential expression in SLE that are also genetically associated with the disease.
While this is an interesting result, these hits are likely to be already well-known as potential SLE targets given their clear genetic association.

We will store essential information about these hits in a results \texttt{data.frame}:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{prioritised_hits <-}\StringTok{ }\KeywordTok{unique}\NormalTok{(}\KeywordTok{data.frame}\NormalTok{(}
  \DataTypeTok{snp_id =}\NormalTok{ snps_easy_in_degs}\OperatorTok{$}\NormalTok{SNPS,}
  \DataTypeTok{snp_pvalue =}\NormalTok{ snps_easy_in_degs}\OperatorTok{$}\NormalTok{P.VALUE,}
  \DataTypeTok{snp_location =}\NormalTok{ snps_easy_in_degs}\OperatorTok{$}\NormalTok{LOCATION,}
  \DataTypeTok{gene_id =}\NormalTok{ snps_easy_in_degs}\OperatorTok{$}\NormalTok{gene_id,}
  \DataTypeTok{gene_symbol =}\NormalTok{ snps_easy_in_degs}\OperatorTok{$}\NormalTok{symbol,}
  \DataTypeTok{gene_pvalue =}\NormalTok{ snps_easy_in_degs}\OperatorTok{$}\NormalTok{padj,}
  \DataTypeTok{gene_log2foldchange =}\NormalTok{ snps_easy_in_degs}\OperatorTok{$}\NormalTok{log2FoldChange,}
  \DataTypeTok{method =} \StringTok{"Direct overlap"}\NormalTok{,}
  \DataTypeTok{row.names =} \OtherTok{NULL}\NormalTok{))}
\KeywordTok{head}\NormalTok{(prioritised_hits, }\DecValTok{3}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##       snp_id snp_pvalue snp_location            gene_id gene_symbol
## 1  rs1887428      1e-06      fiveUTR ENSG00000096968.13        JAK2
## 2 rs58688157      5e-13     promoter ENSG00000099834.18       CDHR5
## 3  rs1990760      4e-08       coding  ENSG00000115267.5       IFIH1
##    gene_pvalue gene_log2foldchange         method
## 1 1.951160e-04            0.636590 Direct overlap
## 2 1.455662e-05            1.033372 Direct overlap
## 3 2.719420e-10            1.745324 Direct overlap
\end{verbatim}

\subsection{Use of regulatory genomic data to map intronic and intergenic SNPs to target genes}\label{use-of-regulatory-genomic-data-to-map-intronic-and-intergenic-snps-to-target-genes}

But what about all the SNPs in introns and intergenic regions?
Some of those might be regulatory variants affecting the expression level of their target gene(s) through a distal enhancer.
Let's create a dataset of candidate regulatory SNPs that are either intronic or intergenic and remove the annotation obtained with \texttt{VariantAnnotation}:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{snps_hard <-}\StringTok{ }\KeywordTok{subset}\NormalTok{(snps_anno, LOCATION }\OperatorTok{==}\StringTok{ "intron"} \OperatorTok{|}\StringTok{ }\NormalTok{LOCATION }\OperatorTok{==}\StringTok{ "intergenic"}\NormalTok{, }\DataTypeTok{select =} \KeywordTok{c}\NormalTok{(}\StringTok{"SNPS"}\NormalTok{, }\StringTok{"P.VALUE"}\NormalTok{, }\StringTok{"LOCATION"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\subsubsection{eQTL data}\label{eqtl-data}

A well-established way to gain insights into target genes of regulatory SNPs is to use eQTL data, where correlations between genetic variants and expression of genes are computed across different tissues or cell types \citep{Albert2015}.
Here, we will simply match GWAS SNPs and eQTLs according to their genomic locations, which is a rather crude way to integrate these two types of data.
More robust alternatives such as PrediXcan \citep{Gamazon2015}, TWAS \citep{Gusev2016} and SMR \citep{Zhu2016} exist and should be adopted if possible.
One downside of these methods is that they require subject-level or complete summary data, making them less practical in some circumstances.

We will use blood eQTL data from the GTEx consortium \citep{GTEx2017}.
To get the data, you will have to register and download the file \texttt{GTEx\_Analysis\_v7\_eQTL.tar.gz} from the \href{https://www.gtexportal.org}{GTEx portal} to the current working directory:

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# uncomment the following line to extract the gzipped archive file}
\CommentTok{#untar("GTEx_Analysis_v7_eQTL.tar.gz")}
\NormalTok{gtex_blood <-}\StringTok{ }\KeywordTok{read.delim}\NormalTok{(}\KeywordTok{gzfile}\NormalTok{(}\StringTok{"GTEx_Analysis_v7_eQTL/Whole_Blood.v7.signif_variant_gene_pairs.txt.gz"}\NormalTok{), }\DataTypeTok{stringsAsFactors =} \OtherTok{FALSE}\NormalTok{)}
\KeywordTok{head}\NormalTok{(gtex_blood, }\DecValTok{3}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##           variant_id           gene_id tss_distance ma_samples ma_count
## 1 1_231153_CTT_C_b37 ENSG00000223972.4       219284         13       13
## 2    1_61920_G_A_b37 ENSG00000238009.2       -67303         18       20
## 3    1_64649_A_C_b37 ENSG00000238009.2       -64574         16       16
##         maf pval_nominal    slope slope_se pval_nominal_threshold
## 1 0.0191740  3.69025e-08 1.319720 0.233538            1.35366e-04
## 2 0.0281690  7.00836e-07 0.903786 0.178322            8.26088e-05
## 3 0.0220386  5.72066e-07 1.110040 0.217225            8.26088e-05
##   min_pval_nominal   pval_beta
## 1      3.69025e-08 4.67848e-05
## 2      6.50297e-10 1.11312e-06
## 3      6.50297e-10 1.11312e-06
\end{verbatim}

We have to extract the genomic locations of the SNPs from the IDs used by GTEx:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{locs <-}\StringTok{ }\KeywordTok{strsplit}\NormalTok{(gtex_blood}\OperatorTok{$}\NormalTok{variant_id, }\StringTok{"_"}\NormalTok{)}
\NormalTok{gtex_blood}\OperatorTok{$}\NormalTok{chr <-}\StringTok{ }\KeywordTok{sapply}\NormalTok{(locs, }\StringTok{"["}\NormalTok{, }\DecValTok{1}\NormalTok{)}
\NormalTok{gtex_blood}\OperatorTok{$}\NormalTok{start <-}\StringTok{ }\KeywordTok{sapply}\NormalTok{(locs, }\StringTok{"["}\NormalTok{, }\DecValTok{2}\NormalTok{)}
\NormalTok{gtex_blood}\OperatorTok{$}\NormalTok{end <-}\StringTok{ }\KeywordTok{sapply}\NormalTok{(locs, }\StringTok{"["}\NormalTok{, }\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

We can then convert the \texttt{data.frame} into a \texttt{GRanges} object:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{gtex_blood <-}\StringTok{ }\KeywordTok{makeGRangesFromDataFrame}\NormalTok{(gtex_blood, }\DataTypeTok{keep.extra.columns =} \OtherTok{TRUE}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

We also need to ensure that the chromosome notation is consistent with the previous objects:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{seqlevelsStyle}\NormalTok{(gtex_blood)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] "NCBI"    "Ensembl"
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{seqlevelsStyle}\NormalTok{(gtex_blood) <-}\StringTok{ "UCSC"}
\end{Highlighting}
\end{Shaded}

From the publication \citep{GTEx2017}, we know the genomic coordinates are mapped to genome reference GRCh37, so we will have to uplift them to GRCh38 using \texttt{rtracklayer} \citep{Lawrence2009} and a mapping (``chain'') file.
The \href{https://cran.r-project.org/package=R.utils}{\texttt{R.utils}} package is only required to extract the gzipped file:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(rtracklayer)}
\KeywordTok{library}\NormalTok{(R.utils)}
\CommentTok{# uncomment the following line to download file}
\CommentTok{#download.file("http://hgdownload.cse.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz", destfile = "hg19ToHg38.over.chain.gz")}
\CommentTok{# uncomment the following line to extract gzipped file}
\CommentTok{#gunzip("hg19ToHg38.over.chain.gz")}
\NormalTok{ch <-}\StringTok{ }\KeywordTok{import.chain}\NormalTok{(}\StringTok{"hg19ToHg38.over.chain"}\NormalTok{)}
\NormalTok{gtex_blood <-}\StringTok{ }\KeywordTok{unlist}\NormalTok{(}\KeywordTok{liftOver}\NormalTok{(gtex_blood, ch))}
\end{Highlighting}
\end{Shaded}

We will use the \texttt{GenomicRanges} package \citep{Lawrence2013} to compute the overlap between GWAS SNPs and blood eQTLs:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(GenomicRanges)}
\NormalTok{hits <-}\StringTok{ }\KeywordTok{findOverlaps}\NormalTok{(snps_hard, gtex_blood)}
\NormalTok{snps_hard_in_gtex_blood =}\StringTok{ }\NormalTok{snps_hard[}\KeywordTok{queryHits}\NormalTok{(hits)]}
\NormalTok{gtex_blood_with_snps_hard =}\StringTok{ }\NormalTok{gtex_blood[}\KeywordTok{subjectHits}\NormalTok{(hits)]}
\KeywordTok{mcols}\NormalTok{(snps_hard_in_gtex_blood) <-}\StringTok{ }\KeywordTok{cbind}\NormalTok{(}\KeywordTok{mcols}\NormalTok{(snps_hard_in_gtex_blood), }\KeywordTok{mcols}\NormalTok{(gtex_blood_with_snps_hard))}
\NormalTok{snps_hard_in_gtex_blood <-}\StringTok{ }\KeywordTok{as.data.frame}\NormalTok{(snps_hard_in_gtex_blood)}
\end{Highlighting}
\end{Shaded}

We have 59 blood eQTL variants that are associated with SLE.
We can now check whether any of the genes differentially expressed in SLE is an \emph{eGene}, a gene whose expression is influenced by an eQTL.
Note that gene IDs in GTEx are mapped to GENCODE v19 \citep{GTEx2017}, while we are using the newer v25 for the DEGs.
To match the gene IDs in the two objects, we will simply strip the last bit containing the GENCODE gene version, which effectively gives us Ensembl gene IDs:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{snps_hard_in_gtex_blood}\OperatorTok{$}\NormalTok{ensembl_id <-}\StringTok{ }\KeywordTok{sub}\NormalTok{(}\StringTok{"(ENSG[0-9]+)}\CharTok{\textbackslash{}\textbackslash{}}\StringTok{.[0-9]+"}\NormalTok{, }\StringTok{"}\CharTok{\textbackslash{}\textbackslash{}}\StringTok{1"}\NormalTok{, snps_hard_in_gtex_blood}\OperatorTok{$}\NormalTok{gene_id)}
\NormalTok{degs}\OperatorTok{$}\NormalTok{ensembl_id <-}\StringTok{ }\KeywordTok{sub}\NormalTok{(}\StringTok{"(ENSG[0-9]+)}\CharTok{\textbackslash{}\textbackslash{}}\StringTok{.[0-9]+"}\NormalTok{, }\StringTok{"}\CharTok{\textbackslash{}\textbackslash{}}\StringTok{1"}\NormalTok{, degs}\OperatorTok{$}\NormalTok{gene_id)}
\NormalTok{snps_hard_in_gtex_blood_in_degs <-}\StringTok{ }\KeywordTok{merge}\NormalTok{(snps_hard_in_gtex_blood, degs, }\DataTypeTok{by =} \StringTok{"ensembl_id"}\NormalTok{, }\DataTypeTok{all =} \OtherTok{FALSE}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

We can add these 17 genes to our list:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{prioritised_hits <-}\StringTok{ }\KeywordTok{unique}\NormalTok{(}\KeywordTok{rbind}\NormalTok{(prioritised_hits, }\KeywordTok{data.frame}\NormalTok{(}
  \DataTypeTok{snp_id =}\NormalTok{ snps_hard_in_gtex_blood_in_degs}\OperatorTok{$}\NormalTok{SNPS,}
  \DataTypeTok{snp_pvalue =}\NormalTok{ snps_hard_in_gtex_blood_in_degs}\OperatorTok{$}\NormalTok{P.VALUE,}
  \DataTypeTok{snp_location =}\NormalTok{ snps_hard_in_gtex_blood_in_degs}\OperatorTok{$}\NormalTok{LOCATION,}
  \DataTypeTok{gene_id =}\NormalTok{ snps_hard_in_gtex_blood_in_degs}\OperatorTok{$}\NormalTok{gene_id.y,}
  \DataTypeTok{gene_symbol =}\NormalTok{ snps_hard_in_gtex_blood_in_degs}\OperatorTok{$}\NormalTok{symbol,}
  \DataTypeTok{gene_pvalue =}\NormalTok{ snps_hard_in_gtex_blood_in_degs}\OperatorTok{$}\NormalTok{padj,}
  \DataTypeTok{gene_log2foldchange =}\NormalTok{ snps_hard_in_gtex_blood_in_degs}\OperatorTok{$}\NormalTok{log2FoldChange,}
  \DataTypeTok{method =} \StringTok{"GTEx eQTLs"}\NormalTok{,}
  \DataTypeTok{row.names =} \OtherTok{NULL}\NormalTok{)))}
\end{Highlighting}
\end{Shaded}

\subsubsection{FANTOM5 data}\label{fantom5-data}

The FANTOM consortium profiled gene expression across a large panel of tissues and cell types using CAGE \citep{Fantom2014, Andersson2014}.
This technology allows mapping of transcription start sites and enhancer RNAs genome-wide.
Correlations between these promoter and enhancer elements across a large panel of tissues and cell types can then be calculated to identify significant promoter - enhancer pairs.
In turn, we will use these correlations to map distal regulatory SNPs to target genes.

Let's read in the enhancer - promoter correlation data:

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# uncomment the following line to download the file}
\CommentTok{#download.file("http://enhancer.binf.ku.dk/presets/enhancer_tss_associations.bed", destfile = "enhancer_tss_associations.bed")}
\NormalTok{fantom <-}\StringTok{ }\KeywordTok{read.delim}\NormalTok{(}\StringTok{"enhancer_tss_associations.bed"}\NormalTok{, }\DataTypeTok{skip =} \DecValTok{1}\NormalTok{, }\DataTypeTok{stringsAsFactors =} \OtherTok{FALSE}\NormalTok{)}
\KeywordTok{head}\NormalTok{(fantom, }\DecValTok{3}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##   X.chrom chromStart chromEnd
## 1    chr1     858252   861621
## 2    chr1     894178   956888
## 3    chr1     901376   956888
##                                                                  name
## 1                   chr1:858256-858648;NM_152486;SAMD11;R:0.404;FDR:0
## 2 chr1:956563-956812;NM_015658;NOC2L;R:0.202;FDR:8.01154668254404e-08
## 3     chr1:956563-956812;NM_001160184,NM_032129;PLEKHN1;R:0.422;FDR:0
##   score strand thickStart thickEnd itemRgb blockCount blockSizes
## 1   404      .     858452   858453   0,0,0          2   401,1001
## 2   202      .     956687   956688   0,0,0          2   1001,401
## 3   422      .     956687   956688   0,0,0          2   1001,401
##   chromStarts
## 1      0,2368
## 2     0,62309
## 3     0,55111
\end{verbatim}

Everything we need is in the fourth column, \texttt{name}: genomic location of the enhancer, gene identifiers, Pearson correlation coefficient and significance.
We will use the \href{https://cran.r-project.org/package=splitstackshape}{\texttt{splitstackshape}} package to parse it:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(splitstackshape)}
\NormalTok{fantom <-}\StringTok{ }\KeywordTok{as.data.frame}\NormalTok{(}\KeywordTok{cSplit}\NormalTok{(fantom, }\DataTypeTok{splitCols =} \StringTok{"name"}\NormalTok{, }\DataTypeTok{sep =} \StringTok{";"}\NormalTok{, }\DataTypeTok{direction =} \StringTok{"wide"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

Now we can extract the genomic locations of the enhancers and the correlation values:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{locs <-}\StringTok{ }\KeywordTok{strsplit}\NormalTok{(}\KeywordTok{as.character}\NormalTok{(fantom}\OperatorTok{$}\NormalTok{name_}\DecValTok{1}\NormalTok{), }\StringTok{"[:-]"}\NormalTok{)}
\NormalTok{fantom}\OperatorTok{$}\NormalTok{chr <-}\StringTok{ }\KeywordTok{sapply}\NormalTok{(locs, }\StringTok{"["}\NormalTok{, }\DecValTok{1}\NormalTok{)}
\NormalTok{fantom}\OperatorTok{$}\NormalTok{start <-}\StringTok{ }\KeywordTok{as.numeric}\NormalTok{(}\KeywordTok{sapply}\NormalTok{(locs, }\StringTok{"["}\NormalTok{, }\DecValTok{2}\NormalTok{))}
\NormalTok{fantom}\OperatorTok{$}\NormalTok{end <-}\StringTok{ }\KeywordTok{as.numeric}\NormalTok{(}\KeywordTok{sapply}\NormalTok{(locs, }\StringTok{"["}\NormalTok{, }\DecValTok{3}\NormalTok{))}
\NormalTok{fantom}\OperatorTok{$}\NormalTok{symbol <-}\StringTok{ }\NormalTok{fantom}\OperatorTok{$}\NormalTok{name_}\DecValTok{3}
\NormalTok{fantom}\OperatorTok{$}\NormalTok{corr <-}\StringTok{ }\KeywordTok{sub}\NormalTok{(}\StringTok{"R:"}\NormalTok{, }\StringTok{""}\NormalTok{, fantom}\OperatorTok{$}\NormalTok{name_}\DecValTok{4}\NormalTok{)}
\NormalTok{fantom}\OperatorTok{$}\NormalTok{fdr <-}\StringTok{ }\KeywordTok{sub}\NormalTok{(}\StringTok{"FDR:"}\NormalTok{, }\StringTok{""}\NormalTok{, fantom}\OperatorTok{$}\NormalTok{name_}\DecValTok{5}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

We can select only the enhancer - promoter pairs with a decent level of correlation and significance and tidy the data at the same time:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{fantom <-}\StringTok{ }\KeywordTok{unique}\NormalTok{(}\KeywordTok{subset}\NormalTok{(fantom, corr }\OperatorTok{>=}\StringTok{ }\FloatTok{0.25} \OperatorTok{&}\StringTok{ }\NormalTok{fdr }\OperatorTok{<}\StringTok{ }\FloatTok{1e-5}\NormalTok{, }\DataTypeTok{select =} \KeywordTok{c}\NormalTok{(}\StringTok{"chr"}\NormalTok{, }\StringTok{"start"}\NormalTok{, }\StringTok{"end"}\NormalTok{, }\StringTok{"symbol"}\NormalTok{)))}
\end{Highlighting}
\end{Shaded}

Now we would like to check whether any of our candidate regulatory SNPs are falling in any of these enhancers.
To do this, we have to convert the \texttt{data.frame} into a \texttt{GRanges} object and uplift the GRCh37 coordinates \citep{Fantom2014} to GRCh38:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{fantom <-}\StringTok{ }\KeywordTok{makeGRangesFromDataFrame}\NormalTok{(fantom, }\DataTypeTok{keep.extra.columns =} \OtherTok{TRUE}\NormalTok{)}
\NormalTok{fantom <-}\StringTok{ }\KeywordTok{unlist}\NormalTok{(}\KeywordTok{liftOver}\NormalTok{(fantom, ch))}
\end{Highlighting}
\end{Shaded}

We can now compute the overlap between SNPs and enhancers:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{hits <-}\StringTok{ }\KeywordTok{findOverlaps}\NormalTok{(snps_hard, fantom)}
\NormalTok{snps_hard_in_fantom =}\StringTok{ }\NormalTok{snps_hard[}\KeywordTok{queryHits}\NormalTok{(hits)]}
\NormalTok{fantom_with_snps_hard =}\StringTok{ }\NormalTok{fantom[}\KeywordTok{subjectHits}\NormalTok{(hits)]}
\KeywordTok{mcols}\NormalTok{(snps_hard_in_fantom) <-}\StringTok{ }\KeywordTok{cbind}\NormalTok{(}\KeywordTok{mcols}\NormalTok{(snps_hard_in_fantom), }\KeywordTok{mcols}\NormalTok{(fantom_with_snps_hard))}
\NormalTok{snps_hard_in_fantom <-}\StringTok{ }\KeywordTok{as.data.frame}\NormalTok{(snps_hard_in_fantom)}
\end{Highlighting}
\end{Shaded}

Let's check if any of these genes is differentially expressed in our RNA-seq data:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{snps_hard_in_fantom_in_degs <-}\StringTok{ }\KeywordTok{merge}\NormalTok{(snps_hard_in_fantom, degs, }\DataTypeTok{by =} \StringTok{"symbol"}\NormalTok{, }\DataTypeTok{all =} \OtherTok{FALSE}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

We have identified 7 genes whose putative enhancers contain SLE GWAS SNPs.
Let's add these to our list:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{prioritised_hits <-}\StringTok{ }\KeywordTok{unique}\NormalTok{(}\KeywordTok{rbind}\NormalTok{(prioritised_hits, }\KeywordTok{data.frame}\NormalTok{(}
  \DataTypeTok{snp_id =}\NormalTok{ snps_hard_in_fantom_in_degs}\OperatorTok{$}\NormalTok{SNPS,}
  \DataTypeTok{snp_pvalue =}\NormalTok{ snps_hard_in_fantom_in_degs}\OperatorTok{$}\NormalTok{P.VALUE,}
  \DataTypeTok{snp_location =}\NormalTok{ snps_hard_in_fantom_in_degs}\OperatorTok{$}\NormalTok{LOCATION,}
  \DataTypeTok{gene_id =}\NormalTok{ snps_hard_in_fantom_in_degs}\OperatorTok{$}\NormalTok{gene_id,}
  \DataTypeTok{gene_symbol =}\NormalTok{ snps_hard_in_fantom_in_degs}\OperatorTok{$}\NormalTok{symbol,}
  \DataTypeTok{gene_pvalue =}\NormalTok{ snps_hard_in_fantom_in_degs}\OperatorTok{$}\NormalTok{padj,}
  \DataTypeTok{gene_log2foldchange =}\NormalTok{ snps_hard_in_fantom_in_degs}\OperatorTok{$}\NormalTok{log2FoldChange,}
  \DataTypeTok{method =} \StringTok{"FANTOM5 correlations"}\NormalTok{,}
  \DataTypeTok{row.names =} \OtherTok{NULL}\NormalTok{)))}
\end{Highlighting}
\end{Shaded}

\subsubsection{Promoter Capture Hi-C data}\label{promoter-capture-hi-c-data}

More recently, chromatin interaction data was generated across 17 human primary blood cell types using promoter capture Hi-C \citep{Javierre2016}.
More than 30,000 promoter baits were used to capture promoter-interacting regions genome-wide, which were then mapped to enhancers based on annotation present in the Ensembl Regulatory Build \citep{Zerbino2015}.
This dataset provides a valuable resource for interpreting complex genomic data, especially in the context of autoimmune diseases (and other conditions where immune cells play a role).
Significant interactions between enhancers and promoters can be accessed in the supplementary data of the paper:

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# uncomment the following line to download file}
\CommentTok{#download.file("http://www.cell.com/cms/attachment/2086554122/2074217047/mmc4.zip", destfile = "mmc4.zip")}
\CommentTok{# uncomment the following lines to extract zipped files}
\CommentTok{#unzip("mmc4.zip")}
\CommentTok{#unzip("DATA_S1.zip")}
\NormalTok{pchic <-}\StringTok{ }\KeywordTok{read.delim}\NormalTok{(}\StringTok{"ActivePromoterEnhancerLinks.tsv"}\NormalTok{, }\DataTypeTok{stringsAsFactors =} \OtherTok{FALSE}\NormalTok{)}
\KeywordTok{head}\NormalTok{(pchic, }\DecValTok{3}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##   baitChr  baitSt baitEnd baitID oeChr    oeSt   oeEnd oeID
## 1    chr1 1206873 1212438    254  chr1  943676  957199  228
## 2    chr1 1206873 1212438    254  chr1 1034268 1040208  235
## 3    chr1 1206873 1212438    254  chr1 1040208 1043143  236
##                       cellType.s.
## 1                            nCD8
## 2 nCD4,nCD8,Mac0,Mac1,Mac2,MK,Mon
## 3     nCD4,nCD8,Mac0,Mac1,Mac2,MK
##                                                                                                     sample.s.
## 1                                                                                                    C0066PH1
## 2 S007DDH2,S007G7H4,C0066PH1,S00C2FH1,S00390H1,S001MJH1,S001S7H2,S0022IH2,S00622H1,S00BS4H1,S004BTH2,C000S5H2
## 3          S007DDH2,S007G7H4,C0066PH1,S00C2FH1,S00390H1,S001MJH1,S001S7H2,S0022IH2,S00622H1,S00BS4H1,S004BTH2
\end{verbatim}

We will use the \texttt{InteractionSet} package \citep{Lun2016}, which is specifically designed for the representation of chromatin interaction data.
We start by creating a \texttt{GInteractions} object:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(InteractionSet)}
\NormalTok{promoters <-}\StringTok{ }\KeywordTok{GRanges}\NormalTok{(}\DataTypeTok{seqnames =}\NormalTok{ pchic}\OperatorTok{$}\NormalTok{baitChr, }\DataTypeTok{ranges =} \KeywordTok{IRanges}\NormalTok{(}\DataTypeTok{start =}\NormalTok{ pchic}\OperatorTok{$}\NormalTok{baitSt, }\DataTypeTok{end =}\NormalTok{ pchic}\OperatorTok{$}\NormalTok{baitEnd))}
\NormalTok{enhancers <-}\StringTok{ }\KeywordTok{GRanges}\NormalTok{(}\DataTypeTok{seqnames =}\NormalTok{ pchic}\OperatorTok{$}\NormalTok{oeChr, }\DataTypeTok{ranges =} \KeywordTok{IRanges}\NormalTok{(}\DataTypeTok{start =}\NormalTok{ pchic}\OperatorTok{$}\NormalTok{oeSt, }\DataTypeTok{end =}\NormalTok{ pchic}\OperatorTok{$}\NormalTok{oeEnd))}
\NormalTok{pchic <-}\StringTok{ }\KeywordTok{GInteractions}\NormalTok{(promoters, enhancers)}
\end{Highlighting}
\end{Shaded}

As gene identifiers are not provided, we also have to map promoters to the respective genes so that we know which genes are regulated by which enhancers.
We can do this by using the \texttt{TxDb} object we previously built to extract positions of transcription start sites (TSSs) and then add the GENCODE gene IDs as metadata to the \texttt{pchic} object:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{tsss <-}\StringTok{ }\KeywordTok{promoters}\NormalTok{(txdb, }\DataTypeTok{upstream =} \DecValTok{0}\NormalTok{, }\DataTypeTok{downstream =} \DecValTok{1}\NormalTok{, }\DataTypeTok{columns =} \StringTok{"gene_id"}\NormalTok{)}
\NormalTok{hits <-}\StringTok{ }\KeywordTok{nearest}\NormalTok{(promoters, tsss)}
\NormalTok{pchic}\OperatorTok{$}\NormalTok{gene_id <-}\StringTok{ }\KeywordTok{unlist}\NormalTok{(tsss[hits]}\OperatorTok{$}\NormalTok{gene_id)}
\end{Highlighting}
\end{Shaded}

Next, we calculate the overlaps between SLE GWAS SNPs and enhancers (the \emph{second} region of the \texttt{GInteractions} object) :

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{hits <-}\StringTok{ }\KeywordTok{findOverlaps}\NormalTok{(snps_hard, pchic, }\DataTypeTok{use.region =} \StringTok{"second"}\NormalTok{)}
\NormalTok{snps_hard_in_pchic =}\StringTok{ }\NormalTok{snps_hard[}\KeywordTok{queryHits}\NormalTok{(hits)]}
\NormalTok{pchic_with_snps_hard =}\StringTok{ }\NormalTok{pchic[}\KeywordTok{subjectHits}\NormalTok{(hits)]}
\KeywordTok{mcols}\NormalTok{(snps_hard_in_pchic) <-}\StringTok{ }\KeywordTok{cbind}\NormalTok{(}\KeywordTok{mcols}\NormalTok{(snps_hard_in_pchic), }\KeywordTok{mcols}\NormalTok{(pchic_with_snps_hard))}
\NormalTok{snps_hard_in_pchic <-}\StringTok{ }\KeywordTok{as.data.frame}\NormalTok{(snps_hard_in_pchic)}
\end{Highlighting}
\end{Shaded}

We check if any of these enhancers containing SLE variants are known to putatively regulate genes differentially expressed in SLE:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{snps_hard_in_pchic_in_degs <-}\StringTok{ }\KeywordTok{merge}\NormalTok{(snps_hard_in_pchic, degs, }\DataTypeTok{by =} \StringTok{"gene_id"}\NormalTok{, }\DataTypeTok{all =} \OtherTok{FALSE}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

And finally we add these 13 genes to our list:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{prioritised_hits <-}\StringTok{ }\KeywordTok{unique}\NormalTok{(}\KeywordTok{rbind}\NormalTok{(prioritised_hits, }\KeywordTok{data.frame}\NormalTok{(}
  \DataTypeTok{snp_id =}\NormalTok{ snps_hard_in_pchic_in_degs}\OperatorTok{$}\NormalTok{SNPS,}
  \DataTypeTok{snp_pvalue =}\NormalTok{ snps_hard_in_pchic_in_degs}\OperatorTok{$}\NormalTok{P.VALUE,}
  \DataTypeTok{snp_location =}\NormalTok{ snps_hard_in_pchic_in_degs}\OperatorTok{$}\NormalTok{LOCATION,}
  \DataTypeTok{gene_id =}\NormalTok{ snps_hard_in_pchic_in_degs}\OperatorTok{$}\NormalTok{gene_id,}
  \DataTypeTok{gene_symbol =}\NormalTok{ snps_hard_in_pchic_in_degs}\OperatorTok{$}\NormalTok{symbol,}
  \DataTypeTok{gene_pvalue =}\NormalTok{ snps_hard_in_pchic_in_degs}\OperatorTok{$}\NormalTok{padj,}
  \DataTypeTok{gene_log2foldchange =}\NormalTok{ snps_hard_in_pchic_in_degs}\OperatorTok{$}\NormalTok{log2FoldChange,}
  \DataTypeTok{method =} \StringTok{"Promoter capture Hi-C"}\NormalTok{,}
  \DataTypeTok{row.names =} \OtherTok{NULL}\NormalTok{)))}
\end{Highlighting}
\end{Shaded}

These are the final results of our target identification exercise.
We can have a look at the most significant SNPs mapped with each of the methods:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{top_prioritised_hits <-}\StringTok{ }\NormalTok{prioritised_hits[}\KeywordTok{order}\NormalTok{(prioritised_hits}\OperatorTok{$}\NormalTok{snp_pvalue),]}
\NormalTok{top_prioritised_hits <-}\StringTok{ }\KeywordTok{split}\NormalTok{(top_prioritised_hits, top_prioritised_hits}\OperatorTok{$}\NormalTok{method)}
\KeywordTok{do.call}\NormalTok{(rbind, }\KeywordTok{lapply}\NormalTok{(top_prioritised_hits, head, }\DecValTok{1}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                          snp_id snp_pvalue snp_location            gene_id
## Direct overlap        rs3757387      1e-48     promoter ENSG00000128604.18
## GTEx eQTLs            rs1270942     2e-165       intron ENSG00000166278.14
## FANTOM5 correlations  rs1150754      6e-29       intron  ENSG00000204421.2
## Promoter capture Hi-C rs1270942     2e-165       intron  ENSG00000219797.2
##                       gene_symbol  gene_pvalue gene_log2foldchange
## Direct overlap               IRF5 5.006707e-03           0.4041349
## GTEx eQTLs                     C2 1.625111e-03           0.9269526
## FANTOM5 correlations       LY6G6C 3.575357e-05           1.4327915
## Promoter capture Hi-C          NA 1.919459e-04           0.4556364
##                                      method
## Direct overlap               Direct overlap
## GTEx eQTLs                       GTEx eQTLs
## FANTOM5 correlations   FANTOM5 correlations
## Promoter capture Hi-C Promoter capture Hi-C
\end{verbatim}

We can also visualise the relative contributions from the different approaches we used (Figure 7):

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{prioritised_genes <-}\StringTok{ }\KeywordTok{unique}\NormalTok{(}\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{gene_id =}\NormalTok{ prioritised_hits}\OperatorTok{$}\NormalTok{gene_id, }\DataTypeTok{method =}\NormalTok{ prioritised_hits}\OperatorTok{$}\NormalTok{method))}
\KeywordTok{ggplot}\NormalTok{(}\DataTypeTok{data =}\NormalTok{ prioritised_genes, }\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ method)) }\OperatorTok{+}
\StringTok{  }\KeywordTok{geom_bar}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{fill =}\NormalTok{ method), }\DataTypeTok{stat =} \StringTok{"count"}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{ylab}\NormalTok{(}\StringTok{"Number of genes"}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{theme}\NormalTok{(}\DataTypeTok{axis.title.x =} \KeywordTok{element_blank}\NormalTok{(),}
        \DataTypeTok{axis.text.x =} \KeywordTok{element_blank}\NormalTok{(),}
        \DataTypeTok{axis.ticks.x =} \KeywordTok{element_blank}\NormalTok{())}
\end{Highlighting}
\end{Shaded}

\begin{figure}

{\centering \includegraphics{biocondutor-regulatory-genomics-workflow_files/figure-latex/figure7-1} 

}

\caption{Bar plot showing number of genes identified by each variant mapping method.}\label{fig:figure7}
\end{figure}

We observe that all methods significantly contributed to the identification of genes associated with GWAS SNPs.
The majority of genes were identified through the integration of the GTEx blood eQTL data, followed by the methods based on direct overlap, promoter capture Hi-C data and FANTOM5 correlations.

\subsection{Functional analysis of prioritised hits}\label{functional-analysis-of-prioritised-hits}

We will use biological processes from the Gene Ontology \citep{Ashburner2000} and the \texttt{clusterProfiler} package \citep{Yu2012} to functionally characterise our list of genes:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(clusterProfiler)}
\NormalTok{prioritised_hits_ensembl_ids <-}\StringTok{ }\KeywordTok{unique}\NormalTok{(}\KeywordTok{sub}\NormalTok{(}\StringTok{"(ENSG[0-9]+)}\CharTok{\textbackslash{}\textbackslash{}}\StringTok{.[0-9]+"}\NormalTok{, }\StringTok{"}\CharTok{\textbackslash{}\textbackslash{}}\StringTok{1"}\NormalTok{, prioritised_hits}\OperatorTok{$}\NormalTok{gene_id))}
\NormalTok{all_genes_ensembl_ids <-}\StringTok{ }\KeywordTok{unique}\NormalTok{(}\KeywordTok{sub}\NormalTok{(}\StringTok{"(ENSG[0-9]+)}\CharTok{\textbackslash{}\textbackslash{}}\StringTok{.[0-9]+"}\NormalTok{, }\StringTok{"}\CharTok{\textbackslash{}\textbackslash{}}\StringTok{1"}\NormalTok{, }\KeywordTok{rownames}\NormalTok{(rse)))}
\NormalTok{gobp_enrichment <-}\StringTok{ }\KeywordTok{enrichGO}\NormalTok{(prioritised_hits_ensembl_ids,}
                            \DataTypeTok{universe =}\NormalTok{ all_genes_ensembl_ids,}
                            \DataTypeTok{OrgDb =}\NormalTok{ org.Hs.eg.db,}
                            \DataTypeTok{keyType =} \StringTok{"ENSEMBL"}\NormalTok{,}
                            \DataTypeTok{ont =} \StringTok{"BP"}\NormalTok{,}
                            \DataTypeTok{pAdjustMethod =} \StringTok{"BH"}\NormalTok{,}
                            \DataTypeTok{pvalueCutoff =} \FloatTok{0.05}\NormalTok{,}
                            \DataTypeTok{qvalueCutoff =} \FloatTok{0.05}\NormalTok{,}
                            \DataTypeTok{readable =} \OtherTok{TRUE}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

We can visualise the most enriched terms (Figure 8):

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{dotplot}\NormalTok{(gobp_enrichment, }\DataTypeTok{showCategory =} \DecValTok{20}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{figure}

{\centering \includegraphics{biocondutor-regulatory-genomics-workflow_files/figure-latex/figure8-1} 

}

\caption{Dot plot showing enrichment of Gene Ontology biological processes for the list of prioritised genes.}\label{fig:figure8}
\end{figure}

We observe a significant enrichment for interferon responses, antigen processing and presentation, and T cell stimulation, all processes which are well-known to play key roles in the pathogenesis of SLE \citep{Oon2016, Morris2014, Suarez-Fueyo2016}.

From a drug discovery perspective, JAK2 is probably the most attractive target: rs1887428 (\emph{p}-value = 1 x 10\textsuperscript{-6}) is located in its 5' UTR and the genes is significantly upregulated in disease.
Tofacitinib, a pan-JAK inhibitor, showed promising results in mouse \citep{Furumoto2017} and is currently being tested or safety in a \href{https://clinicaltrials.gov/ct2/show/NCT02535689}{phase I clinical trial}.
We find 7 GWAS SNPs that are blood eQTLs linked to the expression of C2, a protease active in the complement signalling cascade.
The most significant variant is rs1270942 (\emph{p}-value = 2 x 10\textsuperscript{-165}) and is found in an intron of CFB, another component of the complement system.
As with other autoimmune diseases, the complement plays a key role SLE in and has been investigated as a therapeutic approach \citep{Leffler2014}.
Another potentially interesting hit is TAX1BP1: rs849142 (\emph{p}-value = 1 x 9\textsuperscript{-11}) is found within an intron of JAZF1, but can be linked to TAX1BP1 via a chromatin interaction with its promoter.
TAX1BP1 inhibits TNF-induced apoptosis \citep{DeValck1999} and is involved in the IL1 signalling cascade \citep{Ling2000}, another relevant pathway in SLE that could be therapeutically targeted \citep{Ronnblom2010}.

\section{Conclusions}\label{conclusions}

In this Bioconductor workflow we have used several packages and datasets to demonstrate how regulatory genomic data can be used to annotate significant hits from GWASs and prioritise gene lists from expression studies, providing an intermediate layer connecting genetics and transcriptomics.
Overall, we identified 46 SLE-associated SNPs that we mapped to 49 genes differentially expressed in SLE, using eQTL data \citep{GTEx2017} and enhancer - promoter relationships from CAGE \citep{Fantom2014} and promoter capture Hi-C experiments \citep{Javierre2016}.
These genes are involved in key inflammatory signalling pathways and some of them could develop into therapeutic targets for SLE.

The workflow also demonstrates some real-world challenges encountered when working with genomic data from different sources, such as the use of different genome assemblies and gene annotation systems, the parsing of files with custom formats into Bioconductor objects and the mapping of genomic locations to genes.
While options for the visualisations of genomic data and interactions are outside the scope of this workflow, at least three good alternatives exist in Bioconductor: ggbio \citep{Yin2012}, Sushi \citep{Phanstiel2014} and Gviz \citep{Hahne2016} coupled with the GenomicInteractions package \citep{Harmston2015}.
We refer the reader to these publications and package vignettes for examples.

As the sample size and power of GWASs and gene expression studies continue to increase, it will become more and more challenging to identify truly significant hits and interpret them.
The use of regulatory genomics data as presented here can be an important tool to gain insights into large biomedical datasets and help in the identification of biomarkers and therapeutic targets.

\section{Abbreviations}\label{abbreviations}

CAGE: cap analysis of gene expression\\
DHS: DNase I hypersensitive site\\
eQTL: expression quantitative trait locus\\
GWAS: genome-wide association study\\
ISM: interferon signature metric\\
SLE: systemic lupus erythematosus\\
SNP: single nucleotide polymorphism

\section{Data and software availability}\label{data-and-software-availability}

Download links for all datasets are part of the workflow.\\
Software packages required to reproduce the analysis can be installed as part of the workflow.\\
Source code is available at: \url{https://github.com/enricoferrero/bioconductor-regulatory-genomics-workflow}.\\
Archived source code as at the time of publication is available at: \url{https://doi.org/10.5281/zenodo.1154124}.

\section{Competing interests}\label{competing-interests}

EF is a full time employee of GSK.

\section{Grant information}\label{grant-information}

The author declared that no grants were involved in supporting this work.

\renewcommand\refname{References}
{\small\bibliography{references.bib}}

\end{document}
