---
title: "Data Processing with SelectSim"
author: "Arvind Iyer"
email: "arvind.iyer@unil.ch"
date: "`r Sys.Date()`"
output: 
 rmarkdown::html_vignette:
    css: doc.css
vignette: >
  %\VignetteIndexEntry{Data Processing with SelectSim}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

- The goal of SelectSim package is to implement the SelectSim methodology to infer inter-dependencies between functional alterations in cancer.
- `SelectSim` package provides functions to generate the background model and other utility functions.

### Installation

- You can install the development version of SelectSim from [GitHub](https://github.com/CSOgroup/SelectSim) with:

``` r
devtools::install_github("CSOgroup/SelectSim",dependencies = TRUE, build_vignettes = TRUE)
```

- For more details on installation refer to [INSTALLATION](https://github.com/CSOgroup/SelectSim/blob/main/INSTALLATION.md)

### Example

- We will process LUAD MAF dataset from TCGA provided with the package.
- We will also use oncokb v3.9 cancer genes and mutations to filter the maf to create the gam provided with the package.

- NOTE: This is an example to process the MAF file. To know more about how we processed the MAF file to run_data object, please refer to SelectSim_analysis [repository](https://github.com/CSOgroup/SelectSim_analysis).


```{r}
library(SelectSim)
library(dplyr)
if (requireNamespace("tictoc", quietly = TRUE)) library(tictoc)
## Load the data provided with the package
data(luad_maf, package = "SelectSim")
data(oncokb_genes, package = "SelectSim")
data(oncokb_truncating_genes, package = "SelectSim")
data(variant_catalogue, package = "SelectSim")
```

```{r}
# Check the MAF
dim(luad_maf)
```

- Let print number of lines and number of samples

```{r}
input_maf <- luad_maf
print(paste('##### Number of lines ####',nrow(input_maf),sep="->"))
genes_to_consider =  oncokb_genes
print(paste('##### Number of genes ####',length(genes_to_consider),sep="->"))
```

- Let create a table schema of mutations to consider and columns defined in maf file

```{r}
mutation_type = list(
      'ignore' = c("Silent","Intron","RNA","3'UTR","5'UTR","5'Flank","3'Flank","IGR"),
      'truncating'= c('Frame_Shift_Del','Frame_Shift_Ins','In_Frame_Del','In_Frame_Ins','Nonsense_Mutation','Nonstop_Mutation','Splice_Region','Splice_Site','Translation_Start_Site'),
      'missense' = c('Missense_Mutation')
)
custom_maf_schema = list(
    'name' = 'custom_maf',
    'column' = list(
          'gene' = 'Hugo_Symbol'
        , 'gene.name' = 'Hugo_Symbol'
        , 'sample' = 'sample'
        , 'sample.name' = 'sample'
        , 'mutation.type' = 'Variant_Classification'
        , 'mutation' = 'HGVSp_Short'
        ),
        'mutation.type' = mutation_type
)
```

- Number of samples in the maf file

```{r}
mut_samples = unique(input_maf[, custom_maf_schema$column$sample])
print(paste('##### Number of samples ####',length(mut_samples),sep="->"))
```


- Filter the maf file to include oncokb cancer genes

```{r}
maf_genes = filter_maf_gene.name(input_maf, genes = genes_to_consider, gene.col = custom_maf_schema$column$gene)
print(paste('##### Number of lines ####',nrow(maf_genes),sep="->"))
```

#### Generating the GAMs

- Let generate the truncating data
  - We generate the GAM consider the genes to consider truncating mutations.
  - We also create a TMB dataframe which takes all truncating mutation of all genes in consideration.

```{r}
# Creating Truncating GAM
if (requireNamespace("tictoc", quietly = TRUE)) tictoc::tic('##### Creating Truncating GAM ####')
    maf_trunc = filter_maf_truncating(maf_genes,genes=oncokb_truncating_genes, custom_maf_schema)
    input_maf_trunc<-filter_maf_truncating(input_maf, custom_maf_schema)
    truncating_tmb <- data.frame('sample'=mut_samples,'mutation'=rep(0,length(mut_samples)))
    rownames(truncating_tmb)<-mut_samples
    temp <- input_maf_trunc %>% count(sample)
    rownames(temp)<-temp$sample
    truncating_tmb[intersect(truncating_tmb$sample,temp$sample),]$mutation <-temp[intersect(truncating_tmb$sample,temp$sample),'n']
    tcga_truc_gam = maf2gam(maf_trunc,
                     sample.col = custom_maf_schema$column$sample,
                     gene.col = custom_maf_schema$column$gene,
                     value.var = 'Variant_Classification',
                     samples = mut_samples,
                     genes = genes_to_consider,
                     fun.aggregate = length,
                     binarize=TRUE,
                     fill=0)
    truncating_data <- list('gam'=tcga_truc_gam,
                            'tmb'=truncating_tmb)
if (requireNamespace("tictoc", quietly = TRUE)) tictoc::toc()
```

- Let generate the Missense data
  - We generate the GAM with genes and hotspot mutations from oncokb v3.9.
  - We also create a TMB dataframe which takes all Missense mutation of all genes in consideration.

```{r}
# Creating Missense GAM
if (requireNamespace("tictoc", quietly = TRUE)) tictoc::tic('##### Creating Missense GAM ####')
    maf_valid = filter_maf_schema(input_maf,
                             schema = custom_maf_schema,
                             column = 'mutation.type',
                             values = custom_maf_schema[['mutation.type']][['ignore']],
                             inclusive = FALSE)
    missense_maf<-filter_maf_mutation.type(input_maf,
                                      variants = 'Missense_Mutation',
                                      variant.col = custom_maf_schema$column$mutation.type)
    missense_tmb <- data.frame('sample'=mut_samples,'mutation'=rep(0,length(mut_samples)))
    rownames(missense_tmb)<-mut_samples
    temp <- missense_maf %>% count(sample) 
    rownames(temp)<-temp$sample
    missense_tmb[intersect(missense_tmb$sample,temp$sample),]$mutation <-temp[intersect(missense_tmb$sample,temp$sample),'n']
    t_m = substr(maf_valid[[custom_maf_schema$column$mutation]],3,1000)
    t_m1 =  gsub('[A-Z]*$', '', t_m)
    maf_valid$HGVSp_Short_fixed = t_m1
    maf_hotspot = filter_maf_mutations(maf_valid,
                                  variant_catalogue,
                                  maf.col = c(custom_maf_schema$column$gene, 'HGVSp_Short_fixed'),
                                  values.col = c('gene', 'mut'))

    missense_tcga_gam = maf2gam(maf_hotspot,
                     sample.col = custom_maf_schema$column$sample,
                     gene.col = custom_maf_schema$column$gene,
                     value.var = 'Variant_Classification',
                     samples = mut_samples,
                     genes = genes_to_consider,
                     fun.aggregate = length,
                     binarize=TRUE,
                     fill=0)
    missesne_data <- list('gam'=missense_tcga_gam,
                          'tmb'=missense_tmb)

if (requireNamespace("tictoc", quietly = TRUE)) tictoc::toc()
```

#### Generating the run_object to run SelectSim

- We create a run_object data which is list object which consists of 
  - M: a list object of GAMs which is presence absence matrix of alterations
  - tmb: a list object of tumor mutation burden as data frame with column names (should be) as sample and mutation
  - sample.class a named vector of sample annotations
  - alteration.class a named vector of alteration annotations

```{r}
gene_to_take <- colnames(missesne_data$gam)
order <- rownames(missesne_data$gam)

data <-list('M'=list('missense'=t(missesne_data$gam[order,gene_to_take]),
                     'truncating'=t(truncating_data$gam[rownames(missesne_data$gam[order,]),gene_to_take])),
            'tmb'=list('missense'=missesne_data$tmb[order,],
                       'truncating'=truncating_data$tmb[order,]))

alteration_covariates <- rep('MUT',ncol(missesne_data$gam[order,gene_to_take]))
names(alteration_covariates)<-colnames(missesne_data$gam[order,gene_to_take])
sample_covariates<-rep('LUAD',length(order))
names(sample_covariates)<-order
run_data <- list('M'=data,'sample.class' = sample_covariates,'alteration.class' = alteration_covariates)
str(run_data)
```

- Save the `run_data` and check the introduction vignette to see how to run selectX to discover EDs. 

### Session information
```{r}
# Print the sessionInfo
sessionInfo()
```