---
title: "Quantifying similarity between copy number profiles"
author: Astrid Deschênes, Pascal Belleau and Alexander Krasnitz
output:
rmarkdown::html_document:
highlight: pygments
toc: true
fig_width: 5
number_sections: true
pkgdown:
asis: true
bibliography: CNVMetrics.bibtex
vignette: >
%\VignetteIndexEntry{Copy number variant metrics}
%\VignettePackage{CNVMetrics}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r style, echo = FALSE, results = 'asis', warning=FALSE, message=FALSE}
BiocStyle::markdown()
suppressPackageStartupMessages({
library(knitr)
library(GenomicRanges)
library(CNVMetrics)
})
set.seed(121444)
```
<br />
**Package**: `r Rpackage("CNVMetrics")`<br />
**Authors**: `r packageDescription("CNVMetrics")[["Author"]]`<br />
**Version**: `r packageDescription("CNVMetrics")$Version`<br />
**Compiled date**: `r Sys.Date()`<br />
**License**: `r packageDescription("CNVMetrics")[["License"]]`<br />
# Licensing
The `r Githubpkg("KrasnitzLab/CNVMetrics")` package and the underlying
`r Githubpkg("KrasnitzLab/CNVMetrics")` code are distributed under the
Artistic license 2.0. You are free to use and redistribute this software.
<br>
<br>
# Citing
If you use this package for a publication, we would ask you
to cite the following:
>Belleau P, Deschênes A, Beyaz S et al. CNVMetrics package: Quantifying similarity between copy number profiles [version 1; not peer reviewed]. F1000Research 2021, 10:737 (slides) (doi: 10.7490/f1000research.1118704.1)
[F1000Research poster](http://www.doi.org/10.7490/f1000research.1118704.1)
<br>
<br>
# Introduction
Copy number variation (CNV) includes multiplication and deletion of DNA
segment. Copy number variations have been shown to be associated
with a wide spectrum of pathological conditions and complex traits, such
as developmental neuropsychiatric disorders [@Hiroi2013] and especially
cancer [@Stratton2009].
CNVs are usually reported, for each sample, as genomic regions that are
duplicated or deleted with respect to a reference. Those regions are denoted
as _CNV status calls_. The level of amplification or deletion can also be
reported, usually in log2 ratio values or normalized read depth [@Zhao2013].
As an example, the Figure 1 shows the copy number profiles from sequencing
data of two mouse pancreatic organoids [@Oni2020], calculated with
`r Githubpkg("KrasnitzLab/CNprep")` [@Belleau2020] and plot with
`r Biocpkg("gtrellis")` [@Gu2016a].
```{r graphCNVpro, echo = FALSE, fig.align="center", fig.cap="Copy number profiles of two mouse metastatic pancreatic organoids (M10 and M30).", out.width = '95%', results='asis'}
knitr::include_graphics("CNV_mM30_mM10_v03_Feb_08_2021_small.jpeg")
```
While visual representation is a practical way to qualitatively compare copy
number profiles, metrics are useful statistical tools for quantitatively
measuring similarity and dissimilarity between profiles. Similarity metrics
can be employed to compare CNV profiles of genetically unrelated samples.
Moreover, those metrics can as well be put to use on samples with common
genetic background. As an example, a comparison between primary and metastatic
tumor CNV profiles may reveal genomic determinants of metastasis. Similarly,
patient-derived xenograft or organoid models of cancer are expected to
recapitulate CNV patterns of the tumor tissue of origin [@Gendoo2019].
The `r Githubpkg("KrasnitzLab/CNVMetrics")` package calculates metrics to
estimate the level of similarity between copy number profiles. Some metrics
are calculated using the _CNV status calls_ (amplification/deletion/LOH status
or any user specific status) while others are based on the level of
amplification/deletion in log2 ratio. Finally, a visualization tool is
provided to explore resulting metrics in the form of
sample-to-sample heatmaps.
<br>
<br>
# Installation
To install this package from [Bioconductor](https://bioconductor.org/packages/CNVMetrics), start R (version "4.2") and enter:
```{r installDemo01, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
# The following initializes usage of Bioc devel
BiocManager::install(version='devel')
BiocManager::install("CNVMetrics")
```
<br>
<br>
# Workflow for metrics calculated using CNV status calls
The following workflow gives an overview of the capabilities of
`r Githubpkg("KrasnitzLab/CNVMetrics")` to calculate metrics using the
_CNV status calls_ (amplification/deletion status or any user specific
status):
The key functions for each step of the workflow are:
Step | Function
----------------------- | ---------------------------------------------
Data Importation | `GenomicRanges::makeGRangesListFromDataFrame()`
Metric Calculation | `calculateOverlapMetric()`
Metric Visualization | `plotMetric()`
The `package::function()` notation is used for functions from other packages.
<br>
## Data Input - Copy number file containing the CNV status calls
_CNV status calls_ are represented as segments with a copy number state. The
state be general, such as "amplification", "deletion" or "neutral", or more
specific such as of loss of heterozygosity (LOH), 1-copy gain, 2-copy gain,
1-copy loss and so on.
A basic five-column input file containing genomic position
(chromosome, start, end), sample identification and _CNV status calls_ is
required. All samples that need to be analyzed together have to be combined
into one file.
A column named **state** is required. In this column, The _CNV status call_ of
each segment must be specified using a string. By default, the states that are
analyzed by this package are the amplification/deletion states with
this specific notation:
* **<span style="color:darkred">AMPLIFICATION</span>**
* **<span style="color:darkred">DELETION</span>**
Segments with other **state** values can be present in the file. However,
those segments won't be retain for the calculation of the metrics.
However, the user can define is how notation and decided which **state** will
be used to calculate the similarity metrics. The user defined states can be
in upper or lower cases. Examples of possible states:
* **<span style="color:darkred">LOH</span>**
* **<span style="color:darkred">loh</span>**
* **<span style="color:darkred">1-copy gain</span>**
* **<span style="color:darkred">GAIN</span>**
* **<span style="color:darkred">loss</span>**
* and so on...
Beware that states with different spelling or upper/lower case nomenclature
are considered as distinct states and are analyzed separately.
```{r figureCNVFile, echo = FALSE, fig.align="center", fig.cap="Example of a copy number file containing CNV calls.", out.width = '100%'}
knitr::include_graphics("Input_CNV_call_300ppi_v02_low_quality.jpg")
```
<br>
## Data Importation - GRangesList
The input format for the copy number information, as needed by the
`calculateOverlapMetric()` function, is a `GRangesList` object.
The easiest way to generate a `GRangesList` object is to first load the
copy number information into an R `data.frame` and then, use the
`GenomicRanges::makeGRangesListFromDataFrame()` function to convert them
to a `GRangesList`.
For this demonstration, we consider _CNV status calls_ as obtained with
`r Githubpkg("KrasnitzLab/CNprep")` [@Belleau2020],
from ten mouse pancreatic organoids [@Oni2020].
```{r demoImport01}
## Load required libraries
library(GenomicRanges)
library(CNVMetrics)
## Load file containing CNV calls for 10 mouse organoids
data.dir <- system.file("extdata", package="CNVMetrics")
cnv.file <- file.path(data.dir, "mousePairedOrganoids.txt")
calls <- read.table(cnv.file, header=TRUE, sep="\t")
## The CNV status calls for all samples are present in one file
## The 'state' column is required
## The chromosome Y has been removed
head(calls)
## The ID column identifies the 10 samples
unique(calls[,"ID"])
## The ID column is used to split the samples into different GRanges
## inside a GRangesList
## The 'keep.extra.column=TRUE' parameter is needed to retained the extra
## column 'state' that is needed for the calculation of the metrics
grl <- GenomicRanges::makeGRangesListFromDataFrame(calls,
split.field="ID", keep.extra.columns=TRUE)
grl
```
<br>
## Metric Calculation
The calculation of the similarity metrics is done with the
`calculateOverlapMetric()` function.
```{r demoCalculateMetric01}
## In this case, the default states (AMPLIFICATION, DELETION) are used.
## So, the 'states' parameter doesn't have to be specified
## The 'states' parameter needs to be adjusted for user-specific states
## Ex: states=c("LOH", "gain")
metric <- calculateOverlapMetric(segmentData=grl, method="sorensen", nJobs=1)
metric
```
<br>
## Metric Visualization
A heatmap of this similarity metrics can be a useful tool to get an overview
over similarities and dissimilarities between samples.
The `plotMetric()` function generates a graphical representation of
the similarity metrics in the form of a sample-to-sample heatmap. By default,
an hierarchical clustering based on the sample distances
(1-metric) is used. When NA values are present in the metric matrix, those
are replaced by zero.
```{r demoPlot01, fig.align="center", fig.height=4.5, fig.width=4.5}
## Create graph for the metrics related to amplified regions
plotMetric(metric, type="AMPLIFICATION")
```
<br>
<br>
The `plotMetric()` function uses the `r CRANpkg("pheatmap")` package to
generate the graph. All arguments accepted by `pheatmap::pheatmap()` function
are valid arguments.
```{r demoPlot02, fig.align="center", fig.height=4.8, fig.width=4.8}
## Create graph for the metrics related to deleted regions
## Metric values are printed as 'display_numbers' and 'number_format' are
## arguments recognized by pheatmap() function
plotMetric(metric, type="DELETION",
colorRange=c("white", "darkorange"),
show_colnames=TRUE,
display_numbers=TRUE,
number_format="%.2f")
```
<br>
<br>
Row and/or column annotation is often useful and can easily be done
by using the `annotation_row` or `annotation_col` arguments, as described in
the `pheatmap::pheatmap` method.
```{r demoPlot03, fig.align="center", fig.height=4.8, fig.width=6}
## Load file containing annotations for the mouse organoids
## The mouse ID identifying the source of the sample
## The stage identifying the state (tumor vs metastasis) of the sample
data.dir <- system.file("extdata", package="CNVMetrics")
annotation.file <- file.path(data.dir, "mousePairedOrganoidsInfo.txt")
annotOrg <- read.table(annotation.file, header=TRUE, sep="\t")
## The row names must correspond to the names assigned to the rows/columns
## in the CNVMetric object
rownames(annotOrg) <- annotOrg$ID
annotOrg$ID <- NULL
all(rownames(annotOrg) == rownames(metric$AMPLIFICATION))
## Create graph for the metrics related to amplified regions
## Rows are annotated with the stage and mouse information
plotMetric(metric, type="AMPLIFICATION",
colorRange=c("white", "steelblue"),
annotation_row=annotOrg)
```
<br>
<br>
# Metrics using the _CNV status calls_
This survey represents the overlap metrics that are implemented in
`r Githubpkg("KrasnitzLab/CNVMetrics")` package. Those metrics are calculated
using the _CNV status calls_. The size of the amplified/deleted regions as
well as the size of the overlapping of regions are always in base paired.
<br>
## Sørensen
The Sørensen coefficient [@Sorensen48] is calculated by dividing twice the
size of the intersection by the sum of the size of the two sets:
\begin{equation}
\frac{2\times \left| X \cap Y \right| }{\left| X \right| + \left| Y \right|}
(\#eq:sorensen)
\end{equation}
where $X$ and $Y$ represent the regions of each sample in base paired.
<br>
## Szymkiewicz–Simpson
The Szymkiewicz–Simpson coefficient [@Vijaymeena2016], also known as the
overlap coefficient, is calculated by dividing the size of the intersection
by the smaller of the size of the two sets:
\begin{equation}
\frac{\left| X \cap Y \right|}{min \left(\left| X \right|,\left| Y \right|\right)}
(\#eq:szymkiewicz)
\end{equation}
where $X$ and $Y$ represent the regions of each sample in base paired. If
set $X$ is a subset of $Y$ or vice versa, the overlap coefficient
value is 1.
<br>
## Jaccard
The Jaccard coefficient [@Jaccard1912], also known as coefficient of
community, is calculated by dividing the size
of the intersection by the smaller of the size of the two sets:
\begin{equation}
\frac{\left| X \cap Y \right| }{ \left| X \cup Y \right|}
(\#eq:jaccard)
\end{equation}
where $X$ and $Y$ represent the regions of each sample in base paired.
<br>
<br>
# Workflow for metrics calculated using the level of amplification/deletion
The following section gives an overview of the capabilities of
`r Githubpkg("KrasnitzLab/CNVMetrics")` to calculate metrics using the
_the level of amplification/deletion_ (log2 ratio values):
The key functions for each step of the workflow are:
<div align="center">
<div style="width: 60%; margin:0 left;text-align: center;">
Step | Function
----------------------- | ---------------------------------------------
Data Importation | `GenomicRanges::makeGRangesListFromDataFrame()`
Metric Calculation | `calculateLog2ratioMetric()`
Metric Visualization | `plotMetric()`
</div>
</div>
The `package::function()` notation is used for functions from other packages.
<br>
## Data Input - Copy number file containing the level of amplification/deletion
Copy number are often represented as segments with a copy number state and/or
the level of amplification/deletion. One usual unit to quantify the level
of amplification or deletion is in log2 ratio.
A basic five-column input file containing genomic position
(chromosome, start, end), sample identification and
_the level of amplification/deletion_ is
required. All samples that need to be analyzed together have to be combined
into one file.
A column named **log2ratio** is required. In this column, the amplified and
deleted segments must be assigned a numerical value representing the log2ratio
or NA.
```{r figureCNVFileLog2Ratio, echo = FALSE, fig.align="center", fig.cap="Example of a copy number file containing log2ratio values.", out.width = '100%'}
knitr::include_graphics("Input_CNV_log2ratio_v01_low_quality.jpg")
```
<br>
## Data Importation - GRangesList
The input format for the copy number information, as needed by the
`calculateLog2ratioMetric()` function, is a `GRangesList` object.
The easiest way to generate a `GRangesList` object is to first load the
copy number information into an R `data.frame` and then, use the
`GenomicRanges::makeGRangesListFromDataFrame()` function to convert them
to a `GRangesList`.
For this demonstration, we consider _the level of amplification/deletion_ as
obtained with `r Githubpkg("KrasnitzLab/CNprep")` [@Belleau2020],
from ten mouse pancreatic organoids [@Oni2020].
```{r demoImport02}
## Load required libraries
library(GenomicRanges)
library(CNVMetrics)
## Load file containing CNV calls for 10 mouse organoids
data.dir <- system.file("extdata", package="CNVMetrics")
cnv.file <- file.path(data.dir, "mousePairedOrganoids.txt")
calls <- read.table(cnv.file, header=TRUE, sep="\t")
## The CNV status calls for all samples are present in one file
## The 'log2ratio' column is required
## The chromosome Y has been removed
head(calls)
## The ID column identifies the 10 samples
unique(calls[,"ID"])
## The ID column is used to split the samples into different GRanges
## inside a GRangesList
## The 'keep.extra.column=TRUE' parameter is needed to retained the extra
## column 'state' that is needed for the calculation of the metrics
grlog <- GenomicRanges::makeGRangesListFromDataFrame(df=calls,
split.field="ID", keep.extra.columns=TRUE)
grlog
```
<br>
## Metric Calculation
The calculation of the similarity metrics is done with the
`calculateOverlapMetric()` function.
```{r demoCalculateMetricLog}
metricLog <- calculateLog2ratioMetric(segmentData=grlog,
method="weightedEuclideanDistance", nJobs=1)
metricLog
```
<br>
## Metric Visualization
A heatmap of this similarity metrics can be a useful tool to get an overview
over similarities and dissimilarities between samples.
The `plotMetric()` function generates a graphical representation of
the similarity metrics in the form of a sample-to-sample heatmap. By default,
an hierarchical clustering based on the sample distances
(1-metric) is used. When NA values are present in the metric matrix, those
are replaced by zero.
```{r demoPlotLog01, fig.align="center", fig.height=4.5, fig.width=4.5}
## Create graph for the metrics related to weighted Euclidean distance-based
plotMetric(metricLog)
```
<br>
The `plotMetric()` function uses the `r CRANpkg("pheatmap")` package to
generate the graph. All arguments accepted by `pheatmap::pheatmap` function
are valid arguments.
```{r demoPlotLog02, fig.align="center", fig.height=4.8, fig.width=4.8}
## Create graph for the weighted Euclidean distance-based metrics
## Remove title (argument main="")
## Metric values are printed as 'display_numbers' and 'number_format' are
## arguments recognized by pheatmap() function
plotMetric(metricLog, colorRange=c("white", "darkorange"),
show_colnames=TRUE,
display_numbers=TRUE,
number_format="%.2f",
main="")
```
<br>
<br>
# Metrics using the level of amplification/deletion
This section presents the similarity measure that is implemented in
`r Githubpkg("KrasnitzLab/CNVMetrics")` package. This metric are calculated
using the level of amplification/deletion. The level of
amplification/deletion is in log2 ratio while the size of the regions
is in base paired.
<br>
## Weighted Euclidean Distance-Based
The Weighted Euclidean Distance corresponds to the euclidean distance between
the log2 values of the two samples multiplied by the natural logarithm
of the number of bases of the analyzed segments. The final metric is 1 over
1 added to the squared sum of the values obtained for all segments included
in the calculation.
The Weighted Euclidean Distance corresponds to the euclidean distance between
the log2 values of the two samples multiplied by the natural logarithm
of the number of bases of the analyzed segments. The final metric is 1 over
1 added to the squared sum of the values obtained for all segments included
in the calculation.
\begin{equation}
\frac{1}{1 + \sqrt{\sum_{i=1} log_{2}(w_{i}) (A_{i} - B_{i})^{2}}}
(\#eq:euclidean)
\end{equation}
where $A_{i}$ and $B_{i}$ represent the log2 ratio values of samples $A$ and
$B$ for the region $i$ while $w_{i}$ is the length of region $i$ in base
paired.
<br>
<br>
# Supplementary information
<br>
## Using parallelization
When the number of samples is limited, the above steps should be processed
in a few minutes. However, for datasets with a high number of samples, the
combinatorial calculation of the metrics can lead to longer processing time.
In this context, take advantage of parallelized computation is a viable
option. Both *calculateOverlapMetric()* and *calculateLog2ratioMetric()*
functions have paralleled implementation done with the
`r Biocpkg("BiocParallel")` package [@Morgan2021].
The copy number data from The Cancer Genome Atlas (TCGA) Uterine
Carcinosarcoma (UCS) study generated by the TCGA Research Network
(https://www.cancer.gov/tcga) is used as an demonstration. The copy number
variation information, as obtained from the DNACopy workflow [@DNAcopy]
is available for 53 patients.
The following table highlights the time differences for processing the
Sørensen metric for all samples (metrics for all the 1378 possible
combinations) using rbenchmark [@rbenchmark] with 100 replications. This
comparison has been done on a high performance computing (HPC) server:
<div align="center">
<div style="width: 60%; margin:0 left;text-align: center;">
| Number of threads<br>(*nJobs* parameter) | Average elapsed time<br> |
| -------- | ----------- |
| 24 | 1 min 26 sec |
| 16 | 1 min 27 sec |
| 8 | 1 min 51 sec |
| 4 | 3 min 15 sec |
| 1 | 7 min 37 sec |
</div>
</div>
<br>
## Creating your own GRangesList
The `GenomicRanges::makeGRangesListFromDataFrame()` function enables the
creation of a list of `GRangesList` objects from a `data.frame`. However,
`GRangesList` can also be generated and filled manually.
```{r demoGR}
## First, create the GRanges objects; one per sample
gr1 <- GRanges(seqnames="chr2", ranges=IRanges(3, 6000),
strand="+", state="AMPLIFICATION",
log2ratio=0.45)
gr2 <- GRanges(seqnames=c("chr1", "chr2"),
ranges=IRanges(c(7,5555), width=c(1200, 40)),
strand=c("+", "-"), state=c("NEUTRAL", "AMPLIFICATION"),
log2ratio=c(0.034, 0.5))
gr3 <- GRanges(seqnames=c("chr1", "chr2"),
ranges=IRanges(c(1, 5577), c(3, 5666)),
strand=c("-", "-"), state=c("NEUTRAL", "AMPLIFICATION"),
log2ratio=c(0.04, 0.31))
## Then, construct a GRangesList() using all the GRanges objects
grl <- GRangesList("sample01"=gr1, "sample02"=gr2, "sample03"=gr3)
```
<br>
## Reproducible research
To ensure reproducible results, *set.seed()* function should be call
before *calculateOverlapMetric()* and *calculateLog2ratioMetric()*.
Beware that the *nJobs* parameter must also be fixed; change in the value of
the *nJobs* parameter might lead to different results.
```{r demoSeed, eval=FALSE}
## First, fixe the seed value
set.seed(121234)
## Run the method to calculated the desired metrics
## The number of jobs (*nJobs* parameter) can be higher than one but
## have to remain the same then the calculation is redone to ensure
## reproducitble results
metricLog <- calculateLog2ratioMetric(segmentData=grlog,
method="weightedEuclideanDistance", nJobs=1)
```
<br>
<br>
# Acknowledgments
This work was supported by the Lustgarten Foundation, where David A. Tuveson
is a distinguished scholar and Director of the Lustgarten Foundation–designated
Laboratory of Pancreatic Cancer Research.
<br>
<br>
# Session info
Here is the output of `sessionInfo()` on the system on which this document was
compiled:
```{r sessionInfo, echo=FALSE}
sessionInfo()
```
<br>
<br>
# References