%\VignetteIndexEntry{coMET users guide} %\VignetteDepends{coMET} %\VignetteKeywords{Software, DifferentialMethylation, Visualization, Sequencing, Genetics, FunctionalGenomics, Microarray, MethylationArray, MethylSeq, ChIPSeq, DNASeq, RIBOSeq, RNASeq, ExomeSeq, DNAMethylation, GenomeWideAssociation } %\VignettePackage{coMET} %\VignetteEngine{knitr::knitr} \documentclass[11pt]{article} % A bunch of styles and package requirements for the Bioconductor vignette branding %<>= <>= #library("BiocStyle") BiocStyle::latex() @ <>= library(knitr) opts_chunk$set(fig.path='figure/minimal-', fig.align='center', fig.show='hold') options(replace.assign=TRUE,width=90) @ \RequirePackage[utf8]{inputenc} \RequirePackage{hyperref} \RequirePackage{url} \RequirePackage{listings} %\RequirePackage[numbers]{natbib} %\bibliographystyle{plainnat} %\bibpunct{(}{)}{;}{a}{,}{,} % \RequirePackage[text={7.2in,9in},centering]{geometry} %\setkeys{Gin}{width=0.95\textwidth} \RequirePackage{longtable} \RequirePackage{graphicx} \newcommand{\code}[1]{{\texttt{#1}}} \newcommand{\term}[1]{{\emph{#1}}} \newcommand{\Rmethod}[1]{{\textit{#1}}} \newcommand{\Rfunarg}[1]{{\textit{#1}}} \newcommand{\scscst}{\scriptscriptstyle} \newcommand{\scst}{\scriptstyle} \newcommand{\mgg}[0]{\Rpackage{coMET} } \newcommand{\Reference}[1]{{\texttt{#1}}} \newcommand{\link}[1]{{#1}} \newcommand{\RR}[0]{{\texttt{R}}} \newcommand\BiocStyle{\Rpackage{BiocStyle}} \newcommand\knitr{\Rpackage{knitr}} %\footnote{thomas.hardiman@kcl.ac.uk}, \footnote{idil.yet@kcl.ac.uk}, \footnote{peichien.tsai@kcl.ac.uk}, \footnote{jordana.bell@kcl.ac.uk} \bioctitle[The coMET User Guide]{The coMET User Guide} \author{Tiphaine C. Martin \footnote{\email{tiphaine.martin@mssm.edu}}, Tom Hardiman, Idil Yet, Pei-Chien Tsai, Jordana T. Bell} \date{Edited: December 2018; Compiled: \today} \begin{document} \maketitle \section{Citation} <>= citation(package='coMET') @ \clearpage \tableofcontents \clearpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{abstract} The \Biocpkg{coMET} package is a web-based plotting tool and R-based package to visualize omic-WAS results in a genomic region of interest, such as EWAS (epigenome-wide association scan). \Biocpkg{coMET} provides a plot of the EWAS association signal and visualisation of the methylation correlation between CpG sites (co-methylation). The \Biocpkg{coMET} package also provides the option to annotate the region using functional genomic information, including both user-defined features and pre-selected features based on the Encode project. The plot can be customized with different parameters, such as plot labels, colours, symbols, heatmap colour scheme, significance thresholds, and including reference CpG sites. Finally, the tool can also be applied to display the correlation patterns of other genomic data in any species, e.g. gene expression array data. \Biocpkg{coMET} generates a multi-panel plot to visualize EWAS results, co-methylation patterns, and annotation tracks in a genomic region of interest. A \Biocpkg{coMET} figure (cf. Fig. 1) includes three components: \begin{enumerate} \item the upper plot shows the strength and extent of EWAS association signal; \item the middle panel provides customized annotation tracks; \item the lower panel shows the correlation between selected CpG sites in the genomic region. \end{enumerate} The structure of the plots builds on \CRANpkg{snp.plotter}\cite{Luna2007}, with extensions to incorporate genomic annotation tracks and customized functions. \Biocpkg{coMET} produces plots in PDF and Encapsulated Postscript (EPS) format. The current version of \Biocpkg{coMET} can visualise EWAS results and annotations from a genomic region up to an entire chromosome in the upper and middle panels of the coMET plot. However, the lower panel (co-methylation) is restricted to visualising a maximum of 120 single-CpG or region-based datapoints. This limitation is due to limitations in the size of a standard A4 plot, and may be updated in the near future. However, the user can use the function \Rfunction{comet.list} to extracts all significant correlations beyond a given threshold in the dataset from either a genomic region or from an entire chromosome if required. \end{abstract} \clearpage \section{Usage} \Biocpkg{coMET} requires the installation of R, the statistical computing software, freely available for Linux, Windows, or MacOS. \Biocpkg{coMET} can be downloaded from bioconductor. Packages can be installed using the install.packages command in R. The \Biocpkg{coMET} R package includes two major functions \Rfunction{comet.web} and \Rfunction{comet} to visualise omci-WAS results. \begin{itemize} \item The function \Rfunction{comet.web} generates output plot with the same settings of genomic annotation tracks as that of the webservice (\url{http://epigen.kcl.ac.uk/comet} or direcly \url{http://comet.epigen.kcl.ac.uk:3838/coMET/}). \item The function \Rfunction{comet} generates output plots with the customized annotation tracks defined by user. \end{itemize} <>= if (!requireNamespace("BiocManager", quietly=TRUE)) install.packages("BiocManager") BiocManager::install("coMET") @ \Biocpkg{coMET} uses the packages called \CRANpkg{psych}, \CRANpkg{corrplot} and \CRANpkg{colortools}, which are not available from bioconductor. This must be installed before the installation of \Biocpkg{coMET} <>= install.packages("psych") install.packages("corrplot") install.packages("colortools") @ \Biocpkg{coMET} has a development version on gitHub, go to the section "Install the development version of \Biocpkg{coMET} from \Bioconductor{}". You can install also on the version R 3.2.2 via the master version of package on gitHub. The same steps must be followed as described in the section "Install the development version of \Biocpkg{coMET} from \Bioconductor{}". After downloading from \Bioconductor{} or gitHUB, and installing on your computer, \Biocpkg{coMET} can be loaded into a R session using this command: <>= require("hash") require("grid") require("grDevices") require("biomaRt") require("Gviz") require("ggbio") require("rtracklayer") require("GenomicRanges") require("colortools") require("gridExtra") require("psych") rdir <- system.file("R", package="coMET",mustWork=TRUE) source(file.path(rdir, "AnalyseFile.R")) source(file.path(rdir, "BiofeatureGraphics.R")) source(file.path(rdir, "comet.R")) source(file.path(rdir, "cometWeb.R")) source(file.path(rdir, "DrawPlot.R")) source(file.path(rdir, "GeneralMethodComet.R")) @ <>= library("coMET") @ The configuration file specifies the options for the coMET plot. Example configuration and input files are also provided on \url{http://epigen.kcl.ac.uk/comet}. Information about the package can viewed from within R using this command: <>= ?comet ?comet.web ?comet.list @ \subsection{Install the development version of coMET from Bioconductor} To install \Biocpkg{coMET} from the development version of Bioconductor, the user must install the appropriate R version. See \url{http://www.bioconductor.org/developers/how-to/useDevel/} for more details. Following this installation, use the standard Bioconductor command: <>= if (!requireNamespace("BiocManager", quietly=TRUE)) install.packages("BiocManager") BiocManager::install(version = "devel") BiocManager::install("coMET") @ \subsection{Install the version of coMET from gitHub} Another way to install \Biocpkg{coMET} is to download the master package from gitHUB \url{https://github.com/TiphaineCMartin/coMET} or the devel package \url{https://github.com/TiphaineCMartin/coMET/tree/devel}. Once downloaded use command line: <>= install.packages("YourPath/coMET_YourVersion.tar.gz",repos=NULL,type="source") ##This is an example install.packages("YourPath/coMET_0.99.9.tar.gz",repos=NULL,type="source") @ \clearpage \section{Functions in coMET} Currently, there are 3 main functions: \begin{enumerate} \item \Rfunction{comet.web} is the pre-customized function that allows us to visualise quickly EWAS (or other omic-WAS) results, annotation tracks, and correlations between features. This version is installed in the Shiny web-service. Currently, it is formated only to visualise human data. \item \Rfunction{comet} is the generic function that allows us to visualise quickly EWAS results, annotation tracks, and correlations between features. Users can visualise more personalised annotation tracks and give multiple extra EWAS/omic-WAS results to plot. \item \Rfunction{comet.list} is an additional function that allows us to extract the values of correlations, the pvalues, and estimates and confidence intervals for all datapoints that surpass a particular threshold. \end{enumerate} The functions can read the data input files, but it is also possible to use data frames within R for all data input except for the configuration file. The latter can be achieved with the two functions \Rfunction{comet} and \Rfunction{comet.list}. The structure of the data frames (number of columns, type, format) follows the same rules as for the data input files (cf. section "File formats"). \clearpage \section{File formats} There are five types of files that can be given by the user to produce the plot: \begin{enumerate} \item Info file is defined in the option \Robject{mydata.file}. \warning{This is mandatory and has to be in tabular format with a header}. \item Correlation file is defined in the option \Robject{cormatrix.file}. \warning{This is mandatory and has to be in tabular format with a header}. \item Extra info files are defined in the option \Robject{mydata.file.large}. \warning{This is optional, and if provided has to be in tabular format with a header}. \item Annotation info file is defined in the option \Robject{biofeat.user.file}. This option exists only in the function \Rfunction{comet.web} and the user should inform also the format to visualise this data with the options \Robject{biofeat.user.type} and \Robject{biofeat.user.type.plot}. \item Configuration file contains the values of these options instead of defining these by command line. \warning{Each line in the file is one option. The name of the option is in capital letters and is separated by its value by "="}. If there are multiple values such as for the option \Robject{list.tracks} or the options for additional data, you need to separated them by a "comma". \end{enumerate} \subsection{Format of the info file (for option: \Robject{mydata.file}, mandatory)} \warning{This file is mandatory and has to be in tabular format with a header. The name of features has to start by a letter}. Info files can be a list of CpG sites with/without Beta value (for example DNA methylation level) or direction sign. If it is a site file then it is mandatory to have the 4 columns as shown below with headers in the same order. Beta can be the 5th column(optional) and can be either a numeric value (positive or negative values) or only direction sign ("+", "-"). The number of columns and their types are defined by the option \Robject{mydata.format}. <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) infofile <- file.path(extdata, "cyp1b1_infofile.txt") data_info <-read.csv(infofile, header = TRUE, sep = "\t", quote = "") head(data_info) @ Alternatively, the info file can be region-based and if so, the region-based info file must have the 5 columns (see below) with headers in this order. The beta or direction can be included in the 6th column (optional). <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) infoexp <- file.path(extdata, "cyp1b1_infofile_exprGene_region.txt") data_infoexp <-read.csv(infoexp, header = TRUE, sep = "\t", quote = "") head(data_infoexp) @ In summary, there are 4 possible formats for the info file: \begin{enumerate} \item \textbf{\emph{site}}: 4 columns with a header: \begin{enumerate} \item Name of omic feature \item Name of chromosome \item Position of omic feature \item P-value of omic feature \end{enumerate} \item \textbf{\emph{region}}: 5 columns with a header: \begin{enumerate} \item Name of omic feature \item Name of chromosome \item Start position of omic feature \item End position of omic feature \item P-value of omic feature \end{enumerate} \item \textbf{\emph{site\_asso}}: 5 columns with a header: \begin{enumerate} \item Name of omic feature \item Name of chromosome \item Position of omic feature \item P-value of omic feature \item Direction of association related to this omic feature. This can be the sign or an actual value of association effect size. \end{enumerate} \item \textbf{\emph{region\_asso}}: 6 columns with a header: \begin{enumerate} \item Name of omic feature \item Name of chromosome \item Start position of omic feature \item End position of omic feature \item P-value of omic feature \item Direction of association related to this omic feature. This can be the sign or an actual value of association effect size. \end{enumerate} \end{enumerate} \subsection{Format of correlation matrix (for option: \Robject{cormatrix.file}, mandatory)} \warning{This file is mandatory and has to be in tabular format with an header}. The data file used for the correlation matrix is described in the option \Robject{cormatrix.file}. This tab-delimited file can take 3 formats described in the option \Robject{cormatrix.format}: \begin{enumerate} \item \textbf{\emph{cormatrix}}: pre-computed correlation matrix provided by the user; Dimension of matrix : CpG\_number X CpG\_number. Need to put the CpG sites/regions in the ascending order of positions and to have a header with the name of CpG sites/regions; \item \textbf{\emph{raw}}: Raw data format. Correlations of these can be computed by one of 3 methods Spearman, Pearson, Kendall (option \Robject{cormatrix.method}). Dimension of matrix : sample\_size X CpG\_number. Need to have a header with the name of CpG sites/regions ; \item \textbf{\emph{raw\_rev}}: Raw data format. Correlations of these can be computed by one of 3 methods Spearman, Pearson, Kendall (option \Robject{cormatrix.method}). Dimension of matrix : CpG\_number X sample\_size. Need to have the row names of CpG sites/regions and a header with the name of samples ; \end{enumerate} <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) corfile <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") data_cor <-read.csv(corfile, header = TRUE, sep = "\t", quote = "") data_cor[1:6,1:6] @ \subsection{Format of extra info file (for option: \Robject{mydata.large.file})} \warning{This file is optional file and if provided has to be in tabular format with an header. The name of features has to start by a letter}. The extra info files can be described in the option \Robject{mydata.large.file} and their format in \Robject{mydata.large.format}. More than one extra info file can be used, each should be separated by a comma. This can be another type of info file (e.g expression or replication data) and should follow the same rules as the standard info file. \subsection{Format of annotation file (for option \Robject{biofeat.user.file})} The file is defined in the option \Robject{biofeat.user.file} and the format of file is the format accepted by \Biocpkg{Gviz} (BED, GTF, and GFF3). \subsection{Option of config.file} \warning{Each line in the file is one option. The name of the option is in lowercase letters and is separated by its value by "=" without space. If there are multiple values such as for the option \Robject{list.tracks} or options for additional data, these need to be separated them by a "comma" withou space}. If you would like to make your own changes to the plot you can download the configuration file, make changes to it, and upload it into R as shown in the example below. The important options of a \Biocpkg{coMET} figure include three components: \begin{enumerate} \item The \Robject{upper plot} shows the strength and extent of EWAS association signal on a regional Manhattan plot. \begin{itemize} \item \Robject{pval.threshold}: Significance threshold to be displayed as a red dashed line \item \Robject{pval.threshold2}: Another Significance threshold (optional) \item \Robject{disp.pvalueplot}: Value can be TRUE or FALSE. Used to either display or hide Manhattan plot. \item \Robject{disp.beta.association}: Value can be TRUE or FALSE. Used to show the effect size. \item \Robject{disp.association}: This logical option works only if \Robject{mydata.file} contains the effect direction (\Robject{mydata.format}=\textbf{\emph{site\_asso}} or \textbf{\emph{region\_asso}}). The value can be TRUE or FALSE: if FALSE (default), for each point of data in the p-value plot, the colour of symbol is the colour of co-methylation pattern between the point and the reference site; if TRUE, the effect direction is shown. If the association is positive, the colour is the one defined with the option \Rfunction{color.list}. On the other hand, if the association is negative, the colour is the inverse to that selected. \item \Rfunction{disp.region} : This logical option works only if the option \Rfunction{mydata.file} contains regions (\Rfunction{mydata.format}=\textbf{\emph{region}} or \textbf{\emph{region\_asso}}). The value can be TRUE or FALSE (default). If TRUE, the genomic element will be shown as a continuous line with the colour of the element, in addition to the symbol at the center of the region. If FALSE, only the symbol is shown. \end{itemize} \item The \textbf{\emph{middle panel}} provides customized annotation tracks; \begin{itemize} \item \Robject{list.tracks} (for \Rfunction{comet.web} function): List of annotation tracks to be visualised. Tracks currently available: geneENSEMBL, CGI, ChromHMM, DNAse, RegENSEMBL, SNP, transcriptENSEMBL, SNPstoma, SNPstru, SNPstrustoma, GAD, ClinVar, GeneReviews, GWAS, ClinVarCNV, GCcontent, genesUCSC, xenogenesUCSC, metQTL, eQTL, BindingMotifsBiomart, chromHMM\_RoadMap, miRNATargetRegionsBiomart, OtherRegulatoryRegionsBiomart, RegulatoryEvidenceBiomart and segmentalDupsUCSC. The elements are separated by a comma. \item \Robject{tracks.gviz} (for \Rfunction{comet} function): For each option, it is possible to give a list of annotation tracks that is created by the \Biocpkg{Gviz} bioconductor packages. \warning{It is noted that the new version of coMET does not support more the visualisation using \Robject{tracks.ggbio} and \Robject{tracks.trackviewer} from \Biocpkg{GGBio} and \Biocpkg{TrackViewer} because The integration of plots from ggbio and trackviewer can be sometimes not really perfect. So only now, it is possible to create plots from \Biocpkg{Gviz} and use \Robject{tracks.gviz}} \end{itemize} \item The \textbf{\emph{lower panel}} shows the correlation between selected CpG sites in the genomic region (heatmap). \begin{itemize} \item \Robject{cormatrix.format} : Format of the input file \textbf{\emph{cormatrix.file}}: either raw data (option RAW if CpG sites are by column and samples by row or option RAW\_REV if CpG site are by row and samples by column) or correlation matrix (option CORMATRIX) \item \Robject{cormatrix.method} : If raw data are provided it will be necessary to produce the correlation matrix using one of 3 methods (spearman, pearson and kendall). \item \Robject{cormatrix.color.scheme} : There are 5 colour schemes (heat, bluewhitered, cm, topo, gray, bluetored) \item \Robject{disp.cormatrixmap} : logical option TRUE or FALSE. TRUE (default), if FALSE correlation matrix is not shown \item \Robject{cormatrix.conf.level} : Alpha level for the confidence interval. Default value= 0.05. CI will be the alpha/2 lower and upper values. \item \Robject{cormatrix.sig.level} : Significant level to visualise the correlation. If the correlation has a pvalue under the significant level, the correlation will be colored in "goshwhite", else the color is related to the correlation level and the color scheme choosen.Default value =1. \item \Robject{cormatrix.adjust} : Indicates which adjustment for multiple tests should be used. "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none".Default value="none". \end{itemize} \end{enumerate} <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) configfile <- file.path(extdata, "config_cyp1b1_zoom_4webserver_Grch38.txt") data_config <-read.csv(configfile, quote = "", sep="\t", header=FALSE) data_config @ \clearpage \section{Creating a plot like the webservice: comet.web} \subsection{\Biocpkg{coMET} plot: usage and plot like in the webservice} The user can create a \Biocpkg{coMET} plot via the \Biocpkg{coMET} website (\url{http://epigen.kcl.ac.uk/comet}). It is possible to reproduce the web service plotting defaults by using the function \Rfunction{comet.web}, for example see Figure \ref{fig:cometweb_simple}. <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) myinfofile <- file.path(extdata, "cyp1b1_infofile_Grch38.txt") myexpressfile <- file.path(extdata, "cyp1b1_infofile_exprGene_region_Grch38.txt") mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") configfile <- file.path(extdata, "config_cyp1b1_zoom_4webserver_Grch38.txt") comet.web(config.file=configfile, mydata.file=myinfofile, cormatrix.file=mycorrelation ,mydata.large.file=myexpressfile, print.image=FALSE,verbose=FALSE) @ \begin{figure} <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) myinfofile <- file.path(extdata, "cyp1b1_infofile_Grch38.txt") myexpressfile <- file.path(extdata, "cyp1b1_infofile_exprGene_region_Grch38.txt") mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") configfile <- file.path(extdata, "config_cyp1b1_zoom_4webserver_Grch38.txt") comet.web(config.file=configfile, mydata.file=myinfofile, cormatrix.file=mycorrelation ,mydata.large.file=myexpressfile, print.image=FALSE,verbose=FALSE) @ %\incfig{fig:cometweb_simple}{\textwidth}{Plot with comet.web function.} \caption{Plot with \Rfunction{comet.web} function.\label{fig:cometweb_simple}} \end{figure} \subsection{Hidden values of \Rfunction{comet.web} function} Hidden values of \Rfunction{comet.web} function are shown in the section. If these values do not correspond to what you want to visualise, you need to use the function \Rfunction{comet}, as a more generic option. \begin{longtable}{|c|c|} \hline \multicolumn{1}{|c|}{Option} & \multicolumn{1}{c|}{Value} \\ \hline \endfirsthead \multicolumn{2}{c}% {\tablename\ \thetable\ -- continued from previous page} \\ \hline \multicolumn{1}{|c|}{Option} & \multicolumn{1}{c|}{Value} \\ \hline \endhead \hline \multicolumn{2}{|r|}{{Continued on next page}} \\ \hline \endfoot \hline \hline \endlastfoot mydata.type & FILE \\ mydata.large.type & LISTFILE \\ cormatrix.type & LISTFILE \\ disp.cormatrixmap & TRUE\\ disp.pvalueplot & TRUE\\ disp.mydata.names & TRUE\\ disp.connecting.lines & TRUE\\ disp.mydata & TRUE\\ disp.type & symbol \\ biofeat.user.type.plot & histogram \\ tracks.gviz & NULL\\ tracks.ggbio & NULL\\ tracks.trackviewer & NULL\\ biofeat.user.file & NULL\\ palette.file & NULL\\ disp.color.bar & TRUE\\ disp.phys.dist & TRUE\\ disp.legend & TRUE\\ disp.marker.lines & TRUE\\ disp.mult.lab.X & FALSE\\ connecting.lines.factor & 1.5\\ connecting.lines.adj & 0.01\\ connecting.lines.vert.adj & -1\\ connecting.lines.flex & 0\\ color.list & red \\ font.factor & NULL\\ dataset.gene & hsapiens\_gene\_ensembl\\ DATASET.SNP & hsapiens\_snp\\ VERSION.DBSNP & snp142Common\\ DATASET.SNP.STOMA & hsapiens\_snp\_som\\ DATASET.REGULATION & hsapiens\_feature\_set\\ DATASET.STRU & hsapiens\_structvar\\ DATASET.STRU.STOMA & hsapiens\_structvar\_som\\ BROWSER.SESSION & UCSC\\ \end{longtable} \clearpage \section{Creating a plot with the generic function: \Rfunction{comet}} It is possible to create the annotation tracks by \Biocpkg{Gviz}, \Biocpkg{trackviewer} or \Biocpkg{ggbio}, for example see Figure \ref{fig:cometPlotfile}. Currently, the \Biocpkg{Gviz} option for annotation tracks, in combination with the heatmap of correlation values between genomic elements, provides the most informative and easy approach to visualize graphics. \subsection{\Biocpkg{coMET} plot: pvalue plot, annotation tracks, and correlation matrix} \subsubsection{Input from data files} In this figure \ref{fig:cometPlotfile}, we create different tracks outside to \Biocpkg{coMET} with \Biocpkg{Gviz}. The list of annotation tracks and different files are given to the function \Rfunction{coMET}. <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) configfile <- file.path(extdata, "config_cyp1b1_zoom_4comet.txt") myinfofile <- file.path(extdata, "cyp1b1_infofile.txt") myexpressfile <- file.path(extdata, "cyp1b1_infofile_exprGene_region.txt") mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") chrom <- "chr2" start <- 38290160 end <- 38303219 gen <- "hg19" strand <- "*" BROWSER.SESSION="UCSC" mySession <- browserSession(BROWSER.SESSION) genome(mySession) <- gen genetrack <-genes_ENSEMBL(gen,chrom,start,end,showId=TRUE) snptrack <- snpBiomart_ENSEMBL(gen,chrom, start, end, dataset="hsapiens_snp_som",showId=FALSE) cpgIstrack <- cpgIslands_UCSC(gen,chrom,start,end) prombedFilePath <- file.path(extdata, "/RoadMap/regions_prom_E063.bed") promRMtrackE063<- DNaseI_RoadMap(gen,chrom,start, end, prombedFilePath, featureDisplay='promotor', type_stacking="squish") bedFilePath <- file.path(extdata, "RoadMap/E063_15_coreMarks_mnemonics.bed") chromHMM_RoadMapAllE063 <- chromHMM_RoadMap(gen,chrom,start, end, bedFilePath, featureDisplay = "all", colorcase='roadmap15' ) listgviz <- list(genetrack,snptrack,cpgIstrack,promRMtrackE063,chromHMM_RoadMapAllE063) comet(config.file=configfile, mydata.file=myinfofile, mydata.type="file", cormatrix.file=mycorrelation, cormatrix.type="listfile", mydata.large.file=myexpressfile, mydata.large.type="listfile", tracks.gviz=listgviz, verbose=FALSE, print.image=FALSE) @ \begin{figure} <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) configfile <- file.path(extdata, "config_cyp1b1_zoom_4comet.txt") myinfofile <- file.path(extdata, "cyp1b1_infofile.txt") myexpressfile <- file.path(extdata, "cyp1b1_infofile_exprGene_region.txt") mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") chrom <- "chr2" start <- 38290160 end <- 38303219 gen <- "hg19" strand <- "*" data(geneENSEMBLtrack) data(snpBiomarttrack) data(cpgIslandtrack) data(promRMtrackE063) data(chromHMM_RoadMapAllE063) listgviz <- list(genetrack,snptrack,cpgIstrack,promRMtrackE063,chromHMM_RoadMapAllE063) comet(config.file=configfile, mydata.file=myinfofile, mydata.type="file", cormatrix.file=mycorrelation, cormatrix.type="listfile", mydata.large.file=myexpressfile, mydata.large.type="listfile", tracks.gviz=listgviz, verbose=FALSE, print.image=FALSE) @ \caption{Plot with \Rfunction{comet} function from files.\label{fig:cometPlotfile}} \end{figure} \clearpage \subsubsection{\Biocpkg{coMET} plot using input from a data frame} In this figure \ref{fig:cometPlotMatrix}, we visualize the same data as in figure \ref{fig:cometPlotfile}, but the data is in data frame format and not read in from an input file. In addition, if the user would like to visualise only the correlations between CpG sites with P-value less than or equal to 0.05 in the upper plot, this option can be included. The correlations with a P-value greater than 0.05 can have the colour "goshwhite" whereas the other correlations will be displayed using a colour related to the correlation level. Conversely, in the P-value plot (upper plot), the points of each omic feature have their colours related to their correlations with the reference omic feature without taking into account the P-value associated with the correlation matrix. Eventually, we increase the size of font using the option \Robject{fontsize.gviz} <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) configfile <- file.path(extdata, "config_cyp1b1_zoom_4comet.txt") myinfofile <- file.path(extdata, "cyp1b1_infofile.txt") myexpressfile <- file.path(extdata, "cyp1b1_infofile_exprGene_region.txt") mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") chrom <- "chr2" start <- 38290160 end <- 38303219 gen <- "hg19" strand <- "*" BROWSER.SESSION="UCSC" mySession <- browserSession(BROWSER.SESSION) genome(mySession) <- gen genetrack <-genes_ENSEMBL(gen,chrom,start,end,showId=TRUE) snptrack <- snpBiomart_ENSEMBL(gen,chrom, start, end, dataset="hsapiens_snp_som",showId=FALSE) #Data no more available in UCSC (from September 2015) iscatrack <-ISCA_UCSC(gen,chrom,start,end,mySession, table="iscaPathogenic") listgviz <- list(genetrack,snptrack,iscatrack) matrix.dnamethylation <- read.delim(myinfofile, header=TRUE, sep="\t", as.is=TRUE, blank.lines.skip = TRUE, fill=TRUE) matrix.expression <- read.delim(myexpressfile, header=TRUE, sep="\t", as.is=TRUE, blank.lines.skip = TRUE, fill=TRUE) cormatrix.data.raw <- read.delim(mycorrelation, sep="\t", header=TRUE, as.is=TRUE, blank.lines.skip = TRUE, fill=TRUE) listmatrix.expression <- list(matrix.expression) listcormatrix.data.raw <- list(cormatrix.data.raw) comet(config.file=configfile, mydata.file=matrix.dnamethylation, mydata.type="dataframe",cormatrix.file=listcormatrix.data.raw, cormatrix.type="listdataframe",cormatrix.sig.level=0.05, cormatrix.conf.level=0.05, cormatrix.adjust="BH", mydata.large.file=listmatrix.expression, mydata.large.type="listdataframe", fontsize.gviz =12, tracks.gviz=listgviz,verbose=FALSE, print.image=FALSE) @ \begin{figure} <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) configfile <- file.path(extdata, "config_cyp1b1_zoom_4comet.txt") myinfofile <- file.path(extdata, "cyp1b1_infofile.txt") myexpressfile <- file.path(extdata, "cyp1b1_infofile_exprGene_region.txt") mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") #configfile <- "../inst/extdata/config_cyp1b1_zoom_4comet.txt" chrom <- "chr2" start <- 38290160 end <- 38303219 gen <- "hg19" strand <- "*" data(geneENSEMBLtrack) data(snpBiomarttrack) data(ISCAtrack) listgviz <- list(genetrack,snptrack,iscatrack) matrix.dnamethylation <- read.delim(myinfofile, header=TRUE, sep="\t", as.is=TRUE, blank.lines.skip = TRUE, fill=TRUE) matrix.expression <- read.delim(myexpressfile, header=TRUE, sep="\t", as.is=TRUE, blank.lines.skip = TRUE, fill=TRUE) cormatrix.data.raw <- read.delim(mycorrelation, sep="\t", header=TRUE, as.is=TRUE, blank.lines.skip = TRUE, fill=TRUE) listmatrix.expression <- list(matrix.expression) listcormatrix.data.raw <- list(cormatrix.data.raw) comet(config.file=configfile, mydata.file=matrix.dnamethylation, mydata.type="dataframe",cormatrix.file=listcormatrix.data.raw, cormatrix.type="listdataframe",cormatrix.sig.level=0.05, cormatrix.conf.level=0.05, cormatrix.adjust="BH", mydata.large.file=listmatrix.expression, mydata.large.type="listdataframe", fontsize.gviz =12, tracks.gviz=listgviz,verbose=FALSE, print.image=FALSE) @ \caption{Plot with \Rfunction{comet} function from matrix data and with a pvalue threshold for the correlation between omics features (here CpG sites).\label{fig:cometPlotMatrix}} \end{figure} \clearpage \subsection{\Biocpkg{coMET} plot: annotation tracks and correlation matrix} It is possible to visualise only annotation tracks and the correlation between genetic elements. In this case, we need to use the option \Robject{disp.pvalueplot=FALSE}, for example see Figure \ref{fig:cometPlotNopval}. <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) configfile <- file.path(extdata, "config_cyp1b1_zoom_4cometnopval.txt") myinfofile <- file.path(extdata, "cyp1b1_infofile.txt") mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") chrom <- "chr2" start <- 38290160 end <- 38303219 gen <- "hg19" strand <- "*" genetrack <-genes_ENSEMBL(gen,chrom,start,end,showId=FALSE) snptrack <- snpBiomart_ENSEMBL(gen, chrom, start, end, dataset="hsapiens_snp_som",showId=FALSE) strutrack <- structureBiomart_ENSEMBL(chrom, start, end, strand, dataset="hsapiens_structvar_som") clinVariant<-ClinVarMain_UCSC(gen,chrom,start,end) clinCNV<-ClinVarCnv_UCSC(gen,chrom,start,end) gwastrack <-GWAScatalog_UCSC(gen,chrom,start,end) geneRtrack <-GeneReviews_UCSC(gen,chrom,start,end) listgviz <- list(genetrack,snptrack,strutrack,clinVariant, clinCNV,gwastrack,geneRtrack) comet(config.file=configfile, mydata.file=myinfofile, mydata.type="file", cormatrix.file=mycorrelation, cormatrix.type="listfile", fontsize.gviz =12, tracks.gviz=listgviz, verbose=FALSE, print.image=FALSE, disp.pvalueplot=FALSE) @ \begin{figure} <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) configfile <- file.path(extdata, "config_cyp1b1_zoom_4cometnopval.txt") #configfile <- "../inst/extdata/config_cyp1b1_zoom_4comet.txt" myinfofile <- file.path(extdata, "cyp1b1_infofile.txt") mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") chrom <- "chr2" start <- 38290160 end <- 38303219 gen <- "hg19" strand <- "*" data(geneENSEMBLtrack) data(snpBiomarttrack) data(strucBiomarttrack) data(ClinVarCnvTrack) data(clinVarMaintrack) data(GWASTrack) data(GeneReviewTrack) listgviz <- list(genetrack,snptrack,strutrack,clinVariant, clinCNV,gwastrack,geneRtrack) comet(config.file=configfile, mydata.file=myinfofile, mydata.type="file", cormatrix.file=mycorrelation, cormatrix.type="listfile", fontsize.gviz =12, tracks.gviz=listgviz, verbose=FALSE, print.image=FALSE, disp.pvalueplot=FALSE) @ \caption{Plot with \Rfunction{comet} function without pvalue plot.\label{fig:cometPlotNopval}} \end{figure} \clearpage \subsection{\Biocpkg{coMET} plot: Manhattan plot and anonation track} It is possible to visualise only The Manhattan plot and the annotation tracks. In this case, we need to use the option \Robject{disp.cormatrixmap = FALSE}, for example see Figure \ref{fig:cometPlotNomatrix}. <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) configfile <- file.path(extdata, "config_cyp1b1_zoom_4nomatrix.txt") myinfofile <- file.path(extdata, "cyp1b1_infofile.txt") mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") chrom <- "chr2" start <- 38290160 end <- 38303219 gen <- "hg19" strand <- "*" genetrack <-genes_ENSEMBL(gen,chrom,start,end,showId=FALSE) snptrack <- snpBiomart_ENSEMBL(gen, chrom, start, end, dataset="hsapiens_snp_som",showId=FALSE) strutrack <- structureBiomart_ENSEMBL(chrom, start, end, strand, dataset="hsapiens_structvar_som") clinVariant<-ClinVarMain_UCSC(gen,chrom,start,end) clinCNV<-ClinVarCnv_UCSC(gen,chrom,start,end) gwastrack <-GWAScatalog_UCSC(gen,chrom,start,end) geneRtrack <-GeneReviews_UCSC(gen,chrom,start,end) listgviz <- list(genetrack,snptrack,strutrack,clinVariant, clinCNV,gwastrack,geneRtrack) comet(config.file=configfile, mydata.file=myinfofile, mydata.type="file", cormatrix.file=mycorrelation, cormatrix.type="listfile", fontsize.gviz =12, font.factor=3, tracks.gviz=listgviz, verbose=FALSE, print.image=FALSE) @ \begin{figure*} <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) configfile <- file.path(extdata, "config_cyp1b1_zoom_4nomatrix.txt") myinfofile <- file.path(extdata, "cyp1b1_infofile.txt") mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") chrom <- "chr2" start <- 38290160 end <- 38303219 gen <- "hg19" strand <- "*" data(geneENSEMBLtrack) data(snpBiomarttrack) data(strucBiomarttrack) data(ClinVarCnvTrack) data(clinVarMaintrack) data(GWASTrack) data(GeneReviewTrack) listgviz <- list(genetrack,snptrack,strutrack,clinVariant, clinCNV,gwastrack,geneRtrack) comet(config.file=configfile, mydata.file=myinfofile, mydata.type="file", cormatrix.file=mycorrelation, cormatrix.type="listfile", fontsize.gviz =12, font.factor=3, tracks.gviz=listgviz, verbose=FALSE, print.image=FALSE) @ \caption{Plot with \Rfunction{comet} function without the correlation matrix.\label{fig:cometPlotNomatrix}} \end{figure*} \clearpage \section{Extract the significant correlations between omic features} \Biocpkg{coMET} can help to visualise the correlations between omic features with EWAS results and other omic data. In addition, a function \Rfunction{comet.list} can extract the significant correlations according the method (options: \Robject{cormatrix.method}) and significance level (option: \Robject{cormatrix.sig.level}). The output file has 7 columns: \begin{enumerate} \item the name of the first omic feature \item the name of the second omic feature \item the correlation between the omic features \item the alpha/2 lower value (e.g. 0.05 (option \Robject{cormatrix.conf.level})) \item the alpha/2 upper value (e.g. 0.05 (option \Robject{cormatrix.conf.level})) \item the pvalue \item the pvalue adjusted with the method selected (e.g. Benjamin and Hochberg) (option \Robject{cormatrix.adjust}) \end{enumerate} <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) mycorrelation <- file.path(extdata, "cyp1b1_res37_rawMatrix.txt") myoutput <- file.path(extdata, "cyp1b1_res37_cormatrix_list_BH05.txt") comet.list(cormatrix.file=mycorrelation,cormatrix.method = "spearman", cormatrix.format= "raw", cormatrix.conf.level=0.05, cormatrix.sig.level= 0.05, cormatrix.adjust="BH", cormatrix.type = "listfile", cormatrix.output=myoutput, verbose=FALSE) listcorr <- read.csv(myoutput, header = TRUE, sep = "\t", quote = "") dim(listcorr) head(listcorr) @ \clearpage \section{Annotation tracks} Annotation tracks can be created with \Biocpkg{Gviz} using four different functions: \begin{enumerate} \item \Rfunction{UcscTrack}. Different UCSC tracks can be selected for visualisation from the table Browser of UCSC \url{http://genome-euro.ucsc.edu/cgi-bin/hgTables?hgsid=202842745\_Dlvit14QO0G6ZPpLoEVABG8aqfrm&clade=mammal&org=Human&db=hg19&hgta_group=varRep&hgta_track=cpgIslandExt&hgta_table=0&hgta\_regionType=genome&position=chr6%3A32726553\-32727053&hgta\_outputType=primaryTable&hgta\_outFileName=} \item \Rfunction{BiomartGeneRegionTrack}. A connection should be established to the Biomart database to visualise the genetic elements. \item \Rfunction{DataTrack}. This allows the visualisation of numerical data. \item \Rfunction{AnnotationTrack}. This allows the visualisation of any annotation data. \end{enumerate} For more information consult the user guide for \Biocpkg{Gviz}. \subsection{Ensembl} The Ensembl project \cite{ENSEMBL2015} produces genome databases for vertebrates and other eukaryotic species, and makes this information freely available online \url{http://www.ensembl.org/index.html}. A set of wrap R functions were created to extract data from Ensembl BioMart for human genome using Ensembl REST \cite{ENSEMBL2014}, but they can be extended to other genomes. You can ask help to \email{tiphaine.martin@mssm.edu}. This is the list of R functions created in \Biocpkg{coMET} to visualise ENSEMBL data. Below described the colors of tracks and specific characteristics of some annotation tracks. \begin{itemize} \item \Rfunction{bindingMotifsBiomart\_ENSEBML} : Visualise the binding motifs in the genomic region of interest \item \Rfunction{genes\_ENSEBML} : Visualise the genes from ENCODE in the genomic region of interest \item \Rfunction{genesName\_ENSEBML} : Visualise the name of genes from ENCODE in the genomic region of interest \item \Rfunction{interestGenes\_ENSEBML} : Visualise the genes from ENCODE in the genomic region of interest with a specific color for genes of interest \item \Rfunction{interestTranscript\_ENSEBML} : Visualise the transcripts from ENCODE in the genomic region of interest with a specific color for exons of interest \item \Rfunction{miRNATargetRegionsBiomart\_ENSEBML} : Visualise the miRNA target regions in the genomic region of interest \item \Rfunction{otherRegulatoryRegionsBiomart\_ENSEBML} : Visualise the other regulatory regions in the genomic region of interest \item \Rfunction{regulationBiomart\_ENSEBML} (obselet function): Visualise the other regulatory regions in the genomic region of interest \item \Rfunction{regulatoryEvidenceBiomart\_ENSEBML} : Visualise the regulatory evidence regions in the genomic region of interest \item \Rfunction{regulatoryFeaturesBiomart\_ENSEBML} : Visualise the regulatory features regions in the genomic region of interest \item \Rfunction{regulatorySegmentsBiomart\_ENSEBML} : Visualise the regulatory segment regions in the genomic region of interest. \warning{no more available} \item \Rfunction{snpBiomart\_ENSEBML} : Visualise the SNPs in the genomic region of interest \item \Rfunction{structureBiomart\_ENSEBML} : Visualise the structural variations in the genomic region of interest \item \Rfunction{transcript\_ENSEBML} : Visualise the transcripts in the genomic region of interest \end{itemize} Below described the colors of tracks and specific characteristics of some annotation tracks. \subsubsection{Genes and transcripts from Ensembl} The color of the genetic elements is defined by the R package \Biocpkg{Gviz}. It is possible to chagne the colour of some exsons by using the function \Rfunction{interestGenesENSEMBL} or \Rfunction{interestTranscriptENSEMBL}. The elements and the colours to be displayed must be given as list. An example is given below: <>= gen <- "hg38" chr <- "chr15" start <- 75011669 end <- 75019876 interestfeatures <- rbind(c("75011883","75013394","bad"), c("75013932","75014410","good")) interestcolor <- list("bad"="red", "good"="green") interestgenesENSMBLtrack<-interestGenes_ENSEMBL(gen,chr,start,end,interestfeatures, interestcolor,showId=TRUE) plotTracks(interestgenesENSMBLtrack, from=start, to=end) @ \begin{figure*} <>= gen <- "hg38" chr <- "chr15" start <- 75011669 end <- 75019876 data(interestgenesENSMBLtrack) plotTracks(interestgenesENSMBLtrack, from=start, to=end) @ \caption{Plot genes with different colors according user's choice.\label{fig:interestgenesENSEMBLtrack}} \end{figure*} \subsubsection{Regulatory elements from Ensembl} This function is now obselet in \Biocpkg{coMET} as Ensembl have restructured their databases due to the new version of the genome GRCh38. The same data is now available by using the function \Rfunction{RegulatoryFeaturesBiomart}. The colors were : \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/RegulatoryElementsENSEMBL}} \subsubsection{structureBiomart from Ensembl} Listed below are the colours for somatic structural variation and structural variation. \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/structureBiomart}} \subsubsection{miRNA Target Regions from Ensembl} The colour of the miRNA target regions is set to Plum4 (hex code: \#8B668B) \subsubsection{Binding Motif Biomart from Ensembl} Listed on the next page are the colours used for the different types of binding motifs. The frequency shown is that found in GRCh38 (hg38). Motifs with red text are found only in GRCh37 (hg19), motifs with blue text are found only in GRCh38 (hg38) \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/BindingMotifsBiomart}} \subsubsection{Other Regulatory Regions Biomart from Ensembl} Listed below are the colours used for the different types of regulatory regions. The frequency shown is that found in GRCh38 (hg38). \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/OtherRegulatoryRegions}} \subsubsection{Regulatory Features Biomart from Ensembl} Listed below are the colours used for the different types of regulatory features The frequency shown is that found in GRCh38 (hg38). \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/RegulatoryFeaturesBiomart}} \subsubsection{Other Regulatory Segments Biomart from Ensembl} \warning(No more available) Listed below are the colours used for the different types of regulatory segments. The frequency shown is that found in GRCh38 (hg38). Segments with red text are found only in GRCh37 (hg19) \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/RegulatorySegmentsBiomart}} \subsubsection{Regulatory Evidence Elements Biomart from Ensembl} Listed on the next 3 pages are the colours used for the different types of regulatory evidence elements. The frequency shown is that found in GRCh37 (hg19). At the current time this track has not been optimised for GRCh38 (hg38) meaning any elements found exclusively in GRCh38 do not have an assigned colour and will be displayed in the default track colour of \Biocpkg{Gviz}. \clearpage \centerline{\includegraphics{../inst/extdata/JpegTables/RegulatoryEvidence_1}} \clearpage \centerline{\includegraphics{../inst/extdata/JpegTables/RegulatoryEvidence_2}} \clearpage \centerline{\includegraphics{../inst/extdata/JpegTables/RegulatoryEvidence_3}} \clearpage \subsection{UCSC} the UCSC Genome Browser \cite{UCSC2002} website \url{http://genome-euro.ucsc.edu/} contains the reference sequence and working draft assemblies for a large collection of genomes. This is the list of \R{} wrapping functions of some tracks found in UCSC genome browser. Below described the colors of tracks and specific characteristics of some annotation tracks. \begin{itemize} \item \Rfunction{chromatinHMMAll\_UCSC} : Visualise the chromHMM Broad found in UCSC genome browser of all tissues in the genomic region of interest. \item \Rfunction{chromatinHMMOne\_UCSC} : Visualise the chromHMM Broad found in UCSC genome browser of the tissue of interest in the genomic region of interest. \item \Rfunction{ClinVarCnv\_UCSC} : Visualise clinical CNVs found in ClinVar tracks of UCSC genome browser in the genomic region of interest. \item \Rfunction{ClinVarMain\_UCSC} : Visualise clinical SNPs found in ClinVar tracks of UCSC genome browser in the genomic region of interest. \item \Rfunction{CoreillCNV\_UCSC} : Visualise CNV found in Coreil tracks of UCSC genome browser in the genomic region of interest. \item \Rfunction{COSMIC\_UCSC} : Visualise SNPs found in COSMIC tracks of UCSC genome browser in the genomic region of interest.\warning{We could not more access to COSMIC data from UCSC genome browser, people needs to extract data from COSMIC directly}. \item \Rfunction{cpgIslands\_UCSC} : Visualise CpG Island found in CpGIsland tracks of UCSC genome browser in the genomic region of interest. \item \Rfunction{DNAse\_UCSC} : Visualise clinical CNV found in ClinVar tracks of UCSC genome browser in the genomic region of interest. \item \Rfunction{GAD\_UCSC} : Visualise genes found in GAD tracks of UCSC genome browser in the genomic region of interest. \item \Rfunction{gcContent\_UCSC} : Visualise GC content found in UCSC genome browser in the genomic region of interest. \item \Rfunction{GeneReviews\_UCSC} : Visualise clinical genes found in GeneReviews tracks of UCSC genome browser in the genomic region of interest. \item \Rfunction{GWAScatalog\_UCSC} : Visualise SNPS found in GWAS catalog tracks of UCSC genome browser in the genomic region of interest. \item \Rfunction{HistoneAll\_UCSC} : Visualise histone patterns found in UCSC genome browser of all tissues in the genomic region of interest. \item \Rfunction{HistoneOne\_UCSC} : Visualise histone patterns found in UCSC genome browser of one tissue of interest in the genomic region of interest. \item \Rfunction{ISCA\_UCSC} (obselete) : Visualise clinical CNV found in UCSC genome browser in the genomic region of interest. \item \Rfunction{knownGenes\_UCSC} : Visualise known genes found in UCSC genome browser in the genomic region of interest. \item \Rfunction{refGenes\_UCSC} : Visualise reference genes found in UCSC genome browser in the genomic region of interest. \item \Rfunction{repeatMasker\_UCSC} : Visualise repeat elements found in UCSC genome browser in the genomic region of interest. \item \Rfunction{segmentalDups\_UCSC} : Visualise segmental duplcations found in UCSC genome browser in the genomic region of interest. \item \Rfunction{snpLocations\_UCSC} : Visualise SNPs found in UCSC genome browser in the genomic region of interest. \item \Rfunction{xenorefGenes\_UCSC} : Visualise xeno reference genes found in UCSC genome browser in the genomic region of interest. \end{itemize} \subsubsection{ChromHMM from UCSC} For this function there are two possible colour schemes to choose from. The selection between schemes is made with the variable \Robject{colour}. The default scheme is \Biocpkg{coMET}, the colours chosen have been selected so that different elements can be easily distinguished. The second scheme is \emph{UCSC}, these are the set colours used by UCSC, in certain plots it may be difficult to distinguish elements apart. These UCSC colours are correct at the time this document was writtern however if these change in the future and this is not reflected here please contact us. the colours used in both schemes are listed below: \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/ChromHMM_coMET}} \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/ChromHMM_UCSC}} % % \begin{longtable}{|c|c|} % \hline \multicolumn{1}{|c|}{Omic feature} & \multicolumn{1}{c|}{Color} \\ \hline % \endfirsthead % % \multicolumn{2}{c}% % {\tablename\ \thetable\ -- continued from previous page} \\ % \hline \multicolumn{1}{|c|}{Omic feature} & % \multicolumn{1}{c|}{Color} \\ \hline % \endhead % % \hline \multicolumn{2}{|r|}{{Continued on next page}} \\ \hline % \endfoot % % \hline \hline % \endlastfoot % 1\_Active\_Promoter & firebrick1 \\ % 2\_Weak\_Promoter & darksalmon \\ % 3\_Poised\_Promoter & blueviolet \\ % 4\_Strong\_Enhancer & Orange \\ % 5\_Strong\_Enhancer & coral \\ % 6\_Weak\_Enhancer & yellow \\ % 7\_Weak\_Enhancer & gold \\ % 8\_Insulator & cornflowerblue \\ % 9\_Txn\_Transition & darkolivegreen \\ % 10\_Txn\_Elongation & forestgreen \\ % 11\_Weak\_Txn & darkseagreen1 \\ % 12\_Repressed & gainsboro \\ % 13\_Heterochrom/lo & gray74 \\ % 14\_Repetitive/CNV & gray77 \\ % 15\_Repetitive/CNV & gray86 \\ % \end{longtable} \subsubsection{ISCA track (obselete database)} International Standards of Cytogenomic Arrays Consortium defined a set of phenotypes for CNVs. Different colours are defined to represent them. This database is not more accessible from UCSC (from September 2015). \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/ISCATrack}} \subsubsection{Other potential data from UCSC} You can access to other data via UCSC track hub \cite{UCSC2013} : \begin{itemize} \item Other tracks and table accessible to UCSC genome browser \url{https://genome.ucsc.edu/cgi-bin/hgTables?hgsid=444062899_lxuSrw4J9exVt1OafMuY4LDbVs1F&clade=mammal&org=Human&db=hg19&hgta_group=allTracks&hgta_track=knownGene&hgta_table=0&hgta_regionType=genome&position=chr21%3A33031597-33041570&hgta_outputType=primaryTable&hgta_outFileName=} \item Track HUB of UCSC genome browser \url{https://genome-euro.ucsc.edu/cgi-bin/hgHubConnect?hubUrl=http%3A%2F%2Ffantom.gsc.riken.jp%2F5%2Fdatahub%2Fhub.txt&hgHubConnect.remakeTrackHub=on&redirect=manual&source=genome.ucsc.edu} \end{itemize} and use \Rfunction{DataTrack} or \Rfunction{AnnotationTrack} or \Rfunction{UCSCTrack} of \Biocpkg{Gviz} to visuaslise them. \clearpage \subsection{NIH Roadmap epigenomics project} NIH Roadmap epigenomics projects \url{http://www.roadmapepigenomics.org/} \cite{ROADMAPConsortium2015} aims to produce a public resource of human epigenomic data to catalyze basic biology and disease-oriented research. The project has generated high-quality, genome-wide maps of several key histone modifications, chromatin accessibility, DNA methylation and mRNA expression across 100s of human cell types and tissues (111 consolidated epigenomes from the NIH Roadmap Epigenomics Project and 16 epigenomes from The Encyclopedia of DNA Elements (ENCODE) project). Release 9 of the compendium contains uniformly pre-processed and mapped data from multiple profiling experiments (technical and biological replicates from multiple individuals and/or datasets from multiple centers) spanning 183 biological samples and 127 consolidated epigenomes. More information on each type data are on the site of NIH Roadmap Epigenomics Program \url{http://egg2.wustl.edu/roadmap/web_portal/index.html} and the meta-data on different tissues (more for correspondance between Epigenome ID (EID) and the standartized epigenome name), you need to look at this spreadsheet \url{https://docs.google.com/spreadsheets/d/1yikGx4MsO9Ei36b64yOy9Vb6oPC5IBGlFbYEt-N6gOM/edit#gid=15} The current data are done on Release 9. The data are mapped on the reference genome \textbf{hg19}. Below described the colors of tracks and specific characteristics of some annotation tracks. \begin{itemize} \item \Rfunction{chromHMM\_RoadMap} : Visualisation of chromatin states defined in NIH Roadmap project \item \Rfunction{dgfootprints\_RoadMap}: Visualisation of DNA motif positional bias in digital genomic Footprinting Sites \item \Rfunction{DNaseI\_RoadMap} : Visualisation of promoter/enhancer regions \end{itemize} \subsubsection{Chromatin state} There are 3 chromatin states defined in NIH Roadmap project (15 states, 18 states and 25 states). For 18 and 25 states, there are the choice beteen 2 set of colors. First, the colors defined by NIH Roadmap and second, the colors defined by us for a better differentiation between states. you can use \Rfunction{chromHMM\_RoadMap} to visualise chromatin state in : \begin{itemize} \item 15-states, go to \url{http://egg2.wustl.edu/roadmap/data/byFileType/chromhmmSegmentations/ChmmModels/coreMarks/jointModel/final/} and select the MNEMONICS BED FILES, where bins with the same state label are merged and a label is assigned to the entire merged regions, related to your tissue of interest. \item 18-states, go to \url{http://egg2.wustl.edu/roadmap/data/byFileType/chromhmmSegmentations/ChmmModels/core_K27ac/jointModel/final/} and select the MNEMONICS BED FILES, where bins with the same state label are merged and a label is assigned to the entire merged regions, related to your tissue of interest . \item 25-states, go to \url{http://egg2.wustl.edu/roadmap/data/byFileType/chromhmmSegmentations/ChmmModels/imputed12marks/jointModel/final/} and select your tissue of interest. \end{itemize} You can have more information about these data from NIH Roadmap website \url{http://egg2.wustl.edu/roadmap/web_portal/chr_state_learning.html#core_15state}. You can visualise this bed using the function \Rfunction{chromHMM\_RoadMap} and you can choice the color between \emph{roadmap15}, \emph{roadmap18}, \emph{comet18}, \emph{roadmap25} and \emph{comet25}. Below you can find the color code for each state depending if 15-,18- or 25-state Listed below are the colours used for the different elements contained in NIH Roadmap data with 15 states. \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/RoadMap15_RoadMap}} Listed below are the colours used for the different elements contained in NIH Roadmap data with 18 states with NIH Roadmap colors. \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/RoadMap18_RoadMap}} Listed below are the colours used for the different elements contained in NIH Roadmap data with 18 states with \Biocpkg{coMET} colors. \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/RoadMap18_coMET}} Listed below are the colours used for the different elements contained in NIH Roadmap data with 25 states with NIH Roadmap colors. \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/RoadMap25_RoadMap}} Listed below are the colours used for the different elements contained in NIH Roadmap data with 25 states with \Biocpkg{coMET} colors. \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/RoadMap25_coMET}} \subsubsection{DNA Motif Positional Bias in Digital Genomic Footprinting Sites} The Digital Genomic Footprinting (DGF) sites in each cell type can be visualised using the function \Rfunction{dgfootprints\_RoadMap} using the file of DNase/DGF Footprint calls \url{http://egg2.wustl.edu/roadmap/data/byDataType/dgfootprints/} \subsubsection{DNaseI-accessible regulatory regions} Using the core 15-state chromatin state model across any of the 111 NIH Roadmap reference epigenomes, and focusing on states TssA, TssAFlnk, and TssBiv for promoters, and EnhG, Enh, and EnhBiv for enhancers, and state BivFlnk (flanking bivalent Enh/Tss) for ambiguous regions, 3 set of data were constructed. The data can be visualised using the function \Rfunction{DNaseI\_RoadMap} with the good name of data (variable \Robject{featureDisplay}) like in Fig. \ref{fig:cometPlotfile}: \begin{itemize} \item for \textbf{promoter} regions the file of tissue of interest \url{http://egg2.wustl.edu/roadmap/data/byDataType/dnase/BED\_files\_prom/} or RData files containing matrice of chromatin state call for promoter. Thus, user can select for different tissues. \item for \textbf{enhancer} regions the file of tissue of interest \url{http://egg2.wustl.edu/roadmap/data/byDataType/dnase/BED\_files\_enh/} \item for \textbf{dyadic} promoter/enhancer region the file of tissue of interest \url{http://egg2.wustl.edu/roadmap/data/byDataType/dnase/BED\_files\_dyadic/} \end{itemize} <>= chr<-"chr2" start <- 38290160 end <- 38303219 gen<-"hg19" extdata <- system.file("extdata", package="coMET",mustWork=TRUE) prombedFilePath <- file.path(extdata, "/RoadMap/regions_prom_E001.bed") promRMtrack<- DNaseI_RoadMap(gen,chr,start, end, prombedFilePath, featureDisplay='promotor', type_stacking="squish") enhbedFilePath <- file.path(extdata, "/RoadMap/regions_enh_E001.bed") enhRMtrack<- DNaseI_RoadMap(gen,chr,start, end, enhbedFilePath, featureDisplay='enhancer', type_stacking="squish") dyabedFilePath <- file.path(extdata, "/RoadMap/regions_dyadic_E001.bed") dyaRMtrack<- DNaseI_RoadMap(gen,chr,start, end, dyabedFilePath, featureDisplay='dyadic', type_stacking="squish") genetrack <-genes_ENSEMBL(gen,chr,start,end,showId=TRUE) listRoadMap <- list(genetrack,promRMtrack,enhRMtrack,dyaRMtrack) plotTracks(listRoadMap, chromosome=chr,from=start,to=end) @ \begin{figure*} <>= chr<-"chr2" start <- 38290160 end <- 38303219 gen<-"hg19" data(promRMtrack) data(enhRMtrack) data(dyaRMtrack) data(genetrack4RM) listRoadMap <- list(genetrack,promRMtrack,enhRMtrack,dyaRMtrack) plotTracks(listRoadMap, chromosome=chr,from=start,to=end) @ \caption{Plot of NIH Roadmap data.\label{fig:RoadMaptrack}} \end{figure*} \subsubsection{Processed data and Imputed data} BED and BigWIG file can be visualised with DataTrack objects from files of \Biocpkg{Gviz} package. The data are in \url{http://www.genboree.org/EdaccData/Release-9/sample-experiment/} and \url{http://www.genboree.org/EdaccData/Release-9/experiment-sample/} or go to \url{http://egg2.wustl.edu/roadmap/web_portal/processed_data.html} for processed data or to \url{http://egg2.wustl.edu/roadmap/web_portal/imputed.html#imp_sig} for imputed data. \clearpage \subsection{ENCODE and GENCODE data} The ENCODE (Encyclopedia of DNA Elements) Consortium is an international collaboration of research groups funded by the National Human Genome Research Institute (NHGRI) \url{https://www.encodeproject.org/}. The goal of ENCODE is to build a comprehensive parts list of functional elements in the human genome, including elements that act at the protein and RNA levels, and regulatory elements that control cells and circumstances in which a gene is active. Genes and transcripts of GENCODE are accessible from ENSEMBL biomart or can be visualised wtith \Rfunction{GeneRegionTrack} of \Biocpkg{Gviz}. Other data are in BED or BAM format that can be visualised with \Biocpkg{Gviz} tracks. <>= #Genes from GENCODE chr<-3 start <- 132239976 end <- 132541303 gen<-"hg19" extdata <- system.file("extdata", package="coMET",mustWork=TRUE) gtfFilePath <- file.path(extdata, "/GTEX/gencode.v19.genes.patched_contigs.gtf") options(ucscChromosomeNames=FALSE) grtrack <- GeneRegionTrack(range=gtfFilePath ,chromosome = chr, start= start, end= end, name = "Gencode V19", collapseTranscripts=TRUE, showId=TRUE,shape="arrow") plotTracks(grtrack, chromosome=chr,from=start,to=end) @ \begin{figure*} <>= #Genes from GENCODE chr<-3 start <- 132239976 end <- 132541303 gen<-"hg19" data(genesGencodetrack) options(ucscChromosomeNames=FALSE) plotTracks(grtrack, chromosome=chr,from=start,to=end) options(ucscChromosomeNames=TRUE) @ \caption{Plot of genes defined by GeneCode.\label{fig:GeneCodetrack}} \end{figure*} \subsubsection{Predicting motifs and active regulators} You can browse known and discovered motifs for the ENCODE TF ChIP-seq datasets. The position of motifs can be visualised using the function \Rfunction{ChIPTF\_ENCODE} using one of files from \url{http://compbio.mit.edu/encode-motifs/} \cite{Kheradpour2014} such as \url{http://compbio.mit.edu/encode-motifs/matches.txt.gz} <>= #TF Chip-seq data gen <- "hg19" chr<-"chr1" start <- 1000 end <- 329000 extdata <- system.file("extdata", package="coMET",mustWork=TRUE) bedFilePath <- file.path(extdata, "ENCODE/motifs1000_matches_ENCODE.txt") motif_color <- file.path(extdata, "ENCODE/TFmotifs_colors.csv") chipTFtrack <- ChIPTF_ENCODE(gen,chr,start, end, bedFilePath, featureDisplay=c("AHR::ARNT::HIF1A_1", "AIRE_1","AIRE_2","AHR::ARNT_1"), motif_color,type_stacking="squish",showId=TRUE) plotTracks(chipTFtrack, chromosome=chr,from=start,to=end) @ \begin{figure*} <>= #TF Chip-seq data gen <- "hg19" chr<-"chr1" start <- 1000 end <- 329000 data(chipTFtrack) plotTracks(chipTFtrack, from = start, to = end) @ \caption{Plot ENCODE TF ChIP-seq datasets of ENCODE.\label{fig:ENCODEtrack}} \end{figure*} \clearpage \subsection{GTEx Portal} The Genotype-Tissue Expression (GTEx) \cite{GTEX2013} project aims to provide to the scientific community a resource with which to study human gene expression and regulation and its relationship to genetic variation. By analyzing global RNA expression within individual tissues and treating the expression levels of genes as quantitative traits, variations in gene expression that are highly correlated with genetic variation can be identified as expression quantitative trait loci, or eQTLs. The data are accessible via \url{http://www.gtexportal.org/}. A set of data are downloadable from \url{http://www.gtexportal.org/home/datasets2} (need to have login). The data were mapped on the reference genome \textbf{hg19}. Below described the colors of tracks and specific characteristics of some annotation tracks. 2 functions were created to visualise data from GTEx version 6: \begin{enumerate} \item \Rfunction{eQTL\_GTEx} visualise eGene and significant snp-gene associations based on permutations in a tissue specific. The name of folder in GTEx version 6 is \emph{GTEx\_Analysis\_V6\_eQTLs.tar.gz}. \item \Rfunction{geneExpression\_GTEx} (need to update) visualise fully processed, normalized and filtered gene expression data, which was used as input into Matrix eQTL for eQTL discovery in a tissue specific. The name of folder in GTEx version 6 is \emph{GTEx\_Analysis\_V6\_eQTLInputFiles\_geneLevelNormalizedExpression.tar.gz} \end{enumerate} One function from \Biocpkg{Gviz}: \begin{enumerate} \item \Rfunction{GeneRegionTrack} can visualise gene level model based on the GENCODE transcript model (cf. example below. Isoforms have been collapsed to single genes. The name of file in GTEx version 6 is \emph{gencode.v19.genes.patched\_contigs.gtf}. \end{enumerate} <>= ## eQTL data chr<-"chr3" start <- 132239976 end <- 132541303 gen<-"hg19" extdata <- system.file("extdata", package="coMET",mustWork=TRUE) bedFilePath <- file.path(extdata, "/GTEX/eQTL_Uterus_Analysis_extract100.snpgenes") eGTex<- eQTL_GTEx(gen,chr, start, end, bedFilePath, featureDisplay = 'all', showId=TRUE, type_stacking="squish", just_group="left" ) eGTex_SNP<- eQTL_GTEx(gen,chr, start, end, bedFilePath, featureDisplay = 'SNP', showId=FALSE, type_stacking="dense", just_group="left") #Genes from gtfFilePath <- file.path(extdata, "/GTEX/gencode.v19.genes.patched_contigs.gtf") options(ucscChromosomeNames=FALSE) grtrack <- GeneRegionTrack(genome="hg19",range=gtfFilePath ,chromosome = chr, start= start, end= end, name = "Gencode V19", collapseTranscripts=TRUE, showId=TRUE,shape="arrow") eGTexTracklist <- list(grtrack,eGTexTrackSNP) plotTracks(eGTexTracklist, chromosome=chr,from=start,to=end) @ \begin{figure*} <>= ## eQTL data chr<-"chr3" start <- 132239976 end <- 132541303 gen<-"hg19" data(eGTexTrackSNP) data(eGTexTrackall) data(grtrack4eGTex) #Genes from eGTexTracklist <- list(grtrack,eGTexTrackSNP) plotTracks(eGTexTracklist, chromosome=chr,from=start,to=end) @ \caption{Plot eQTL from GTex.\label{fig:eQTLGTextrack}} \end{figure*} 2 other functions were created to visualise supplement data from GTEx version 3 \begin{enumerate} \item \Rfunction{psiQTL\_GTEx} visualise results from the protein truncating variants QTL (psiQTL) analysis for mine main tissues, plus brain, plus multi-tissue that averages the exons where data for three or more tissues is available. The name of file in GTEX version 3 is \emph{gtex\_psiqtls.zip}. \item \Rfunction{imprintedGenes\_GTEx} visuaslise gene imprinting genes in different tissues \cite{Baran2015} via url \url{http://www.gtexportal.org/home/imprintingPage}. There are 33 tissues and 5 classification \end{enumerate} <>= ### psiQTL chr<-"chr13" start <- 52713837 end <- 52715894 gen<-"hg19" extdata <- system.file("extdata", package="coMET",mustWork=TRUE) psiQTLFilePath <- file.path(extdata, "/GTEX/psiQTL_Assoc-total.AdiposeTissue.txt") psiGTex<- psiQTL_GTEx(gen,chr,start, end, psiQTLFilePath, featureDisplay = 'all', showId=TRUE, type_stacking="squish",just_group="above" ) genetrack <-genes_ENSEMBL(gen,chr,start,end,showId=TRUE) psiTrack <- list(genetrack,psiGTex) plotTracks(psiTrack, chromosome=chr,from=start,to=end) @ \begin{figure*} <>= ### psiQTL chr<-"chr13" start <- 52713837 end <- 52715894 gen<-"hg19" data(psiGTexTrackall) data(genetrack4psiGTEX) psiTrack <- list(genetrack,psiGTexTrackall) plotTracks(psiTrack, chromosome=chr,from=start,to=end) @ \caption{Plot psiQTL from GTex.\label{fig:psiQTLGTextrack}} \end{figure*} <>= data(imprintedGenesGTEx) as.character(unique(imprintedGenesGTEx$Tissue.Name)) as.character(unique(imprintedGenesGTEx$Classification)) @ <>= ### inprinted genes chr<- "chr1" start <- 7895752 end <- 7914572 gen<-"hg19" genesTrack <- genes_ENSEMBL(gen,chr,start,end,showId=TRUE) allIG <- imprintedGenes_GTEx(gen,chr,start, end, tissues="all", classification="imprinted",showId=TRUE) allimprintedIG <- imprintedGenes_GTEx(gen, chr,start, end, tissues="all", classification="imprinted",showId=TRUE) StomachIG <-imprintedGenes_GTEx(gen,chr,start, end, tissues="Stomach", classification="all",showId=TRUE) PancreasIG <- imprintedGenes_GTEx(gen,chr,start, end, tissues="Pancreas", classification="all",showId=TRUE) PancreasimprintedIG <- imprintedGenes_GTEx(gen,chr,start, end, tissues="Pancreas", classification="imprinted",showId=TRUE) plotTracks(list(genesTrack, allIG, allimprintedIG, StomachIG,PancreasIG,PancreasimprintedIG), chromosome=chr, from=start, to=end) @ \begin{figure*} <>= ### inprinted genes chr<- "chr1" start <- 7895752 end <- 7914572 gen<-"hg19" genesTrack <- genes_ENSEMBL(gen,chr,start,end,showId=TRUE) data(allIGtrack) data(allimprintedIGtrack) data(StomachIGtrack) data(PancreasIGtrack) data(PancreasimprintedIGtrack) imprintinglist <- list(genesTrack,allIGtrack,allimprintedIGtrack,StomachIGtrack, PancreasIGtrack,PancreasimprintedIGtrack) plotTracks(imprintinglist, from = start, to = end) @ \caption{Plot imprinted genes from GTex.\label{fig:IGTextrack}} \end{figure*} % <>= % ## sQTL data % chr<-3 % start <- 132239976 % end <- 132541303 % gen<-"hg19" % % extdata <- system.file("extdata", package="coMET",mustWork=TRUE) % % ### sQTL from Altran methods % sQTLAFilePath <- file.path(extdata, "/GTEX/sQTL_HeartLeftVentricle.Altrans.FDR05.bestPerLink") % % sAGTex<- sQTL_Altrans_GTEx(gen,chr,start, end, % sQTLAFilePath, featureDisplay = 'all', showId=TRUE, % type_stacking="squish",just_group="left" ) % @ \clearpage \subsection{Hi-C data} Below are examples of Hi-C data available for different tissues. \subsubsection{Hi-C data at 1kb resolution at Lieberman Aiden lab} They \cite{Rao2014} used in situ Hi-C to probe the three-dimensional architecture of genomes, constructing haploid and diploid maps of nine cell types. The densest, in human lymphoblastoid cells, contains 4.9 billion contacts, achieving 1-kilobase resolution.The data were mapped on \textbf{hg19} reference genome. You can download intrachromosomal matrice from \url{http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE63525} for the cell-type of interest. <>= library('corrplot') #Hi-C data gen <- "hg19" chr<-"chr1" start <- 5000000 end <- 9000000 extdata <- system.file("extdata", package="coMET",mustWork=TRUE) bedFilePath <- file.path(extdata, "HiC/chr1_1mb.RAWobserved") matrix_HiC <- HiCdata2matrix(chr,start, end, bedFilePath) cor_matrix_HiC <- cor(matrix_HiC) diag(cor_matrix_HiC)<-1 corrplot(cor_matrix_HiC, method = "circle") @ % \begin{figure*} % <>= % library('corrplot') % #Hi-C data % gen <- "hg19" % chr<-"chr1" % start <- 5000000 % end <- 9000000 % % data(matrix_HiC_Rao) % cor_matrix_HiC <- cor(matrix_HiC_Rao) % diag(cor_matrix_HiC)<-1 % corrplot(cor_matrix_HiC, method = "circle") % @ % \caption{plot HiC data.\label{fig:HiCrack}} % \end{figure*} You can quick visualise this data using this HiC-interaction tool \url{http://promoter.bx.psu.edu/hi-c/view.php?species=human&assembly=hg19&source=inside&tissue=GM12878&resolution=1&c_url=&gene=CTXN1&sessionID=} \subsubsection{Hi-C Data Browser} You can download heatmap of your region of interest from two cell-line GM06690 (immortalized lymphoblast) or K562 (leukemia) using their website \url{http://hic.umassmed.edu/heatmap/heatmap.php}. This data was produced by \cite{LiebermanAiden2009}. The region that you want to visualise with this data need to large more than either 100Kb or 1Mb as Heatmaps were generated by dividing the chromosome up into 100 Kb or 1 Mb windows. The data were mapped on \textbf{hg19} reference genome. You need to create info file to define the position of each bin composing your interaction matrice in using the row name of matrice as the name of bin contain the start and end of bin. \subsubsection{Hi-C project at Ren Lab} Interaction matrices for each of the four cell types analysis (mouse ES cell, mouse cortex, human ES cell (H1), and IMR90 fibroblasts) by Ren Lab (to cite them, you need to select the publication for this url \url{http://promoter.bx.psu.edu/hi-c/publications.html}) are accessible via url \url{http://chromosome.sdsc.edu/mouse/hi-c/download.html}. The interaction matrices are created using either a 40kb bin size throughout the genome. So the region that you want to visualise with this data need to large more than 40Kb. The data were mapped on \textbf{hg19} reference genome. You need to : \begin{enumerate} \item Extract from the BED file that contains the locations of each of the topological domains the region of interest \item Extract in either raw or normalised matrice only the sub-matrice of interest \end{enumerate} <>= extdata <- system.file("extdata", package="coMET",mustWork=TRUE) info_HiC <- file.path(extdata, "Human_IMR90_Fibroblast_topological_domains.txt") data_info_HiC <-read.csv(info_HiC, header = FALSE, sep = "\t", quote = "") intrachr_HiC <- file.path(extdata, "Human_IMR90_Fibroblast_Normalized_Matrices.txt") data_intrachr_HiC <- read.csv(intrachr_HiC, header = TRUE, sep = "\t", quote = "") chr_interest <- "chr2" start_interest <- "1" end_interest <- "160000" list_bins <- which(data_info_HiC[,1] == chr_interest & data_info_HiC[,2] >= start_interest & data_info_HiC[,2] <= end_interest ) subdata_info_Hic <- data_info_HiC[list_bins,] subdata_intrachr_HiC <- data_intrachr_HiC[list_bins,list_bins] @ \clearpage \subsection{FANTOM5 database} FANTOM \url{http://fantom.gsc.riken.jp/} established the FANTOM database (transcripts, transcription factors, promoters and enhancers active,TSS) and the FANTOM full-length cDNA clone bank, which are available worldwide for about 400 distinct cell types. Currently, FANTOM is in version FANTOM5 phase 2 where data were mapped on reference genome \textbf{hg19} for human or \textbf{mm9} for mouse \cite{Lizio2015}. To extract data \begin{itemize} \item from \url{http://fantom.gsc.riken.jp/5/} \item from \url{http://fantom.gsc.riken.jp/data/} or \url{http://fantom.gsc.riken.jp/views/} \item from BED file used by UCSC HUB \url{http://fantom.gsc.riken.jp/5/datahub/}, more information here \url{http://fantom.gsc.riken.jp/5/datahub/description.html} \end{itemize} As the data are in classical format such as BED file, you can use easily \Biocpkg{Gviz}'s \Rfunction{DataTrack} function to visuaslise them. However, there are some comment lines that you need to remove in the top of files. 2 functions were created : \begin{itemize} \item \Rfunction{DNaseI\_FANTOM} helps to visualise enhancer regions defined by FANTOM5 \item \Rfunction{TFBS\_FANTOM} helps to visualise TFBS regions defined by FANTOM5 \end{itemize} <>= gen <- "hg19" chr<- "chr1" start <- 6000000 end <- 6500000 extdata <- system.file("extdata", package="coMET",mustWork=TRUE) ##Enhancer enhFantomFile <- file.path(extdata, "/FANTOM/human_permissive_enhancers_phase_1_and_2.bed") enhFANTOMtrack <-DNaseI_FANTOM(gen,chr,start, end, enhFantomFile, featureDisplay='enhancer') ### TFBS motif AP1FantomFile <- file.path(extdata, "/FANTOM/Fantom_hg19.AP1_MA0099.2.sites.txt") tfbsFANTOMtrack <- TFBS_FANTOM(gen,chr,start, end, AP1FantomFile) @ \begin{figure*} <>= gen <- "hg19" chr<- "chr1" start <- 6000000 end <- 6500000 extdata <- system.file("extdata", package="coMET",mustWork=TRUE) ##Enhancer data(enhFANTOMtrack) ### TFBS motif data(tfbsFANTOMtrack) Fantom5list <- list(enhFANTOMtrack,tfbsFANTOMtrack) plotTracks(Fantom5list, from = start, to = end) @ \caption{plot FANTOM5 data.\label{fig:FANTOM5Track}} \end{figure*} \clearpage \subsection{BLUEprint project} BLUEprint \url{http://www.blueprint-epigenome.eu/} aims to further the understanding of how genes are activated or repressed in both healthy and diseased human cells. BLUEPRINT will focus on distinct types of haematopoietic cells from healthy individuals and on their malignant leukaemic counterparts. the data were mapped on reference genome partially on \textbf{GRCh37} and all on \textbf{GRCh38}. As the data are in classical format such as BED file, BigWig of GTF, you can use easily \Rfunction{DataTrack} or \Rfunction{AnnotationTrack} of \Biocpkg{Gviz} to visuaslise them. \subsection{Our data} \subsubsection{eQTL data} You can visualise our eQTL using \Rfunction{eQTL} function. Listed below are the colours used for the different elements contained in eQTL data. \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/eQTL}} \subsubsection{metQTL data} You can visualise our eQTL using \Rfunction{metQTL} function. Listed below are the colours used for the different elements contained in metQTL data. \\ \\ \centerline{\includegraphics{../inst/extdata/JpegTables/metQTL}} \clearpage \section{\Biocpkg{coMET}: \CRANpkg{Shiny} web-service} \subsection{How to use the \Biocpkg{coMET} web-service} If you want to use \Biocpkg{coMET} via its webservice, please go to \url{http://epigen.kcl.ac.uk/comet} and select one of different instances or direcly access one of the instances, for example \url{http://comet.epigen.kcl.ac.uk:3838/coMET/}. We have created different instances of \Biocpkg{coMET} because we did not have access to the pro version of \CRANpkg{Shiny}. All instances use the same version of \Biocpkg{coMET}. If you use \Biocpkg{coMET} from a \CRANpkg{Shiny} webservice, you do not need to install the \Biocpkg{coMET} package on your computer. The web service is user friendly and requires input files and configuration of the plot. The creation of the \Biocpkg{coMET} plot can take some time because it makes a live connection to UCSC or/and ENSEMBL for the annotation tracks. First, the plot is created on the webpage, and then it can be saved as an output file. For better quality plots please use the download option and the plot will be recreated in a file in pdf or eps format. \subsection{How to install the \Biocpkg{coMET} web-service} These are different steps to install \Biocpkg{coMET} on your \CRANpkg{Shiny} web-service and you need to be root to install it. \begin{enumerate} \item You need to install \R{}, \Bioconductor{} and the \Biocpkg{coMET} package under root. \item You need first to install the \CRANpkg{Shiny} and \CRANpkg{rmarkdown} R package before \CRANpkg{Shiny} Server. \begin{verbatim} sudo su - -c "R -e \"install.packages('shiny', repos='http://cran.rstudio.com/')\"" sudo su - -c "R -e \"install.packages('rmarkdown', repos='http://cran.rstudio.com/')\"" \end{verbatim} \item You can install \CRANpkg{Shiny} Server \url{http://shiny.rstudio.com/}, go to \url{https://www.rstudio.com/products/shiny/download-server/}. \begin{verbatim} sudo apt-get install gdebi-core wget \url{https://download3.rstudio.org/ubuntu-12.04/x86_64/shiny-server-1.4.2.786-amd64.deb} sudo gdebi shiny-server-1.4.2.786-amd64.deb \end{verbatim} \item \CRANpkg{Shiny} Server should now be installed and running on port 3838. You should be able to see a default welcome screen at \url{http://your\_server\_ip:3838/}. You can make sure your \CRANpkg{Shiny} Server is working properly by going to \url{http://your\_server\_ip:3838/sample-apps/hello/}. \item You now have a functioning \CRANpkg{Shiny} Server that can host \CRANpkg{Shiny} applications or interactive documents. The configuration file for \CRANpkg{Shiny} Server is at /etc/shiny-server/shiny-server.conf. By default it is configured to serve applications in the /srv/shiny-server/ directory. This means that any \CRANpkg{Shiny} application that is placed at /srv/shiny-server/app\_name will be available to the public at \url{http://your\_server\_ip:3838/app\_name/}. \item In a \CRANpkg{Shiny}'s folder (e.g. /var/shiny-server/www), you can create a folder called "COMET". \item Following this, you can install the two \Biocpkg{coMET} scripts in www of the \Biocpkg{coMET} package, within this new folder. \item You need to change owner and permissions to access this folder. Only the user called \CRANpkg{Shiny} can access it. \begin{verbatim} mkdir -p /var/shiny-server/www/COMET chmod -R 755 /var/shiny-server/www/COMET chown -R shiny:shiny /var/shiny-server/www/COMET mkdir -p /var/shiny-server/log chmod -R 755 /var/shiny-server/log chown -R shiny:shiny /var/shiny-server/log \end{verbatim} \item You need now to update the configuration file of \CRANpkg{Shiny} (e.g. /etc/shiny-server/shiny-server.conf). \item You need to change owner and the permission to access this file \begin{verbatim} chmod 744 /etc/shiny-server/shiny-server.conf chown shiny:shiny /etc/shiny-server/shiny-server.conf \end{verbatim} \item At the end, you should restart the service \CRANpkg{Shiny} via the command line: \begin{verbatim} ###2.13.0.1 systemd (RedHat 7, Ubuntu 15.04+, SLES 12+) #File to change: /etc/systemd/system/shiny-server.service #How to define the environment variable: [Service] Environment="SHINY\_LOG\_LEVEL=TRACE" #Commands to run for the changes to take effect: sudo systemctl stop shiny-server sudo systemctl daemon-reload sudo systemctl start shiny-server ###2.13.0.2 Upstart (Ubuntu 12.04 through 14.10 and RedHat 6) #File to change: /etc/init/shiny-server.conf #How to define the environment variable: env SHINY\_LOG\_LEVEL=TRACE #Commands to run for the changes to take effect: sudo restart shiny-server \end{verbatim} \end{enumerate} \clearpage Your Shiny's configuration file: \begin{verbatim} run_as shiny; # Define a top-level server which will listen on a port server { # Instruct this server to listen on port 3838 listen 3838; # Define the location available at the base URL location / { # Run this location in 'site_dir' mode, which hosts the entire directory # tree at '/srv/shiny-server' site_dir /var/shiny-server/www; # Define where we should put the log files for this location log_dir /var/shiny-server/log; # Should we list the contents of a (non-Shiny-App) directory when the user # visits the corresponding URL? directory_index off; # app_init_timeout 3600; # app_idle_timeout 3600; } } \end{verbatim} \clearpage \section{FAQs} \begin{itemize} \item\textbf{I cannot see my plot after running \Rfunction{comet} or \Rfunction{comet.web}. What should I do?} \\ If the previous time \Rfunction{comet} or \Rfunction{comet.web} ran and error was produced it prevents the plot from being closed. to fix this use the command '\emph{dev.off()}' as many times as necessary. \leavevmode \\ \item\textbf{How do we know if my track has data? and what the data is?} \\ Type the name of your track, visualise the track with \Rfunction{plotTrack} or read different parameters with \Rfunction{str} function. \begin{lstlisting} genetrack <-genesENSEMBL(gen,chrom,start, end,showId=TRUE) plotTracks(genetrack) str(genetrack) \end{lstlisting} \leavevmode \\ \item\textbf{How do you increase the size of the font of the name of an object}? \leavevmode \\ To enlarge the name of gene, as the object is \Biocpkg{Gviz} object, you can use the option from \Biocpkg{Gviz}. \\ You can see the value of different parameters via this command line: \begin{lstlisting} genetrack <-genesENSEMBL(gen,chrom,start, end,showId=TRUE) displayPars(genetrack) \end{lstlisting} \leavevmode \\ So if you want to enlarge the name of gene, you need to do use the option \Robject{fontsize.gviz} in the \Rfunction{comet} function, an example is given below: \\ \begin{lstlisting} comet(config.file = configfile, mydata.file = myinfofile, mydata.format = "file", cormatrix.file = mycorrelation, cormatrix.type = "listfile", mydata.large.file = mylargedata, mydata.large.type = "listfile", tracks.gviz = listGviz, verbose = TRUE, print.image=TRUE,fontsize.gviz=10) \end{lstlisting} \leavevmode \\ \item\textbf{Can I make a selection of which genes or transcripts to display}? \\ To make a selection of genes to display first create the track like you would if you were displaying all genes. From this track create another with only the genes you want to display like in the example below. Please note it is not possible to select genes based on their names unless the option to display gene names instead of gene reference is used, in other cases it is possible to make a selection based on the genes reference number. \\ \begin{lstlisting} geneTrack <- refGenesUCSC(gen, chr, start, end, IdType ="name", showId = TRUE) geneTrackShow <- geneTrack[gene(geneTrack) %in% c("AHRR")] \end{lstlisting} \leavevmode \\ \item\textbf{How can I better understand where the \Rfunction{comet} function stopped}? \\ Use option \emph{VERBOSE=TRUE} in the function \Rfunction{coMET} or \Rfunction{coMET.web} \leavevmode \\ If this does not help resolve the issue, please to send your command line with \emph{VERBOSE=TRUE} and its error message to \email{tiphaine.martin@mssm.edu}. Do not forget to give alsoinformation about the session by using \Rfunction{sessionInfo()}. \\ \item\textbf{How do you visualise coMET plots working within a R \CRANpkg{Markdown} or \CRANpkg{knitr} framework}? \leavevmode \\ When coMET writes to a PDF, it is writing out to a 7X7 square area. So, it turns out that one can 'force; the R \CRANpkg{Markdown} block as well as \CRANpkg{knitr} block to also write to a 7 x 7 square PDF using option for chunck \Rcode{fig.height=7,fig.width=7} , as follows: \begin{lstlisting} '''{r plot_ex1,fig.keep='last', fig.height=7, fig.width=7,dev='pdf'} comet(config.gile=configfile, mydata.file=myinfofile, mydata_type="file", cormatrix.file-mycorrelation, cormatrix.type="listfile", mydata.large.file=myexpressfile, mydata.large.type="listfile", tracks.gviz=listgviz, verbose=FALSE, print.image=FALSE, disp.pvalueplot=FALSE) ''' \end{lstlisting} \leavevmode \\ \end{itemize} \clearpage \section{Acknowledgement} T.C.M would like to thank Bioconductor team for their help and advice in the development of a R Bioconductor package. Moreover, T.C.M would like to thank different users for their feedback that help to improve this present R package. \begin{itemize} \item Prof Daniel Weeks and Dr Annie Infancia Arockiaraj to share with us how to visualise correctly coMET plot in R Markdown code. \end{itemize} \clearpage \section{SessionInfo} The following is the session info that generated this vignette: <>= toLatex(sessionInfo()) @ <>= # #Doc for Shiny Server # # https://www.digitalocean.com/community/tutorials/how-to-set-up-shiny-server-on-ubuntu-14-04 # # #Need to have the last version of R associated with Bioconductor in parallele of your original version # #install from source http://bioconductor.org/developers/how-to/useDevel/ # mkdir /home/tiphaine/Rdevel2.2/ # cd Rdevel2.2 # mkdir tar # #Extract R source in this latter folder # ./configure -prefix=/home/tiphaine/Rdevel2.2/ --enable-R-shlib=yes # # sudo make # # sudo make check # # sudo make install # # if need to recompile, do "sudo make uninstall" "sudo make clean" and rerun from ./configure line https://stackoverflow.com/questions/28096239/how-to-configure-r-3-1-2-with-enable-r-shlib # # sudo R # # #Need to have the last version of Bioconductor and BiocCheck # #Need to update under Rdev # remove.packages("BiocManager") # if (!requireNamespace("BiocManager", quietly=TRUE)) # install.packages("BiocManager") # # #Need to update different packages associated with coMET # BiocManager::install("coMET") # # #Go to the parent folder of the package # # To create the manual documentation, need to run # # need to be the folder parent of coMET # cd ../ # R-devel CMD Rd2pdf --pdf coMET # # # Need to build the vignette # cd coMET/vignettes # R-devel CMD Sweave --engine=knitr::knitr --pdf coMET.Rnw # R-devel CMD texi2dvi --pdf coMET.tex # R-devel CMD pdflatex coMET.tex # # #Need to build the new package # to addd some parameters for R CMD build/ check/install https://support.rstudio.com/hc/en-us/articles/200486518-Customizing-Package-Build-Options # # Need to be in the folder parent of coMET # cd ../../ # R-devel CMD build coMET --resave-data --no-build-vignettes # # #Need to check if the new package follow the rules of R # R-devel CMD check coMET # # #Need to install if the new package follow the rules of R # R-devel CMD INSTALL coMET # # #Need to check the examples in R after building and checking # # Need to correct any examples that fail # example for the function comet # R>library(coMET) # R>example(comet) # # #To check if the new package follow the rules from bioconductor # R-devel CMD BiocCheckdev3.3 coMET # To install the new package # R-devel CMD INSTALL coMET # Authors@R: c(person("Tiphaine C.","Martin", role=c("aut","cre"),email="tiphaine.martin@mssm.edu), # person("Thomas","Hardiman",role=c("aut")), # person("Idil","Yet",role=c("aut")), # person("Pei-Chien","Tsai",role=c("aut")), # person("Jordana T.","Bell",role=c("aut"))) #other version of R with Rstudio #https://support.rstudio.com/hc/en-us/articles/200486138-Changing-R-versions-for-RStudio-desktop @ \clearpage %\bibliographystyle{authordate1} \bibliography{biblio} \end{document}