\documentclass{article}
\usepackage[affil-it]{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
\usepackage{url}
\usepackage{hyperref}
\hypersetup{colorlinks=false,pdfborder={0 0 0}}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[english]{babel}




\begin{document}

\title{Supplemental Material for WASP: allele-specific methods for unbiased discovery of molecular quantitative trait loci}


\author{Graham McVicker}
\affil{Affiliation not available}

  
\author{Bryce van de Geijn}
\affil{Affiliation not available}
  



\date{\today}

\maketitle




\documentclass{article}
%\usepackage{scicite}
\bibliographystyle{naturemag}
\usepackage{times}
\usepackage{graphicx}
\usepackage{url}
\usepackage{amsmath}
\usepackage{lscape}

\DeclareMathOperator{\Like}{L}
\DeclareMathOperator{\logit}{logit}
\DeclareMathOperator{\expit}{expit}

\renewcommand\thefigure{S\arabic{figure}}
\renewcommand\thetable{S\arabic{table}}


\topmargin 0.0cm
\oddsidemargin 0.2cm
\textwidth 16cm
\textheight 21cm
\footskip 1.0cm

\title{Supplemental Material for WASP: allele-specific software for unbiased discovery of molecular quantitative trait loci}

\author{Bryce van de Geijn$^{1,2,\dagger}$, Graham McVicker$^{3,\dagger}$,
Yoav Gilad$^1$, Jonathan K. Pritchard$^{3,4,5,\ast}$\\  \\
\normalsize{$^{1}$Department of Human Genetics, University of Chicago}\\
\normalsize{$^{2}$Committee on Genetics, Genomics and Systems Biology, University of Chicago}\\
\normalsize{$^{3}$Department of Genetics, Stanford University}\\
\normalsize{$^{4}$Department of Biology, Stanford University}\\
\normalsize{$^{5}$Howard Hughes Medical Institute, Stanford University}\\
\normalsize{$^\dagger$These authors contributed equally.}\\ 
\normalsize{$^\ast$To whom correspondence should be addressed: pritch@stanford.edu}\\
}
\date{}


\section{Unbiased read mapping with WASP}

To detect differences in molecular phenotypes from sequencing data it is essential to remove read mapping biases, which are a major source of false positives. The WASP read mapping pipeline accomplishes this task by ensuring that the mapping of each individual read is unbiased.

\subsection{How WASP mapping works}

In the WASP mapping pipeline, the user first maps reads to the genome using any mapper that outputs BAM or SAM format. For example, ChIP-seq reads can be mapped by BWA \cite{Li_2009} or Bowtie 2 \cite{Langmead_Salzberg_2012}, and RNA-seq reads can be mapped using tophat \cite{Trapnell_2009}. WASP then identifies mapped reads that overlap with known polymorphisms. For each read that overlaps a polymorphism, all possible allelic combinations that differ from the original read are generated and re-mapped to the genome. For example, when a read overlaps two bi-allelic SNPs four allelic combinations are possible, three of which differ from the original read. The original read is discarded if any of the allelic combinations map non-uniquely or map to another location. 

This simple method has the advantages that it works with almost any existing mapping pipeline and it handles reads with sequencing errors, which are a major source of biased mapping \cite{Degner_2009}.



\subsection{Comparing WASP mapping to N-masked and personal genome mapping}

% TODO: verify that IMPUTE ref panel info correct
To test the accuracy of allelic mapping using WASP, we simulated 100 bp reads from a lymphoblastoid cell line (GM12878) that has been genotyped by the 1000 Genomes and HapMap projects. We additionally imputed and phased genotypes for this cell line with IMPUTE2 \cite{Howie_Donnelly_Marchini_2009} using the 1000 Genomes Phase1 integrated version 3 reference panel. 

For each test, we evaluated the performance of WASP compared to mapping to a personal or N-masked genome. To create an N-masked genome, we created a copy of the hg19 genome with Ns in place of known variants from the GM12878 cell line. We similarly created maternal and paternal copies of GM12878 using the phased genotypes.  We mapped the simulated reads to the original, N-masked, and personal versions of the hg19 genome with BWA \cite{Li_2009} allowing up to 2 mismatches per read ($\verb|-n 2|$), and excluding gapped alignments ($\verb|-o 0|$).  For the personal genome, we kept a read if it mapped uniquely to either genome copy.  If it mapped to both genomes, we kept the location with the highest mapping quality (ties were broken randomly). 

\subsubsection{Quantifying fraction of reads showing imbalance}

We first identified each base where a read starting at that base would overlap a heterozygous site. We generated reads from each haplotype while introducing identical sequencing errors at a predefined rate. For each mapping type, we considered the mapping of a read to be biased if the read from one haplotype mapped to the correct location but the other did not.

\subsubsection{Assessing the effects of mapping bias on an allele specific study}

For each heterozygous site, we simulated 100 reads from random bases that overlap the chosen SNP.  We chose the haplotype of each simulated read at random.  Reads from peaks without effects came from haplotype 1 vs haplotype 2 with a 1:1 ratio.  Reads from peaks with effects were simulated with ratios ranging from 1.3:1 to 2.5:1 to test a range of effect sizes.

For each effect size, we simulated sets of peaks that were composed of $90\%$ null peaks and $10\%$ peaks with effects.  We mapped the sets with each mapping scheme and performed a binomial test for imbalance on each peak, calling significantly imbalanced loci at 10\% FDR.  We finally assessed the fraction of significant loci that came from the null peaks.  If there is no mapping induced imbalance, this should be ~10\%.

\section{Discovery of quantitative trait loci with WASP}

To discover molecular quantitative trait loci (QTLs) WASP uses a statistical test, which we call the Combined Haplotype Test (CHT). As input, CHT takes genotype probabilities at known SNPs as well as mapped reads from sequencing-based experiments such as ChIP-seq or RNA-seq. CHT combines two types of information: the depth of mapped reads and the allelic imbalance of mapped reads that overlap heterozygous sites. CHT models the overdispersion of read counts (both across regions and across individuals) and accounts for variability introduced by GC content and the fraction of reads that fall within peaks. 

We used an early version of CHT to discover histone mark QTLs \cite{McVicker_2013}. We have since extended the statistical model and improved the software implementation. These improvements include:
\begin{itemize}
\item Automatic estimation of dispersion parameters from the data,
\item Adjustment for unknown covariates by allowing principal component loadings to be provided,
\item Allowing tested regions to be split across multiple genomic segments, such as exons,
\item Greater efficiency so the model can now be run with hundreds of individuals.
\end{itemize}\selectlanguage{english}
\begin{table}[htdp]
\caption{{Description of mathematical variables used in the combined haplotype test}}
%\begin{center}
\small
\begin{tabular}{l l}
\textbf{Variable} & \textbf{Description}\\
 & \textbf{Index variables}\\
$h$ & test number (one per test SNP / target region pair) \\
$i$ & individual\\
$j$ & target region\\
$k$ & SNP within target region\\
$m$ & test SNP\\
& \textbf{Latent variables}\\
$\alpha_h$ & molecular phenotype level of the reference allele for test $h$\\
$\beta_h$ & molecular phenotype level of the alternative allele for test $h$\\
$p_h$ & fraction of allele-specific reads expected from reference allele $(p_h = \frac{\alpha_h}{\alpha_h + \beta_h})$\\
$T^{*}_{i,j}$ & genotype-independent expected total read count for individual $i$, target region $j$\\
$\lambda_{hi}$ & expected total read count for test $h$, individual $i$\\
$\Omega_i$ & overdispersion of read counts for individual $i$ (across all tests)\\
$\phi_j$ & overdispersion of read counts for test number $j$ (across all individuals)\\
$\Upsilon_i$ & overdispersion of allele-specific reads for individual $i$\\
& \textbf{Observed variables}\\
$x_{ij}$ & number of reads for individual $i$, target region $j$\\
$G_{im}$ & genotype call for individual $i$, test SNP $m$\\ 
$T_{i}$ & total number of genome-wide mapped reads for individual $i$\\
$n_{ik}$ & total number of allele-specific reads for individual $i$, target SNP $k$\\
$y_{ik}$ & number of allele-specific reads from reference haplotype for individual $i$, target SNP $k$\\
$H_{ik}$ & probability individual $i$ is heterozygous for target SNP $k$\\
%\hline \hline
\end{tabular}
%\end{center}
\label{tab:variables}
\end{table}%\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.42\columnwidth]{figures/WASP-mapping-pipeline/WASP-mapping-pipeline}

\end{center}
\end{figure}

\section{The combined haplotype test}

The combined haplotype test (CHT) determines whether the genotype of a ``test SNP'', $m$, is associated with read depth and allelic imbalance within a nearby ``target region'', $j$, on the same chromosome. A single target region may be discontiguous and span multiple genomic loci. For example, the exons of a gene can be used as a single target region, which is useful when searching for expression QTLs using RNA-seq reads. The test SNP is not required to be within the target region, but is assumed to be nearby and cis-acting. This allows us to combine information from across phased heterozygous SNPs and assign reads to one haplotype or the other.


\subsection{Basic model}

The CHT is a likelihood ratio test with two components: one to model the depth and the other to model the allelic imbalance of the mapped reads. Both components of the test are parameterized by $\alpha_h$ and $\beta_h$, which define the expected read depth from chromosomes with the reference and alternative alleles.  Since variants are assumed to be additive and cis-acting, the expected allelic imbalance in heterozygotes, $p_h$, is $p_h = \frac{\alpha_h}{\alpha_h + \beta_h}$ \cite{Sun_2011}.


\subsection{Modeling the read depths}

The number of reads mapping to a target region is often modelled using a poisson distribution \cite{Marioni_2008}.  However, the poisson assumption that the variance is equal to the mean is often violated because read counts from target regions are overdispersed. Part of this overdispersion can be accommodated by modelling the data with a negative-binomial distribution with a variance parameter for each test \cite{Anders_2010}. However, the negative binomial distribution assumes that the mean and variance have a quadratic relationship that is consistent across individuals. We have found that this assumption is violated by sequencing data and causes poor calibration of the tests, particularly when sample sizes are small, due to technical differences between experiments.  The CHT therefore models a negative binomial overdispersion parameter, $\Omega_i$, for each individual.  It also adds an additional overdispersion parameter, $\phi_h$, which models the variance at each site.  After adding this additional dispersion parameter, the data are modelled with a beta-negative-binomial (BNB) distribution.  The expected number of counts, $\lambda_{hi}$, is calculated based on $\alpha_h$, $\beta_h$, and the genotype, $G_{im}$ for individual $i$ at test SNP $m$.  The estimate is scaled by the total number of mapped reads for individual $i$, $T_i$.

\begin{equation}
\lambda_{hi} = \left\{
\begin{array}{ll}
2 \alpha_h T_i & \textrm{if } G_{im} = 0 \textrm{ (homozygous allele 1)} \\
\\
\left( \alpha_h + \beta_h \right) T_i & \textrm{if } G_{im} = 1 \textrm{ (heterozygous)} \\
\\
2 \beta_h T_i & \textrm{if } G_{im} = 2 \textrm{ (homozygous allele 2)}
\end{array} \right.
\label{eqn:read_depth_model}
\end{equation}

The likelihood of the data is then given by the equation

% Is it worth writing out the full PDF, since the parameterization of the BNB may not be clear 
% to people?

\[
\textrm{L}\left( D_h \left| \alpha_h, \beta_h, \Omega_\bullet, \phi_j \right. \right) = \prod_i \Pr_{\mathrm{BNB}} \left( X = x_{ij} \left| \lambda_{hi}, \Omega_i, \phi_j \right. \right) \\
\]

where $x_{ij}$ is the number of reads for individual $i$ in target region $j$.

\subsection{Correcting for GC content and other effects on expected read depth}

Since the number of mapped reads can differ between sequencing lanes and runs, we initially model the expected number of counts, $\lambda_{hi}$, as a linear function of the total number of mapped reads for each individual, $T_i$.  However, technical variation between experiments can change this relationship and reduce the power to detect true differences in read depths between samples or cause spurious associations. We directly model some known sources of technical variation and estimate adjusted total read depths, $T^{*}_{ij}$, for each individual and target region. We then replace $T_i$ in equation \ref{eqn:read_depth_model} with $T^{*}_{ij}$. This gives us a more accurate estimate of the expected number of reads and boosts our ability to detect true QTLs.

\subsubsection{Adjusting total read depth}

In RNA-seq experiments, a large fraction of mapped reads can come from a small number of highly expressed genes. Variation in the expression level of these genes can therefore have a large effect on the number of reads that map to all other genes \cite{Wagner_2012}. In ChIP-seq experiments, the fraction of reads that come from peaks varies between experiments, likely due to differences in the antibody efficacy.

To account for these types of variation, we calculate an adjusted total read depth, $T^{*}_{ij}$ for each region and individual. The adjusted read depth is defined by a quartic function of the total read depth (summed across individuals) for each target region. The coefficients of the quartic function are estimated using a maximum likelihood approach described below.

\subsubsection{GC content}

GC content also affects read depth, with a relationship that varies across samples \cite{Pickrell_2010,Benjamini_2012}. For example, in some samples, high GC content regions have high read depth, while in other samples they have low read depth. To account for this variation, we add GC content terms to the model of adjusted total read depth. These terms are modeled with a log linker so that $T^{*}_{ij}$ is guaranteed to be positive.  After fitting this model we can calculate an adjusted total read depth for each region that takes into account both the GC content variation and the total read depth variation.

\subsubsection{Fitting adjustment coefficients}

For each target region $j$, we count the total number of reads $v_j = \sum_i x_{ij}$ and calculate the GC content $w_j$.  Then, for each individual $i$, we fit coefficients $a_{0i}, a_{1i}, \ldots, b_{4i} $ that maximize the likelihood of the observed counts given the adjusted expected counts, $T^{*}_{ij}$:

\[
\textrm{L}\left( D_i \left| a_{0i}, a_{1i}, \ldots, b_{4i} \right. \right) = \prod_j \Pr_{\mathrm{Pois}} \left(X_{ij} = x_{ij} \left| T^{*}_{ij} \right. \right)\\
\]
\[
T^{*}_{ij} = \exp\left(a_{0i} + a_{1i} w_j  +   a_{2i} w_j^2 +  a_{3i} w_j^3  + a_{4i} w_j^4 \right) \left(b_{1i} v_j + b_{2i} v_j^2  + b_{3i} v_j^3  + b_{4i} v_j^4 \right)
\]\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/all-examples/all-examples}
\caption{{Replace this text with your caption%
}}
\end{center}
\end{figure}

\subsection{Modeling the allelic imbalances}

Allele-specific read counts are sometimes modelled using the binomial distribution \cite{Reddy_2012}, however, we have found that allele-specific read counts are overdispersed.  We instead model allele-specific read counts with a beta-binomial (BB) distribution and include a parameter $\Upsilon_i$ (estimated separately) that captures the overdispersion for each individual. The likelihood of the parameters given the data is then:

\[
\textrm{L}\left(\alpha_h, \beta_h \left| D \right. \right) = \prod_i \prod_k \Pr_{\mathrm{BB}} \left( Y = y_{ik} \left| n_{ik}, p_h, \Upsilon_i \right. \right) \\
\]

where $y_{ik}$ is the number of allele-specific reads from the reference haplotype and $n_{ik}$ is the total number of allele-specific reads for individual $i$ at target SNP $k$. The expected fraction of allele-specific reads from the reference allele is $p_h = \frac{\alpha_h}{\alpha_h + \beta_h}$. 

\subsection{Correcting for incorrect genotype calls}

SNP genotypes that are incorrectly called as heterozygous are a major source of false positives, since reads that overlap them appear to come from only one allele. To account for this issue, we assume that allele specific reads are drawn from a mixture of two beta-binomials, with probabilities $H_{ik}$ and $1-H_{ik}$, where $H_{ik}$ is the probability that individual $i$ is heterozygous for SNP $k$. Reads from heterozygous individuals contain the reference allele with probability $p_{h}$. We assume that reads from homozygous individuals still have a small probability of coming from the other allele due to sequencing errors, which occur with probability, $p_{\textrm{err}}$. The probability of observing $y_{ik}$ reads from the reference allele for individual $i$ at SNP $k$ then becomes:

\begin{eqnarray*}
& \Pr_{\mathrm{BB-mix}}\left(Y = y_{ik} \left| p_{h}, n_{ik}, \Upsilon_i \right. \right) =
H_{ik} \Pr_{\mathrm{BB}} \left(Y = y_{ik} \left| p_{h}, n_{ik}, \Upsilon_i \right. \right) &\\
& + (1 - H_{ik}) 
\left[ \Pr_{\mathrm{BB}} \left(Y = y_{ik} \left| p_{\textrm{err}}, n_{ik}, \Upsilon_i \right. \right) 
+ \Pr_{\mathrm{BB}} \left(Y = y_{ik} \left| 1-p_{\textrm{err}}, n_{ik}, \Upsilon_i \right. \right) \right] &
\end{eqnarray*}

We found that even SNPs with heterozygous probabilities of 1.0 are occasionally miscalled so we set heterozygous probabilities to a maximum value of 0.99. We then update this heterozygous probability using sequencing data obtained from the same individual. Sequencing data may consist of DNA sequencing reads or reads aggregated across multiple types of experiments performed on the same individual (e.g. RNA-seq and ChIP-seq reads).

For a SNP with heterozygous probability $H_{ik} = \min(0.99, H_{ik}^{\textrm{obs}})$, we define the updated heterozygous probability, $\hat{H}_{ik}$ as:

\[
\hat{H}_{ik} = \frac{H_{ik} \Pr_{\mathrm{Bin}} \left( D \left| p=0.5 \right. \right)}
{H_{ik} \Pr_{\mathrm{Bin}} \left( D \left| p=0.5 \right. \right)
+ (1 - H_{ik}) \left[ \Pr_{\mathrm{Bin}} \left( D \left| p=p_{err} \right. \right) + \Pr_{\mathrm{Bin}} \left( D \left| p=1-p_{\textrm{err}} \right. \right) \right]}
\]



\subsection{The Combined Likelihood Ratio Test}

The combined likelihood of both components of the model is:

\[
\textrm{L}\left(\alpha_h, \beta_h, \phi_j \left| D \right. \right) =  \prod_i \left[  \Pr_{\mathrm{BNB}} \left( X = x_{ij} \left| \lambda_{hi}, \Omega_i, \phi_j \right. \right) \prod_k \Pr_{\mathrm{BB-mix}}\left(Y = y_{ik} \left| p_{h}, n_{ik}, \Upsilon_i \right. \right) \right]
\]

To test for an association with genotype we perform a likelhood ratio test that compares the alternative hypothesis $\alpha_h \neq \beta_h$ to the null hypothesis $\alpha_h = \beta_h$. The CHT returns a likelihood ratio statistic $\Lambda = \frac{L(\hat{\theta}_1|D)}{L(\hat{\theta}_0|D)}$ where $\hat{\theta}_1$ and $\hat{\theta}_0$ are maximum likelihood estimates of the parameters under the alternative and null hypotheses. P-values can be calculated from the the test statistic under the asymptotic assumption that $-2log(\Lambda)$ is $\chi^2$ distributed with one degree of freedom.

\subsection{Estimating overdispersion parameters}

In order to estimate the genome-wide overdispersion parameters $\Phi_i$ and $\Upsilon_i$, we use the same likelihood equations as in the CHT, but assume that there are no genetic effects.  This means that for the read depth part of the test, $\lambda_{hi}$, is equal to the expected counts, $T^*_{ij}$, and for the allele-specific part of the test, $p_h$ is equal to $0.5$. Since the allele specific and read depth parts of the likelihood equation are independent, we can fit the overdispersion parameters separately.

\subsubsection{Beta-Negative-Binomial}

To find the maximum likelihood estimate of $\Omega_i$ we need to sum the log likelihood across all regions. This presents a problem, as $\phi_j$ must also be estimated for each region. We therefore interatively estimate $\phi_j$ by first finding a maximum likelihood estimate for $\phi_j$ for each region using the equation:

\[
\textrm{L}\left(\phi_j \left| D \right. \right) =  \prod_i \left[  \Pr_{\mathrm{BNB}} \left( X = x_{ij} \left| \lambda = T^*_{ij}, \Omega_i, \phi_j \right. \right) \right] 
\]

and then finding a maximum likelihood estimate for $\Omega_i$ for each individual using the equation:

\[
\textrm{L}\left(\Omega_i \left| D \right. \right) =  \prod_j \left[  \Pr_{\mathrm{BNB}} \left( X = x_{ij} \left| \lambda = T^*_{ij}, \Omega_i, \phi_j \right. \right) \right] 
\]

We repeat this iterative procedure until the improvement in likelihood becomes negligable.


\subsubsection{Beta-Binomial}

We calculate the genome-wide likelihood of $\Upsilon_i$ by taking the product of likelihoods from all target region SNPs that are heterozygous in individual $i$. We again assume there is no genetic effect, so $p$ = 0.5, and we use the following equation to find the maximum likelihood estimate of $\Upsilon_i$:

\[
\textrm{L}\left(\Upsilon_i \left| D \right. \right) = \prod_k \Pr_{\mathrm{BB-mix}} \left( Y = y_{ik} \left| n_{ik}, p = 0.5, \Upsilon_i \right. \right) \\
\]

\subsection{Correcting for unknown covariates using principal components}

Both known and unknown covariates such as time of experiment, age of sample, etc. can affect molecular trait measurements and confound QTL studies.  Principal component analysis (PCA) is sometimes used to capture and remove these effects \cite{Stephens_Gilad_Pritchard_2010,Degner_2012}. To leverage PCA while maintaining the discrete nature of the count data, the CHT directly models the covariate effects. To do this we include a user-defined number of PCA loadings $u_{i\bullet}$ and fit coefficients $c_{h\bullet}$ when calculating $\lambda_{hi}$.

\[
\lambda_{hi} = \left\{
\begin{array}{ll}
2 \alpha_h (1 + c_{h1} u_{i1} + c_{h2} u_{i2} + \ldots) T_i & \textrm{if } G_{im} = 0 \textrm{ (homozygous allele 1)} \\
\\
\left( \alpha_h + \beta_h \right) (1 + c_{h1} u_{i1} + c_{h2} u_{i2} + \ldots) T_i & \textrm{if } G_{im} = 1 \textrm{ (heterozygous)} \\
\\
2 \beta_h (1 + c_{h1} u_{i1} + c_{h2} u_{i2} + \ldots) T_i & \textrm{if } G_{im} = 2 \textrm{ (homozygous allele 2)}
\end{array} \right.
\]

Fitting many coefficients simultaneously can be quite slow, but since the principal components are by definition orthogonal, we can optimize their coefficients one at a time without losing accuracy. We then use the fitted coefficients to calculate $\lambda_{hi}$ for the null and alternative models.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/WASP-H3K27ac-permuted/WASP-H3K27ac-permuted}
\caption{{Replace this text with your caption%
}}
\end{center}
\end{figure}

\subsection{CHT calibration}

Generally the overdispersion parameters estimated by the CHT allow the model to be well calibrated, showing little signal when run on permuted data.  However permuted tests can sometimes diverge from the null, particularly when small sample sizes are used.  This may occur because by chance the permutation are unable to completely break up the signal when there aren't many samples to permute or because of innacuracy in the overdispersion estimates.  We suggest running the CHT on permuted data using the options we provide and viewing a quantile-quantile plot to ensure that the test is working properly.  If these permutations do not follow the null, the user may manually set overdispersion parameters or adjust the p-values according to the permuted distribution.


\subsection{CHT sensitivity analysis}

To evaluate the sensitivity of CHT, we compared the number eQTLs identified using CHT at 10\% FDR to the number identified using a linear model. For the linear model we divided the observed counts $x_{ij}$ by the expected number of counts, $T^*_{ij}$, which is estimated from the GC content and total read depth of each region. We then used quantile normalization to bring the distribution of counts for each individual to a standard normal distribution. We included principal components as covariates in the linear model and determined the number of principal components to include by maximizing the number of significant eQTLs.

\subsubsection{Identifying known European eQTLs in 70 Yoruba LCLs}

We downloaded a a set of 1895 eQTLs which were identified in 462 European lymphoblastoid cell lines (LCLs) by the GEUVADIS project \cite{24037378}.  We identified the subset of these eQTL SNPs that were segregating in an independent dataset of 69 Yoruba LCLs \cite{Stephens_Gilad_Pritchard_2010} with a minimum minor allele count of 2. We then applied the CHT and linear model to the RNA-seq data from the Yoruba LCLs. As target regions for each test we used the set of non-redundant Ensembl exons from the associated gene.

\subsubsection{Genome-wide QTL discovery in small sample sizes of ChIP-seq data}

We also applied the two models to a dataset of ChIP-seq data for the histone modification H3K27ac on 10 individuals, which we collected in a previous study \cite{McVicker_2013}. Principle components were not included in this analysis because of the small number of dimensions in the dataset. As test SNPs we chose SNPs that were segregating in the 10 individuals and defined the target region as a 2kb region centered on the test SNP. We only tested target regions with at least 100 filtered reads (summed across individuals).

\selectlanguage{english}
\FloatBarrier
\bibliographystyle{plain}
\bibliography{bibliography/converted_to_latex.bib%
}

\end{document}

