% -*- mode: noweb; noweb-default-code-mode: R-mode; -*- %\VignetteIndexEntry{panelcn.mops: Manual for the R package} %\VignetteDepends{panelcn.mops} %\VignettePackage{panelcn.mops} %\VignetteEngine{knitr::knitr} %\VignetteKeywords{copy number analysis, mixture distribution, latent variables, % Poisson distribution, EM algorithm, NGS, CNV, copy number variant} \documentclass[article]{bioinf} \usepackage{amsmath,amssymb} \usepackage{hyperref} \usepackage{float} \usepackage[authoryear]{natbib} \hypersetup{colorlinks=false, pdfborder=0 0 0, pdftitle={panelcn.MOPS - CNV detection tool for targeted NGS panel data}, pdfauthor={Gundula Povysil}} \title{panelcn.MOPS - CNV detection tool for targeted NGS panel data} \author{Verena Haunschmid and Gundula Povysil} \affiliation{Institute of Bioinformatics, Johannes Kepler University Linz\\Altenberger Str. 69, 4040 Linz, Austria\\ \email{povysil@bioinf.jku.at}} \newcommand{\panelcnmops}{\texttt{panelcn.mops}} \newcommand{\cnmops}{\texttt{cn.mops}} \newcommand{\method}[1]{{\fontfamily{phv}\fontshape{rm}\selectfont #1}} \newcommand{\R}{R} \newcommand{\Real}{\mathbb{R}} \renewcommand{\vec}[1]{\mathbf{#1}} \setkeys{Gin}{width=0.55\textwidth} <>= library(knitr) opts_chunk$set( ) @ \begin{document} <>= library(knitr) opts_chunk$set( concordance=TRUE ) @ <>= options(width=65) set.seed(0) library(panelcn.mops) panelcn.mopsVersion <- packageDescription("panelcn.mops")$Version @ \newcommand{\cnmopsVersion}{\Sexpr{panelcn.mopsVersion}} \manualtitlepage[Version \cnmopsVersion, \today] %\section*{Scope and Purpose of this Document} % %This document is a user manual for the \R\ package \panelcn.mops. %It is only meant as a gentle introduction into how to use the basic %functions implemented in this package. Not all features of the \R\ %package are described in full detail. Such details can be obtained %from the documentation enclosed in the \R\ package. Further note %the following: (1) this is neither an introduction to CNV detection from NGS %data; (2) this is not an introduction to \R. %If you lack the background for understanding this manual, you first %have to read introductory literature on these subjects. % \vspace{1cm} \newlength{\auxparskip} \setlength{\auxparskip}{\parskip} \setlength{\parskip}{0pt} \tableofcontents \clearpage \setlength{\parskip}{\auxparskip} \newlength{\Nboxwidth} \setlength{\Nboxwidth}{\textwidth} \addtolength{\Nboxwidth}{-2\fboxrule} \addtolength{\Nboxwidth}{-2\fboxsep} \newcommand{\notebox}[1]{% \begin{center} \fbox{\begin{minipage}{\Nboxwidth} \noindent{\sffamily\bfseries Note:} #1 \end{minipage}} \end{center}} \section{Introduction} The \panelcnmops\ package is based on the \cnmops\ package and allows to detect copy number variations (CNVs) from targeted NGS panel data. Please visit \url{http://www.bioinf.jku.at/software/panelcnmops/index.html} for additional information.\par \section{Getting started and quick start} To load the package, enter the following in your \R\ session: <>= library(panelcn.mops) data(panelcn.mops) @ The whole pipeline will only take a few steps, if BAM files are available (for read count matrices directly go to step 2): \begin{enumerate} \item Getting count windows from the BED file (also see Section \ref{s:input}). <>= bed <- "Genes_part.bed" countWindows <- getWindows(bed) @ \item Getting read counts (RCs) from BAM file (also see Section \ref{s:input}). Note that the BAM file is not included so do not try to run this code. However, the resulting test object is included as part of the data. <>= testbam <- "SAMPLE1.bam" test <- countBamListInGRanges(countWindows = countWindows, bam.files = testbam, read.width = 150) @ \item Running the algorithm (also see Section \ref{s:panelcn.mops}). <>= selectedGenes <- c("ATM") XandCB <- test elementMetadata(XandCB) <- cbind(elementMetadata(XandCB), elementMetadata(control)) resultlist <- runPanelcnMops(XandCB, testiv = 1:ncol(elementMetadata(test)), countWindows = countWindows, selectedGenes = selectedGenes) @ \item Visualization of the detected CNV regions. For more information about the result objects and visualization see Section \ref{s:results} and Section \ref{s:plot}. <>= XandCB <- test elementMetadata(XandCB) <- cbind(elementMetadata(XandCB), elementMetadata(control)) selectedGenes <- "ATM" @ <<>>= sampleNames <- colnames(elementMetadata(test)) resulttable <- createResultTable(resultlist = resultlist, XandCB = XandCB, countWindows = countWindows, selectedGenes = selectedGenes, sampleNames = sampleNames) (tail(resulttable[[1]])) @ <>= plotBoxplot(result = resultlist[[1]], sampleName = sampleNames[1], countWindows = countWindows, selectedGenes = selectedGenes, showGene = 1) @ <>= sampleNames <- colnames(elementMetadata(test)) selectedGenes <- "ATM" pdf("001.pdf", width = 10) plotBoxplot(result = resultlist[[1]], sampleName = sampleNames[1], countWindows = countWindows, selectedGenes = selectedGenes, showGene = 1) dev.off() @ \begin{figure}[H] \begin{center} \includegraphics[angle=0,width= 0.9\columnwidth]{001.pdf} \end{center} \end{figure} \end{enumerate} \section{Input} \label{s:input} % \subsection{Read count matrices as input} % Like \cnmops\ \panelcnmops\ does not require the data samples to be of any % specific kind or structure. \cnmops\ only requires a {\em read count matrix}, % i.e., given $N$ data samples and $m$ genomic segments, this is an $m\times N$ % real- or integer-valued matrix $\mathbf{X}$, % in which an entry $x_{ij}$ corresponds to the read count of sample $j$ in the % $i$-th segment. E.g. in the following read count matrix sample three has % $17$ reads in the second segment: $x_{23}=71$. % % % \newlength{\mylen} % \setlength{\mylen}{0.43cm} % % \[\mathbf{X}= \begin{array}{c} \\ \mathrm{Segment\ 1} \\ \mathrm{Segment\ 2} % \\ \mathrm{Segment\ 3} \\ \mathrm{Segment\ 4}\\ \mathrm{Segment\ 5} \\ % \mathrm{Segment\ 6} \\ \hspace{0.2cm} \end{array} % \begin{array}{c} % \begin{array}{cccc} % \mathrm{Sample\ 1} & \mathrm{Sample\ 2} & \mathrm{Sample\ 3} & % \mathrm{Sample\ 4}\end{array}\\ % \left(\begin{array}{cccc} % \hspace{\mylen}88\hspace{\mylen} & \hspace{\mylen}82\hspace{\mylen} & % \hspace{\mylen}79\hspace{\mylen} & \hspace{\mylen}101\hspace{\mylen}\\ % 83 & 78 & 71 & 99\\ % 43 & 50 & 55 & 37\\ % 47 & 58 & 48 & 42 \\ % 73 & 86 & 95 & 91\\ % 92 & 90 & 80 & 71 % \end{array}\right) \\ \hspace{0.2cm} \end{array} % \] % % % \panelcnmops\ can handle numeric and integer matrices or \verb+GRanges+ % objects, in which the read counts are stored as \verb+values+ of the object. % % % % \subsection{BAM files as input} % \label{s:bam} The most widely used file format for aligned short reads is the Sequence Alignment Map (SAM) format or in the compressed form the Binary Alignment Map (BAM). \panelcnmops\ modifies the read count function \verb+countBamInGRanges+ from the \R\ package \texttt{exomeCopy} to extract read counts for a list of BAM files. The result object of the function can directly be used as input for \panelcnmops. The first step is to extract all regions of interest (ROIs) that define the count windows from a BED file with the function \verb+getWindows+. The BED file that is provided is a subset of the TruSight Cancer Panel BED file. <<>>= bed <- system.file("extdata/Genes_part.bed", package = "panelcn.mops") countWindows <- getWindows(bed) @ The BED file should have the following structure: <>= bed <- system.file("extdata/Genes_part.bed", package = "panelcn.mops") write.table(head(read.table(bed)), row.names = FALSE, col.names = FALSE, quote = FALSE) @ While the first 3 columns list chromosome name, start and end position, the fourth column needs to start with the gene name. Additional information in the fourth column needs to be separated with a dot and may include the exon number and further information. By default the "chr" prefix of the chromosome name is removed if present. This can be changed by setting the {\tt chr} parameter to TRUE. If a mismatch of chromosome naming between the \verb+countWindows+ object and the BAM files is detected, the naming convention of the BAM file is chosen. In the second step RCs are generated from the BAM files. The read.width parameter reflects the typical length of the reads that should be counted. Note that the BAM file is not included so do not try to run this code. However, the resulting test object is included as part of the data. <>= testbam <- "SAMPLE1.bam" test <- countBamListInGRanges(countWindows = countWindows, bam.files = testbam, read.width = 150) @ In \verb+test+ you have now stored the genomic segments (left of the $\mid$'s) and the read counts (right of the $\mid$'s): <<>>= (test) @ If the BED file contains very large ROIs, a higher resolution of the CNV detection algorithm can be achieved by splitting up larger ROIs into smaller overlapping bins. This can be achieved with the funciton \verb+splitROIs+: <>= splitROIs(bed, "newBed.bed") @ By default all ROIs are split into bins of 100 bp with an overlap of 50 bp. The parameter {\tt limit} controls the minimum size of the ROIs that should be split (default = 0). The parameters {\tt bin} and {\tt shift} control the size of the bins and the no. of bp between start positions of adjacent bins. \section{runPanelcnMops} \label{s:panelcn.mops} The actual copy number analysis is done with the function \verb+runPanelcnMops+. The function requires a \verb+GRanges+ object of the RCs of test and control samples as well as the \verb+countWindows+ object used to extract these RCs. Optional parameters include a vector that indicates which samples to regard as test samples (default = c(1)), a vector of the names of the genes of interest (by default all genes are of interest), parameters for normalizing the RCs, a vector of expected fold changes for the copy number classes and a minimal median RC over all samples to exclude low coverage ROIs. <>= selectedGenes <- "ATM" XandCB <- test elementMetadata(XandCB) <- cbind(elementMetadata(XandCB), elementMetadata(control)) resultlist <- runPanelcnMops(XandCB, countWindows = countWindows, selectedGenes = selectedGenes) @ \section{Results} \label{s:results} The function \verb+runPanelcnMops+ returns a list of objects of the S4 class CNVDetectionResult, one CNVDetectionResult object per test sample. The structure of the CNVDetectionResult object can be listed by calling <>= (str(resultlist[[1]])) @ To get detailed information on which data are stored in such objects, enter <>= help(CNVDetectionResult) @ The CNVs per individual are stored in the slot \verb+integerCopyNumber+: <<>>= integerCopyNumber(resultlist[[1]])[1:5] @ The function \verb+createResultTable+ summarizes all relevant information for user selected genes of interest in a list of tables with one table per test sample: <<>>= sampleNames <- colnames(elementMetadata(test)) resulttable <- createResultTable(resultlist = resultlist, XandCB = XandCB, countWindows = countWindows, selectedGenes = selectedGenes, sampleNames = sampleNames) (tail(resulttable[[1]])) @ The table contains one line per Region Of Interest (ROI) with information about the RCs of the test sample ("RC"), the median RCs of all control samples ("medRC"), the normalized RCs of the test sample ("RC.norm"), the median of the normalized RCs of all control samples ("medRC.norm"), as well as the estimated CN ("CN"). Additionally, in the column "lowQual" low quality ROIs are flagged. \section{Visualization of results} \label{s:plot} \panelcnmops\ contains a plotting function that visualizes the normalized RCs of the samples analyzed as boxplots: \begin{center} <>= plotBoxplot(result = resultlist[[1]], sampleName = sampleNames[1], countWindows = countWindows, selectedGenes = selectedGenes, showGene = 1) @ \begin{figure}[H] \begin{center} \includegraphics[angle=0,width= 0.9\columnwidth]{001.pdf} \end{center} \end{figure} \end{center} The function expects a single CNVDetectionResult object as input together with the name of the test sample, the countWindows used, as well as a vector with the names of the genes of interest and an integer specifying which of the genes of interest to plot. \section{Analysis of chromosome X} \label{s:chrX} The analysis of ROIs on chromosome X is only possible if all samples have the same sex and the parameter sex of the function \verb+runPanelcnMops+ is set accordingly. The default "mixed" results in the removal of all X-chromosomal ROIs. Note, that if all samples are males CN2 in the results really corresponds to CN1. \section{Quality control} The panelcn.MOPS algorithm includes different quality control metrics. 1) ROIs are excluded if their median read count (RC) across all samples does not exceed a user defined threshold (default: 30), additionally a warning message is displayed. 2) ROIs are marked as "low quality" in the result table if their RCs show a high variation across all samples. 3) Samples with a median RC across all ROIs lower than 0.55 times the median of all samples are considered as low quality. 4) For each ROI the ratio between the normalized RCs of each sample compared to the median across all samples is calculated. Samples that show a high variation in these RC ratios are also flagged as low quality. Low quality samples are excluded if they are control samples which leads to a warning message. If a test sample is of low quality, only a warning message is displayed. \section{Adjusting sensitivity and specificity} The default parameters of the \panelcnmops\ algorithm were optimized on a data set of targeted NGS panel data with the aim of detecting CNVs ranging in size from part of a ROI to whole genes. However, you might want to adjust sensitivity and specificity to your specific needs. The parameter that influences sensitivity and specificity the most is {\tt I}, the vector of expected fold changes of the copy number classes. The default for \panelcnmops\, c(0.025, 0.57, 1, 1.46, 2), leads to a higher sensitivity compared to the default of \cnmops\ which is c(0.025, 0.5, 1, 1.5, 2). Increasing the values for CN0 and CN1 further and decreasing the values for CN3 and CN4 may help to improve the sensitivity, a change in the other direction may increase the specificity. Additional parameters that can be tuned to improve the results are the different normalization parameters: {\tt normType}, {\tt sizeFactor}, {\tt qu}, {\tt quSizeFactor}, and {\tt norm}. % \clearpage \section{How to cite this package} If you use this package for research that is published later, you are kindly asked to cite it as follows: \citep{Povysil:17}. To obtain Bib\TeX\ entries of the reference, you can enter the following into your R session: <>= toBibtex(citation("panelcn.mops")) @ \bibliographystyle{natbib} \bibliography{cnv} \end{document}