%\VignetteIndexEntry{GenoGAM: Genome-wide generalized additive models} %\VignettePackage{GenoGAM} %\VignetteEngine{knitr::knitr} % To compile this document % library('knitr'); rm(list=ls()); knit('GenoGAM.Rnw'); system("pdflatex GenoGAM") \documentclass[11pt]{article} \usepackage{amsmath} \newcommand{\genogam}{\textit{GenoGAM}} \newcommand{\chipseq}{ChIP-Seq} <>= library("knitr") opts_chunk$set( tidy=FALSE, dev="png", #fig.show="hide", fig.width=4, fig.height=4.5, dpi = 300, cache=TRUE, message=FALSE, warning=FALSE) @ <>= BiocStyle::latex() @ <>= library("GenoGAM") @ \author{Georg Stricker$^{1}$, Julien Gagneur$^{1}$ \\[1em] \small{$^{1}$ Technische Universit\"at M\"unchen, Department of Informatics, Garching, Germany} } \title{GenoGAM: Genome-wide generalized additive models} \begin{document} \maketitle \begin{abstract} Many genomic assays lead to noisy observations of a biological quantity of interest varying along the genome. This is the case for \chipseq{}, for which read counts reflect local protein occupancy of the ChIP-ed protein. The \genogam{} package allows statistical analysis of genome-wide data with smooth functions using generalized additive models. It provides methods for the statistical analysis of \chipseq{} data including inference of protein occupancy, and pointwise and region-wise differential analysis. Estimation of dispersion and smoothing parameters is performed by cross-validation. Scaling of generalized additive model fitting to whole chromosomes is achieved by parallelization over overlapping genomic intervals. This vignette explains the use of the package for typical \chipseq{} analysis tasks. \vspace{1em} \includegraphics[width=0.8\textwidth]{Overview} \textbf{GenoGAM version:} \Sexpr{packageVersion("GenoGAM")} \vspace{1em} \begin{center} \begin{tabular}{ | l | } \hline If you use \genogam{} in published research, please cite: \\ \\ Stricker, \emph{et al.} \textbf{Genome-wide generalized additive models} \\ \emph{bioRxiv}\\ \hline \end{tabular} \end{center} \end{abstract} <>= options(digits=3, width=80, prompt=" ", continue=" ") @ \newpage \tableofcontents \newpage \section{Standard \chipseq{} analysis} This version of \genogam{} only supports smoothing and differential analysis of \chipseq{} data. \subsection{Goal of the analysis} A small dataset is provided to illustrate the \chipseq{} functionalities. This is a subset of the data published by Thornton et al\cite{Thornton2014}, who assayed histone H3 Lysine 4 trimethylation (H3K4me3) by \chipseq{} comparing wild type yeast versus a mutant with a truncated form of Set1, the yeast H3 Lysine 4 methylase. The goal of this analysis is the identification of genomic positions that are significantly differentially methylated in the mutant compared to the wild type strain. To this end, we will build a \genogam{} model that models the logarithm of the expected ChIP-seq fragment counts $y$ as sums of smooth functions of the genomic position $x$. Specifically, we write (with simplified notations) that: \begin{equation} \label{eq:mutantmodel} \log(\operatorname{E}(y)) = f(x) + \text{genotype} \times f_\text{mutant/wt}(x) \end{equation} where genotype is 1 for data from the mutant samples, and 0 for the wild type. Here the function $f(x)$ is the reference level, i.e. the log-rate in the wild type strain. The function $f_\text{mutant/wt}(x)$ is the log-ratio of the mutant over wild-type. We will then statistically test the null hypothesis $f_\text{mutant/wt}(x) = 0$ at each position $x$. In the following we show how to build the dataset, perform the fitting of the model and perform the testing. \subsection{Registering a parallel backend} The parallel backend is registered using the \Biocpkg{BiocParallel} package. See the documentation in \Rclass{BiocParallel} for the correct use. Also note, that \Rclass{BiocParallel} is just an interface to multiple parallel packages. For example in order to use \genogam on a cluster, the \CRANpkg{BatchJobs} package might be required. The parallel backend should be registered before creating the \Rclass{GenoGAM} class, as this setup will be used throughout the analysis. <>= library(GenoGAM) ## On multicore machines by default the number of available cores - 2 are registered BiocParallel::registered()[1] @ \newpage For this small example we would like to assign less workers and activate the progress bar. Check \Biocpkg{BiocParallel} for other possible backends and more options for \Rfunction{MulticoreParam} <>= BiocParallel::register(BiocParallel::MulticoreParam(workers = 4, progressbar = TRUE)) @ If we check the current registered backend, we see that the number of workers has changed. <>= BiocParallel::registered()[1] @ \subsection{Building a \Rclass{GenoGAM} dataset} BAM files restricted to a region of chromosome XIV around the gene \emph{YNL176C} are provided in the \texttt{inst/extdata} folder of the \genogam{} package. This folder also contains a flat file describing the experimental design. We start by loading the experimental design from the tab-separated text file \texttt{experimentDesign.txt} into a data frame: <>= folder <- system.file("extdata/Set1", package='GenoGAM') expDesign <- read.delim( file.path(folder, "experimentDesign.txt") ) expDesign @ Each row of the experiment design corresponds to the alignment files in BAM format of one ChIP sample. In case of multiplexed sequencing, the BAM files must have been demultiplexed. The experiment design have named columns. Three column names have a fixed meaning for \genogam{} and must be provided: \Robject{ID}, \Robject{file}, and \Robject{paired}. The field \Robject{ID} stores a unique identifier for each alignment file. It is recommended to use short and easy to understand identifiers because they are subsequently used for labelling data and plots. The field \Robject{file} stores the BAM file name. The field \Robject{paired} values \Robject{TRUE} for paired-end sequencing data, and \Robject{FALSE} for single-end sequencing data. Further named columns can be added at wish without naming and data type constraints. Here the important one is the \Robject{genotype} column. Note that it is an indicator variable (i.e. valuing 0 or 1). It will allow us modeling the differential occupancy later on. We will now count sequencing fragment centers per genomic position and sample and store these counts into a \Rclass{GenoGAMDataSet}. \genogam{} reduces \chipseq{} data to fragment center counts rather than full base coverage so that each fragment is counted only once. This reduces artificial correlation between adjacent nucleotides. For single-end libraries, the fragment center is estimated by shifting the read end position by a constant (Details in the help on the constructor function \Rfunction{GenoGAMDataSet()}). <>= bpk <- 20 chunkSize <- 1000 overhangSize <- 15*bpk ## build the GenoGAMDataSet ggd <- GenoGAMDataSet( expDesign, directory = folder, chunkSize = chunkSize, overhangSize = overhangSize, design = ~ s(x) + s(x, by = genotype) ) ## restricts the GenoGAM dataset to the positions of interest ## (this step is only required for running this small example) ggd <- subset(ggd, seqnames == "chrXIV" & pos >= 305000 & pos <= 308000) @ A \Rclass{GenoGAMDataSet} stores this count data into a structure that index genomic positions over \textit{tiles}, defined by \Robject{chunkSize} and \Robject{overhangSize}. A bit of background is required to understand these parameters. The smoothing in \genogam{} is based on splines, which are piecewise polynomials. The \textit{knots} are the positions where the polynomials connect. In our experience, one knot every 20 to 50 bp is required for enough resolution of the smooth fits in typical applications. The fitting of generalized additive models involves steps demanding a number of operations proportional to the square of the number of knots, preventing fits along whole chromosomes. To make the fitting of GAMs genome-wide, \genogam{} performs fitting on small overlapping intervals (\textit{tiles}), and join the fit at the midpoint of the overlap of consecutive tiles. The parameters \Robject{chunkSize} and \Robject{overhangSize} defines the tiles, where the chunk is the core part of a tile that does not overlap other tiles, and the overhangs are the two overlapping parts. Overhangs of about 10 times the knot spacing gives reasonable results. The \Robject{design} parameter is explained in the next section. Finally, the last line of code calls the function \Rfunction{subset()} to restrict the \Rclass{GenoGAMDataset} to the positions of interest. This line is necessary for running this small example but would not be present in a standard genome-wide run of \genogam{}. \subsection{Modeling with smooth functions: the \Robject{design} parameters} \genogam{} models the logarithm of the rate of the count data as sums of smooth functions of the genomic position, denoted \Robject{x}. The \Robject{design} parameter is an R formula which allows encoding how the smooth functions depend on the experimental design. \genogam follows formula convention of the R package \Rpackage{mgcv}. A smooth function is denoted \Robject{s()}. For now, \genogam{} only supports smooth function that are cubic splines of the genomic position \Robject{x}. The \Robject{by} variable allows selecting to which samples the smooth contributes to (see also the documentation of \Robject{gam.models} in the \Rpackage{mgcv}). For now, \genogam{} only allows \Robject{by} variables to (value 0 or 1). Here by setting \Robject{'s(x, by=genotype)'} we encode the term "$ \text{genotype} \times f_\text{mutant/wt}(x)$" in Equation ~\ref{eq:mutantmodel}. \textbf{Note:} As for other generalized additive models packages (\Rpackage{mgcv}, \Rpackage{gam}), \genogam{} use the natural logarithm as link function. This is different than other packages of the bioinformatics files such as \Rpackage{DESeq2} which works in base 2 logarithm. \subsection{Size factor estimation} Sequencing libraries typically vary in sequencing depth. Such variations is controlled for in \genogam{} by adding a sample-specific constant to the right term of Equation ~\ref{eq:mutantmodel}. The estimation of these constants is performed by the function computeSizeFactor() as follows: <>= ggd <- computeSizeFactors(ggd) sizeFactors(ggd) @ \textbf{Note:} The size factors in \genogam are in the natural logarithm scale. \subsection{Model fitting} A \genogam{} model requires two further parameters to be fitted: the regularization parameter, $\lambda$, and the dispersion parameter $\theta$. The regularization parameters $\lambda$ controls the amount of smoothing. The larger $\lambda$ is, the smoother the smooth functions are. The dispersion parameter $\theta$ controls how much the observed counts deviate from their expected value modeled by Equation ~\ref{eq:mutantmodel}. The dispersion captures biological and technical variation which one typically sees across replicate samples, but also errors of the model. In \genogam, the dispersion is modeled by assuming the counts to follow a negative binomial distribution with mean $\mu=\operatorname{E}(y)$ whose logarithm is modeled by Equation ~\ref{eq:mutantmodel} and with variance $v = \mu + \mu^2/\theta$. If not provided, the parameters $\lambda$ and $\theta$ are obtained by cross-validation. This step is a bit time-consuming. For sake of going through this example quickly, we provide the values manually: <>= ## fit model without parameters estimation fit <- genogam(ggd, lambda = 40954.1, family = mgcv::nb(theta = 6.927986), bpknots = bpk ) fit @ \textbf{Remark on parameter estimation:} To estimate the parameters $\lambda$ and $\theta$ by cross-validation, call \Rfunction{genogam()} without setting their values. This will perform 10 fold cross-validation on each tile with initial parameter values and iterate until convergence, often for about 50 iterations. We recommend to do it for 20 to 40 different regions representative of your data (of 1.5kb each). This means that estimation of the parameters will require the equivalent of a \genogam{} fit with fixed $\lambda$ and $\theta$ on 30 Mb (1.5kb x10x40x50). For a genome like yeast (12Mb) the cross-validation thus can take more time than a genome-wide fit. <>= fit_CV <- genogam(ggd,bpknots = bpk) @ \textbf{Remark on parallel computing:} \genogam{} run parallel computations on multicore architecture (using the \Rpackage{BiocParallel} package). Computing time reduces almost linearly with the number of cores of the machine. \subsection{Plotting results} Count data and fits for a region of interest can be extracted using the function \Rfunction{view()}. Following the \Rpackage{mgcv} and \Rpackage{gam} convention the names of the fit for the smooth function defined by the \Robject{by} variable follow the pattern \Robject{s(x):\{by-variable\}}. Here, the smooth function of interest $f_\text{mutant/wt}(x)$ is thus named \Robject{s(x):genotype}. <>= # extract count data into a data frame df_data <- view(ggd) head(df_data) # extract fit into a data frame df_fit <- view(fit) head(df_fit) @ The code below then plots the counts and the fitted smooth log-ratio of mutant over wild type together with its 95\% confidence band. In the count data, the peak of methylation in the two replicates of the wild type (first two panels) in the region 306,500-307,000 seems attenuated in the two replicates of the mutant (3rd and 4th panel). There are relatively more counts for the mutant in the region 305,500-306,500. The \genogam{} fit of the log-ratio (last panel, confidence band dotted) indicates that that these differences are significant. This redistribution of the methylation mark from the promoter (wild type) into the gene body (mutant) was reported by the authors of the study ~\cite{Thornton2014}. <>= ## plot function for the count data dataplot <- function(df, col, ...){ x <- df[["pos"]] y <- df[[col]] plot(0, type='n', xlim=range(x), ylab=col, ... ) points(x, y, pch=19, col="#00000015", cex=0.5) } # plot function for a fit with confidence band fitplot <- function(df, smooth, ...){ x <- df[["pos"]] y <- df[[smooth]] se.y <- df[[paste0("se.",smooth)]] plot(0, type='n', xlim=range(x), ylim=range(c(y - 1.96*se.y, y + 1.96*se.y)), ylab=smooth, ... ) lines(x, y) lines(x, y+1.96*se.y, lty='dotted') lines(x, y-1.96*se.y, lty='dotted') abline(h=0) } ## plot par(mfrow=c(5,1)) par(mar=c(1,4,1,1)) for(id in expDesign$ID) dataplot(df_data, id, xlab="", main="", ylim=c(0,6)) par(mar=c(2,4,1,1)) fitplot(df_fit, 's(x):genotype', xlab="", main="") @ \subsection{Statistical testing} We test for each smooth and at each position $x$ the null hypothesis that it values 0 by a call to \Rfunction{computeSignificance()}. False discovery rate can be computed using the Benjamini-Hochberg procedure with the R function \Rfunction{p.adjust()}: <>= fit <- computeSignificance(fit) df_fit <- view(fit) df_fit[["fdr.s(x):genotype"]] <- p.adjust(df_fit[["pvalue.s(x):genotype"]], method="BH") head(df_fit) @ \section{Other functionalities} Other functionalities demonstrated in the manuscript (peak calling and testing, methylation data, see BioRxiv: \url{http://dx.doi.org/10.1101/047464}) will be gradually integrated into the package. Interested users should check the developer version of \Biocpkg{GenoGAM} for updates. \section{Acknowledgments} We thank Alexander Engelhardt, Herv\'e Pag\`es, and Martin Morgan for input in the development of \genogam{}. \section{Session Info} <>= toLatex(sessionInfo()) @ <>= options(prompt="> ", continue="+ ") @ \bibliography{bibliog} \end{document}