---
title: "Secondary analyses of CNV data (HRD and more) with oncoscanR"
author: "Yann Christinat"
date: "`r format(Sys.time(), '%B %d, %Y')`"
output:
    BiocStyle::html_document:
        number_sections: yes
        toc: true
vignette: >
    %\VignetteIndexEntry{oncoscanR vignette}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)
```

# OncoscanR package description
OncoscanR is an R package to handle Copy Number Variation analyses originating 
from the Oncoscan assay (Affymetrix). It allows computation of different 
homologous recombination deficiency (HRD) scores and the tandem duplication plus
score (TDplus) to identify CDK12-mutated tumors [Popova et al., 
Cancer Res 2016]. 
The package also allows for identification of arm-level alterations (e.g. gain 
of chromosome arm 1p).

## Inclusion in Bioconductor
The package allows secondary analysis of biological microarray data. It relies 
on existing bioconductor packages to enable rigourous and reproducible analyses.
The package could, for example, be used to compute an HRD score using CNV data 
from the `conumee` package. 

**IMPORTANT**: The package expects as input the text exported file from ChAS 
(Chromosome Analysis Suite; the Affymetrix software to identify CNV segments 
from the Oncoscan Assay). The package assumes that all segments given in the 
file are correct and true. The ChAS text file has to contain the columns `Type`,
`CN State` and `Full Location` (to setup in ChAS). Any text file that complies 
with this structure should work equally well.

Starting with version 1.3.0, ASCAT output files can also be used as input.

Note that the Oncoscan does not cover the p arms of chromosome 13, 14, 15 and 
22. 
The coverage on the p arm of chromosome 21 is only partial and is not included 
in the standard Oncoscan workflow (function `workflow_oncoscan.chas` or script 
`bin/oncoscan-workflow.R`).

# Getting started
## Installation

OncoscanR is available through Bioconductor. In R, enter the commands:

```{r eval=FALSE}
if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")

BiocManager::install("oncoscanR")
```

You can also install the development version via GitHub.:
```{r eval=FALSE}
# install.packages('devtools')
install_github('yannchristinat/oncoscanR')
```

Note that the installation of the development version requires the prior 
installation of the packages `GenomicRanges` (bioconductor), `magrittr`, 
`jsonlite` and `readr`.

```{r eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GenomicRanges")

install.packages(c("magrittr", "jsonlite", "readr"))
```

## Testing the installation
Open R and type the following commands:
```{r}
library(oncoscanR)
segs.filename <- system.file("extdata", "chas_example.txt", 
                             package = "oncoscanR")
res <- workflow_oncoscan.chas(segs.filename)
print(res)
```

If everything is setup fine, `res` should contain a list with no arm-level 
alterations and a negative HRD score (nLST=1).


# Use cases
## Loading a ChAS export file and do a bit of cleaning
ChAS exports files contain only basic information about the copy number (gain, 
loss or LOH), plus the segment may overlap the centromere. 
When the file is loaded by OncoscanR (`load_chas` function), all segments are 
assigend to a chromosomal arms and split if necessary.
The LOH segments are given by ChAS independently of the copy number variation 
segments. Therefore one may have a LOH segment overlapping with a copy loss. As 
this information is redundant (a copy loss will always have a LOH), we need to 
trim and split these LOH with the `adjust_loh` function.

```{r}
library(magrittr)

# Load the ChAS file
chas.fn <- system.file("extdata", "chas_example.txt", package = "oncoscanR")
segments <- load_chas(chas.fn, oncoscan_na33.cov)

# Clean the segments: restricted to Oncoscan coverage, LOH not overlapping
# with copy loss segments, smooth&merge segments within 300kb and prune
# segments smaller than 300kb.
segs.clean <- trim_to_coverage(segments, oncoscan_na33.cov) %>%
    adjust_loh() %>%
    merge_segments() %>%
    prune_by_size()
```

Of note, the `oncoscan_na33.cov` objects contains the genomic coverage of the 
oncoscan assay (start/end for each chromosomal arm, hg19). One could re-compute 
the latter by downloading the annotation file from the ThermoFisher website
and process it with the `get_oncoscan_coverage_from_probes`function.

A similar function is available for loading a result file from the ASCAT 
program: `load_ascat`. The ASCAT file is expected to have the following column 
names:
- 'chr' (chromosome number, with or withour "chr")
- 'startpos' (first position of CNV segment)
- 'endpos' (last position of CNV segment)
- 'nMajor' (Number of copies of the major allele)
- 'nMinor' (Number of copies of the minor allele)

```{r}
ascat.fn <- system.file("extdata", "ascat_example.txt", package = "oncoscanR")
ascat.segments <- load_ascat(ascat.fn, oncoscan_na33.cov)
head(ascat.segments)
```

## Computation of arm-level alteration
Function `armlevel_alt`

An arm is declared globally altered if more than 80% of its bases are altered 
with a similar CNV type (amplifications [3 extra copies or more], gains [1-2 
extra copies], losses or copy-neutral losses of heterozygozity [LOH]). For 
instance, "gain of 3p" indicates that there is more than 80% of arm with 3 
copies but less than 80% with 5 (otherwise it would be an amplification). Prior 
to computation, segments of same copy number and at a distance <300Kbp (Oncoscan
resolution genome-wide) are merged. The remaining segments are filtered to a 
minimum size of 300Kbp.

For instance if we want to get all arms that have a global LOH alteration, we 
run:
```{r}
chas.fn <- system.file("extdata", "triploide_gene_list_full_location.txt", 
                       package = "oncoscanR")
segments <- load_chas(chas.fn, oncoscan_na33.cov)
armlevel.loh <- get_loh_segments(segments) %>% 
    armlevel_alt(kit.coverage = oncoscan_na33.cov)
```

The variable `armlevel.loh`is a named vector containing the arms that have 
percentage of base with LOH above the threshold (90%). To obtain the percentage 
of LOH bases in all arms, one could set the threshold to zero:

```{r}
armlevel.loh <- get_loh_segments(segments) %>%
        armlevel_alt(kit.coverage = oncoscan_na33.cov, threshold = 0)
```


## Global level of alteration
Several functions are available to perform such computation:

- `score_avgcn`: compute the average copy number across the genome
- `score_estwgd`: computes an estimation of the number of whole-genome doubling 
events
- `score_mbalt`: computes the total number of Mbp that have an alteration (w/o 
LOH segments)

```{r}
mbalt <- score_mbalt(segments, oncoscan_na33.cov)
percent.alt <- mbalt['sample']/mbalt['kit']
message(paste(mbalt['sample'], 'Mbp altered ->', round(percent.alt*100), 
              '% of genome'))

avgcn <- score_avgcn(segments, oncoscan_na33.cov)
wgd <- score_estwgd(segments, oncoscan_na33.cov)
message(paste('Average copy number:', round(avgcn, 2), '->', wgd['WGD'], 
'whole-genome doubling event'))
```


## HRD scores
The package contains several HRD scores described below.

### Score LST 
Function `score_lst`

Procedure based on the paper from Popova et al, Can. Res. 2012 (PMID: 22933060).
First segments smaller than 3Mb are removed, then segments are smoothed with 
respect to copy number at a distance of 3Mb.
The number of LSTs is the number of breakpoints (breakpoints closer than 3Mb are
merged) that have a segment larger or equal to 10Mb on each side. This score was
linked to BRCA1/2-deficient tumors.

### Score HR-LOH 
Function `score_loh`

Procedure based on the paper from Abkevich et al., Br J Cancer 2012 (PMID: 
23047548). 
Number of LOH segments larger than 15Mb but excluding segments on chromosomes 
with a global LOH alteration. This score was linked to BRCA1/2-deficient tumors.

### Score nLST
Function `score_nlst`

HRD score developed at HUG and based on the LST score by Popova et al. but 
normalized by an estimation of the number of whole-genome doubling events.
Of note, copy-neutral LOH segments are removed before computation.

`nLST = LST - 7*W/2` where `W` is the number of whole-genome doubling events.

The score is positive if there are at least 15 nLST.

The nLST score has been validated on 469 high grade ovarian cancer samples from
the PAOLA-1 clinical trial and is used in routine at the Geneva University 
Hospitals for prediction of PARP inhibitors response.

*How to cite*

Christinat Y, Ho L, Clément S, et al. 2022-RA-567-ESGO The Geneva HRD test: 
clinical validation on 469 samples from the PAOLA-1 trial. International Journal
of Gynecologic Cancer 2022;32:A238-A239.


### Score gLOH
Function `score_gloh`

The percentage genomic LOH score is computed as described in the FoundationFocus
CDx BRCA LOH assay; i.e. the percentage of bases covered by the Oncoscan that 
display a loss of heterozygosity independently of the number of copies, 
excluding chromosomal arms that have a global LOH (>=90% of arm length). To 
compute with the armlevel_alt function on LOH segments only). This score was 
linked to BRCA1/2-deficient tumors.


### Example
First we need to load and clean the ChAS export file (from a female patient). We
adjust the  Oncoscan coverage to exclude the 21p arm as it is only partially 
covered.

```{r}
# Load data
chas.fn <- system.file("extdata", "LST_gene_list_full_location.txt", 
                       package = "oncoscanR")
segments <- load_chas(chas.fn, oncoscan_na33.cov)

# Clean the segments: restricted to Oncoscan coverage, LOH not overlapping
# with copy loss segments, smooth&merge segments within 300kb and prune
# segments smaller than 300kb.
segs.clean <- trim_to_coverage(segments, oncoscan_na33.cov) %>%
    adjust_loh() %>%
    merge_segments() %>%
    prune_by_size()

# Then we need to compute the arm-level alteration for loss and LOH since many 
# scores discard arms that are globally altered.
arms.loss <- names(get_loss_segments(segs.clean) %>%
        armlevel_alt(kit.coverage = oncoscan_na33.cov))
arms.loh <- names(get_loh_segments(segs.clean) %>%
        armlevel_alt(kit.coverage = oncoscan_na33.cov))

# Get the number of LST
lst <- score_lst(segs.clean, oncoscan_na33.cov)

# Get the number of HR-LOH
hrloh <- score_loh(segs.clean, arms.loh, arms.loss, oncoscan_na33.cov)

# Get the genomic LOH score
gloh <- score_gloh(segs.clean, arms.loh, arms.loss, oncoscan_na33.cov)

# Get the number of nLST
wgd <- score_estwgd(segs.clean, oncoscan_na33.cov)  # Get the avg CN, including 21p
nlst <- score_nlst(segs.clean, wgd["WGD"], oncoscan_na33.cov)

print(c(LST=lst, `HR-LOH`=hrloh, gLOH=gloh, nLST=nlst))
```


## TDplus score
function `score_td`

Procedure based on the paper from Popova et al., Cancer Res 2016 (PMID: 
26787835). The TDplus score is defined as the number of regions larger than 1Mb 
but smaller or equal to 10Mb with a gain of one or two copies. This score was 
linked to CDK12-deficient tumors. 
They also identified as second category of tandem duplication whose size is 
smaller or equal than 1Mb and around 300Kb but could not link it to a phenotype.
Note that due to its resolution the Oncoscan assay will most likely miss  this 
second category. Nonetheless it is reported by the function but not by the 
standard workflow.

```{r}
# Load data
chas.fn <- system.file("extdata", "TDplus_gene_list_full_location.txt", 
                       package = "oncoscanR")
segments <- load_chas(chas.fn, oncoscan_na33.cov)

# Clean the segments: restricted to Oncoscan coverage, LOH not overlapping
# with copy loss segments, smooth&merge segments within 300kb and prune
# segments smaller than 300kb.
segs.clean <- trim_to_coverage(segments, oncoscan_na33.cov) %>%
    adjust_loh() %>%
    merge_segments() %>%
    prune_by_size()

td <- score_td(segs.clean)
print(td$TDplus)
```


## Main workflow (as used at the Geneva University Hospitals)
The main workflow used for routine analysis can be launched either in R via the 
`workflow_oncoscan.chas(chas.fn, gender)` function or via the script 
`bin/run_oncoscan_workflow.R`:

Usage: 
`Rscript path_to_oncoscanR_package/bin/oncoscan-workflow.R CHAS_FILE`
- `CHAS_FILE`: Path to the text export file from ChAS or a compatible text file.

The script will output a JSON string into the terminal with all the computed 
information. :

```
{
  "armlevel": {
    "AMP": [],
    "LOSS": ["17p", "2q", "4p"],
    "LOH": ["14q", "5q", "8p", "8q"],
    "GAIN": [19p", "19q", "1q", "20p", "20q", "3q", "5p", "6p", "9p", "9q", 
    "Xp", "Xq"]
  },
  "scores": {
    "HRD": "Negative, nLST=12",
    "TDplus": 22,
    "avgCN": "2.43"
  },
  "file": "H19001012_gene_list_full_location.txt"
}
```

Or to launch the workflow within R:
```{r}
segs.filename <- system.file('extdata', 'LST_gene_list_full_location.txt', 
                             package = 'oncoscanR')
dat <- workflow_oncoscan.chas(segs.filename)

message(paste('Arms with copy loss:', 
              paste(dat$armlevel$LOSS, collapse = ', ')))
message(paste('Arms with copy gains:', 
              paste(dat$armlevel$GAIN, collapse = ', ')))
message(paste('HRD score:', dat$scores$HRD))
```

A similar function is available for running the workflow from an ASCAT result
file: `workflow_oncoscan.ascat`.

```{r}
library(jsonlite)
segs.filename <- system.file('extdata', 'ascat_example.txt', 
                             package = 'oncoscanR')
dat <- workflow_oncoscan.ascat(segs.filename)
toJSON(dat, auto_unbox=TRUE, pretty=TRUE)
```


Please read the manual for a description of all available R functions.

## References
1. "Homologous Recombination Deficiency (HRD) Score Predicts Response to 
Platinum-Containing Neoadjuvant Chemotherapy in Patients with Triple-Negative 
Breast Cancer.", M. Telli et al., Clin Cancer Res volume 22(15), august 2016.
2. "Ovarian Cancers Harboring Inactivating Mutations in CDK12 Display a Distinct
Genomic Instability Pattern Characterized by Large Tandem Duplications.", T. 
Popova et al., Cancer Res volume 76(7), april 2016.
3. "Ploidy and large-scale genomic instability consistently identify basal-like 
breast carcinomas with BRCA1/2 inactivation.", T. Popova et al., Cancer Res 
volume 72(21), november 2012.
4. "Patterns of genomic loss of heterozygosity predict homologous recombination 
repair defects in epithelial ovarian cancer.", V. Abkevich et al., Br J Cancer. 
2012 Nov 6;107(10).
5. "Absolute quantification of somatic DNA alterations in human cancer", S. 
Carter et al., Nat Biotech, 2012 volume 30(5).


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

```{r}
sessionInfo()
```