Package: CNVMetrics
Authors: Astrid Deschênes [aut, cre] (https://orcid.org/0000-0001-7846-6749), Pascal Belleau [aut] (https://orcid.org/0000-0002-0802-1071), David A. Tuveson [aut], Alexander Krasnitz [aut]
Version: 1.0.0
Compiled date: 2022-04-26
License: Artistic-2.0

1 Licensing

The CNVMetrics package and the underlying CNVMetrics code are distributed under the Artistic license 2.0. You are free to use and redistribute this software.



2 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



3 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 (Hiroi et al. 2013) and especially cancer (Stratton, Campbell, and Futreal 2009).

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 (Zhao et al. 2013). As an example, the Figure 1 shows the copy number profiles from sequencing data of two mouse pancreatic organoids (Oni et al. 2020), calculated with CNprep (Belleau et al. 2020) and plot with gtrellis (Gu, Eils, and Schlesner 2016).

Copy number profiles of two mouse metastatic pancreatic organoids (M10 and M30).

Copy number profiles of two mouse metastatic pancreatic organoids (M10 and M30).

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 (Gendoo et al. 2019).

The 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.



4 Installation

To install this package from Bioconductor, start R (version “4.2”) and enter:



5 Workflow for metrics calculated using CNV status calls

The following workflow gives an overview of the capabilities of 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.


5.1 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:

  • AMPLIFICATION
  • DELETION

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:

  • LOH
  • loh
  • 1-copy gain
  • GAIN
  • loss
  • and so on…

Beware that states with different spelling or upper/lower case nomenclature are considered as distinct states and are analyzed separately.

Example of a copy number file containing CNV calls.

Example of a copy number file containing CNV calls.


5.2 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 CNprep (Belleau et al. 2020), from ten mouse pancreatic organoids (Oni et al. 2020).

##    ID  chr     start       end    log2ratio         state
## 1 mT3 chr1   3527630 195471971  0.143253158       NEUTRAL
## 2 mT3 chr2   3549740 182113224  0.008638116       NEUTRAL
## 3 mT3 chr3   3536051 160039680  0.032630258       NEUTRAL
## 4 mT3 chr4   3734693 156508116 -0.026741963       NEUTRAL
## 5 mT3 chr5   3575586 111373751 -0.018519856       NEUTRAL
## 6 mT3 chr5 111373751 114262563  0.226961493 AMPLIFICATION
##  [1] "mT3"  "mM6"  "mM10" "mM15" "mT6"  "mT19" "mT23" "mT34" "mM1"  "mM3P"
## GRangesList object of length 10:
## $mM1
## GRanges object with 33 ranges and 2 metadata columns:
##        seqnames             ranges strand |  log2ratio         state
##           <Rle>          <IRanges>  <Rle> |  <numeric>   <character>
##    [1]     chr1   3527630-62533399      * |  0.0351701       NEUTRAL
##    [2]     chr1  62533399-82082589      * | -0.1066207       NEUTRAL
##    [3]     chr1 82082589-195471971      * |  0.0807449       NEUTRAL
##    [4]     chr2  3549740-182113224      * |  0.0527007       NEUTRAL
##    [5]     chr3  3536051-136509203      * |  0.3431661 AMPLIFICATION
##    ...      ...                ...    ... .        ...           ...
##   [29]    chr18   3558909-90702639      * |  0.0736673       NEUTRAL
##   [30]    chr19   3614043-61431566      * |  0.0902238       NEUTRAL
##   [31]     chrX   5577384-62973719      * | -0.2259987       NEUTRAL
##   [32]     chrX  62973719-84082850      * | -0.0641819       NEUTRAL
##   [33]     chrX 84082850-169538352      * | -0.3902963      DELETION
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths
## 
## ...
## <9 more elements>


5.3 Metric Calculation

The calculation of the similarity metrics is done with the calculateOverlapMetric() function.

## CNV Metrics
## Metric:
## sorensen
## 
## AMPLIFICATION:
##            mM1      mM10 mM15      mM3P mM6
## mM1         NA        NA   NA        NA  NA
## mM10 0.4437132        NA   NA        NA  NA
## mM15 0.0000000 0.0000000   NA        NA  NA
## mM3P 0.1713638 0.2380730    0        NA  NA
## mM6  0.4048136 0.5545147    0 0.2952499  NA
## [ -- omitted 5 rows/columns ]
## 
## DELETION:
##            mM1      mM10      mM15       mM3P mM6
## mM1         NA        NA        NA         NA  NA
## mM10 0.0000000        NA        NA         NA  NA
## mM15 0.4595765 0.0000000        NA         NA  NA
## mM3P 0.2568018 0.3336979 0.6481145         NA  NA
## mM6  0.0000000 0.0000000 0.0000000 0.00918594  NA
## [ -- omitted 5 rows/columns ]


5.4 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.



The plotMetric() function uses the pheatmap package to generate the graph. All arguments accepted by pheatmap::pheatmap() function are valid arguments.



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.

## [1] TRUE



6 Metrics using the CNV status calls

This survey represents the overlap metrics that are implemented in 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.


6.1 Sørensen

The Sørensen coefficient (Sørensen 1948) 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.


6.2 Szymkiewicz–Simpson

The Szymkiewicz–Simpson coefficient (Vijaymeena and K. 2016), 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.


6.3 Jaccard

The Jaccard coefficient (Jaccard 1912), 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.



7 Workflow for metrics calculated using the level of amplification/deletion

The following section gives an overview of the capabilities of CNVMetrics to calculate metrics using the the level of amplification/deletion (log2 ratio values):

The key functions for each step of the workflow are:

Step Function
Data Importation GenomicRanges::makeGRangesListFromDataFrame()
Metric Calculation calculateLog2ratioMetric()
Metric Visualization plotMetric()

The package::function() notation is used for functions from other packages.


7.1 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.

Example of a copy number file containing log2ratio values.

Example of a copy number file containing log2ratio values.


7.2 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 CNprep (Belleau et al. 2020), from ten mouse pancreatic organoids (Oni et al. 2020).

##    ID  chr     start       end    log2ratio         state
## 1 mT3 chr1   3527630 195471971  0.143253158       NEUTRAL
## 2 mT3 chr2   3549740 182113224  0.008638116       NEUTRAL
## 3 mT3 chr3   3536051 160039680  0.032630258       NEUTRAL
## 4 mT3 chr4   3734693 156508116 -0.026741963       NEUTRAL
## 5 mT3 chr5   3575586 111373751 -0.018519856       NEUTRAL
## 6 mT3 chr5 111373751 114262563  0.226961493 AMPLIFICATION
##  [1] "mT3"  "mM6"  "mM10" "mM15" "mT6"  "mT19" "mT23" "mT34" "mM1"  "mM3P"
## GRangesList object of length 10:
## $mM1
## GRanges object with 33 ranges and 2 metadata columns:
##        seqnames             ranges strand |  log2ratio         state
##           <Rle>          <IRanges>  <Rle> |  <numeric>   <character>
##    [1]     chr1   3527630-62533399      * |  0.0351701       NEUTRAL
##    [2]     chr1  62533399-82082589      * | -0.1066207       NEUTRAL
##    [3]     chr1 82082589-195471971      * |  0.0807449       NEUTRAL
##    [4]     chr2  3549740-182113224      * |  0.0527007       NEUTRAL
##    [5]     chr3  3536051-136509203      * |  0.3431661 AMPLIFICATION
##    ...      ...                ...    ... .        ...           ...
##   [29]    chr18   3558909-90702639      * |  0.0736673       NEUTRAL
##   [30]    chr19   3614043-61431566      * |  0.0902238       NEUTRAL
##   [31]     chrX   5577384-62973719      * | -0.2259987       NEUTRAL
##   [32]     chrX  62973719-84082850      * | -0.0641819       NEUTRAL
##   [33]     chrX 84082850-169538352      * | -0.3902963      DELETION
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths
## 
## ...
## <9 more elements>


7.3 Metric Calculation

The calculation of the similarity metrics is done with the calculateOverlapMetric() function.

## CNV Metrics
## Metric:
## weightedEuclideanDistance
## 
## LOG2RATIO:
##             mM1       mM10       mM15       mM3P mM6
## mM1          NA         NA         NA         NA  NA
## mM10 0.12503352         NA         NA         NA  NA
## mM15 0.12100344 0.08125775         NA         NA  NA
## mM3P 0.09499932 0.06838924 0.09040692         NA  NA
## mM6  0.10819206 0.08071789 0.07608047 0.07192307  NA
## [ -- omitted 5 rows/columns ]


7.4 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.


The plotMetric() function uses the pheatmap package to generate the graph. All arguments accepted by pheatmap::pheatmap function are valid arguments.



8 Metrics using the level of amplification/deletion

This section presents the similarity measure that is implemented in 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.


8.1 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.



9 Supplementary information


9.1 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 BiocParallel package (Morgan et al. 2021).

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 (Seshan and Olshen 2021) 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 (Kusnierczyk 2012) with 100 replications. This comparison has been done on a high performance computing (HPC) server:

Number of threads
(nJobs parameter)		 Average elapsed time
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


9.3 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.



10 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.



11 Session info

Here is the output of sessionInfo() on the system on which this document was compiled:

## R version 4.2.0 RC (2022-04-19 r82224)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] CNVMetrics_1.0.0     GenomicRanges_1.48.0 GenomeInfoDb_1.32.0 
## [4] IRanges_2.30.0       S4Vectors_0.34.0     BiocGenerics_0.42.0 
## [7] knitr_1.38           BiocStyle_2.24.0    
## 
## loaded via a namespace (and not attached):
##  [1] highr_0.9              bslib_0.3.1            compiler_4.2.0        
##  [4] BiocManager_1.30.17    RColorBrewer_1.1-3     jquerylib_0.1.4       
##  [7] XVector_0.36.0         bitops_1.0-7           tools_4.2.0           
## [10] zlibbioc_1.42.0        digest_0.6.29          lifecycle_1.0.1       
## [13] jsonlite_1.8.0         evaluate_0.15          gtable_0.3.0          
## [16] rlang_1.0.2            cli_3.3.0              yaml_2.3.5            
## [19] parallel_4.2.0         xfun_0.30              fastmap_1.1.0         
## [22] GenomeInfoDbData_1.2.8 gridExtra_2.3          stringr_1.4.0         
## [25] sass_0.4.1             grid_4.2.0             R6_2.5.1              
## [28] BiocParallel_1.30.0    rmarkdown_2.14         pheatmap_1.0.12       
## [31] farver_2.1.0           magrittr_2.0.3         scales_1.2.0          
## [34] htmltools_0.5.2        colorspace_2.0-3       stringi_1.7.6         
## [37] munsell_0.5.0          RCurl_1.98-1.6



References

Belleau, P, A Deschênes, G Sun, David A Tuveson, and Alex Krasnitz. 2020. “CNprep: Copy number event detection [version 1; not peer reviewed].” F1000Research 9 (743): (poster). https://doi.org/10.7490/f1000research.1118065.1.

Gendoo, Deena M. A., Robert E. Denroche, Amy Zhang, Nikolina Radulovich, Gun Ho Jang, Mathieu Lemire, Sandra Fischer, et al. 2019. “Whole genomes define concordance of matched primary, xenograft, and organoid models of pancreas cancer.” Edited by Amos Tanay. PLOS Computational Biology 15 (1): e1006596. https://doi.org/10.1371/journal.pcbi.1006596.

Gu, Zuguang, Roland Eils, and Matthias Schlesner. 2016. “gtrellis: an R/Bioconductor package for making genome-level Trellis graphics.” BMC Bioinformatics 17 (1): 169. https://doi.org/10.1186/s12859-016-1051-4.

Hiroi, N, T Takahashi, A Hishimoto, T Izumi, S Boku, and T Hiramoto. 2013. “Copy number variation at 22q11.2: from rare variants to common mechanisms of developmental neuropsychiatric disorders.” Molecular Psychiatry 18 (11): 1153–65. https://doi.org/10.1038/mp.2013.92.

Jaccard, Paul. 1912. “The distribution of the flora in the alpine zone.” New Phytologist 11 (2): 37–50. https://doi.org/10.1111/j.1469-8137.1912.tb05611.x.

Kusnierczyk, Wacek. 2012. Rbenchmark: Benchmarking Routine for R. https://CRAN.R-project.org/package=rbenchmark.

Morgan, Martin, Valerie Obenchain, Michel Lang, Ryan Thompson, and Nitesh Turaga. 2021. BiocParallel: Bioconductor Facilities for Parallel Evaluation. https://github.com/Bioconductor/BiocParallel.

Oni, Tobiloba E, Giulia Biffi, Lindsey A Baker, Yuan Hao, Claudia Tonelli, Tim D D Somerville, Astrid Deschênes, et al. 2020. “SOAT1 promotes mevalonate pathway dependency in pancreatic cancer.” Journal of Experimental Medicine 217 (9). https://doi.org/10.1084/jem.20192389.

Seshan, Venkatraman E., and Adam Olshen. 2021. DNAcopy: DNA Copy Number Data Analysis.

Sørensen, Thorvald. 1948. “A method of establishing groups of equal amplitude in plant sociology based on similarity of species and its application to analyses of the vegetation on Danish commons.” Biologiske Skrifter 29 (5): 1–34.

Stratton, Michael R, Peter J Campbell, and P Andrew Futreal. 2009. “The cancer genome.” Nature 458 (7239): 719–24. https://doi.org/10.1038/nature07943.

Vijaymeena, M.K, and Kavitha K. 2016. “A Survey on Similarity Measures in Text Mining.” Machine Learning and Applications: An International Journal 3 (1): 19–28. https://doi.org/10.5121/mlaij.2016.3103.

Zhao, Min, Qingguo Wang, Quan Wang, Peilin Jia, and Zhongming Zhao. 2013. “Computational tools for copy number variation (CNV) detection using next-generation sequencing data: features and perspectives.” BMC Bioinformatics 14 (S11): S1. https://doi.org/10.1186/1471-2105-14-S11-S1.