---
title: "pathwayTMB:Pathway Based Tumor Mutational Burden"
author: "Xiangmei Li,Junwei Han"
date: "`r Sys.Date()`"
output: 
  html_document:
    toc: true
    toc_depth: 3
    toc_float: true
    number_sections: true
    self_contained: yes
    css: corp-styles.css
    highlight: pygments
vignette: |
  %\VignetteIndexEntry{pathwayTMB} 
  %\VignetteEngine{knitr::rmarkdown}
   %\VignetteEncoding{UTF-8}
---

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

# Introduction
Somatic mutations play an important role in cancer development and disease progression, and mutational profiling is used far more commonly than other 'omics' analyses in clinical practice. With the development of PD-1 immunotherapy technology, many experiments show that patients with high TMB have been shown more clinical benefits of immunotherapy (PD-1/PD-L1). However, to date, WES remains confined to research settings, due to the high cost of the large genomic space sequenced. In the clinical setting, instead, targeted enrichment panels (gene panels) of various genomic sizes are emerging as the routine technology for TMB assessment. This package attempts to develop a new pathway-based gene panel for TMB assessment (pathway-based tumor mutational burden, **PTMB**), using somatic mutations files in an efficient manner from either TCGA sources or any in-house studies as long as the data is in MAF format. Besides, we develop a multiple machine learning method using sample PTMB profiles to identify cancer-specific dysfunction pathways, which can be a biomarker of prognostic and predictive for cancer immunotherapy.

# MAF field requirements
MAF files contain many fields ranging from chromosome names to cosmic annotations. However, most of the analysis in PTMB uses the following fields.

- Mandatory fields:   
  **Hugo_Symbol,**   
  **Variant_Classification,**  
  **Variant_Type,**  
  **Tumor_Sample_Barcode.**  
  
Complete specification of MAF files can be found on [NCI GDC documentation page](https://docs.gdc.cancer.gov/Data/File_Formats/MAF_Format/).

# Installation
```{r echo = TRUE, results = 'hide',eval=FALSE}
install.packages("pathwayTMB")
library(pathwayTMB)
```

# Overview of the package
The **pathwayTMB** package is a Bioinformatics tool to identify Cancer-specific pathways. And **pathwayTMB** functions can be categorized into mainly Visualization and Analysis modules. Each of these functions and a short description is summarized as shown below:   
<br/>
1.Obtain survival benefit mutations frequency matrix.  
2.Calculate the length of encoded genes.  
3.Infer pathway-based tumor mutational burden profiles (PTMB).  
4.Identify Cancer-specific dysfunction pathways.   
5.Visualization results:   
5.1 Plot Patients' Kaplan-Meier Survival Curves.  
5.2 Plot patient-specific dysfunction pathways and user-interested geneset mutually exclusive and co-occurrence plots.   
5.3 Plot patient-specific dysfunction pathways' waterfall plots.    
5.4 Plot the ROC curve.
<br/>

## Obtain survival-related genes' mutation frequency matrix   
We downloaded patients' mutation data from cBioPortal or other open source databases in MAF format. For each sample's mutation, we converted MAF format data into a mutation matrix, in which every row represents the gene and every column represents the sample. In our study, we only extract the nonsilent somatic mutations (nonsense mutation, missense mutation, frame-shift indels, splice site, nonstop mutation, translation start site, inframe indels) in protein-coding regions. Then a univariate Cox regression test is performed to identify the mutant genes with survival benefit using the sample's mutant and survival data. The function **get_mut_matrix** in the **pathwayTMB** package can implement the above process. Take simulated melanoma data as an example, the command lines are as follows:  
<br/>

```{r out.width=24,echo = TRUE}
#get the path of the mutation annotation file and samples' survival data
maf<-system.file("extdata","data_mutations_extended.txt",package = "pathwayTMB")
sur_path<-system.file("extdata","sur.csv",package = "pathwayTMB")
sur<-read.csv(sur_path,header=TRUE,row.names = 1)
#perform the function 'get_mut_matrix'
#mut_matrix<-get_mut_matrix(maffile=maf,is.TCGA=FALSE,mut_fre=0,nonsynonymous = TRUE,   cut_fisher.pval=1,cut_oddsRatio=1,sur=sur)
#view the first six lines of mutation matrix
head(mut_matrix)[1:6,1:6]
```
<br/>

## Calculate coding genes' length using gene transfer format (GTF) file.
The function **get_gene_length** can extract coding genes' length from gene transfer format (GTF) file. The GTF file can download from [GENCODE](https://www.gencodegenes.org/human/release_29.html).   
<br/>

```{r echo =TRUE, results = 'hide',eval=FALSE}
#calculate coding genes' length, filepath--the path of the GTF file
#get_gene_length(filepath)

```
<br/>

## Infer pathway-based tumor mutational burden profiles--PTMB.
To identify cancer-specific pathways, **PTMB** requires sample mutation data with survival information. For each pathway, **PTMB** is defined as pathway tumor mutational burden corrected by genes' length and number, the calculation formula is as follows:   
<br/>
$$PTMB_{ij}=\frac{\sum_{k=1}^N\frac{mut_{ijk}}{length_{jk}}}{N}\tag{1}$$   
where $mut_{ij}$ is the $i$ th sample's number of mutations of the $k$ th gene in the $j$ th pathway; $length_{jk}$ is the length of $k$ th gene in the $j$ th pathway; $N$ is the total number of genes in pathway $j$.The $PTMB_{ij}$ denotes the pathway $j$ tumor mutational burden for sample $i$, which was calculated in the context of survival-related genes' mutation frequency matrix.   
<br/>

```{r echo =TRUE, results = 'hold'}
#perform the function `get_PTMB`
#PTMB_matrix<-get_PTMB(freq_matrix=mut_matrix,genesmbol=genesmbol,gene_path=gene_path)
#show the first six lines of PTMB matrix
head(PTMB_matrix)[1:6,1:6]
```
<br/>

## Identify cancer-specific pathways
In order to identify cancer-specific dysfunction pathways, we develop a multiple machine learning method using pathway-based tumor mutational burden profiles. The multiple machine learning method consists of three sequential steps, as follows:   
<br/>

(1) Wilcoxon test is used to tentatively identify differentially PTMB pathways between death and alive sample groups.   
(2) Then random forest regression model is performed for further dysfunction pathways screening.   
(3) Lasso+cox regression model  is also performed for feature selection.  
<br/>

The function **get_final_signature**  in the **pathwayTMB** package can implement the above process.   
<br/>

```{r warning=FALSE}
# filter the survival-related pathways
#set.seed(1)
#final_character<-get_final_signature(PTMB=PTMB_matrix,sur=sur)
#view the final_character
final_character
```
<br/>

## Visualization results
(1) The function 'plotKMcurves' uses to draw Kaplan-Meier survival curves based on PTMB-related risk score. The risk score is generated by the signature's PTMB and the coefficient of 'Univariate' or 'Multivariate' cox regression. The command lines are as follows:   
<br/>

```{r fig.height=6, fig.width=8,warning=FALSE,results='hold'}
#Drawing Kaplan Meier Survival Curves using the final survival-related PTMB.
plotKMcurves(t(PTMB_matrix[final_character,]),sur=sur,returnAll = FALSE,risk.table = TRUE)
```

(2) The function 'plotMutInteract' uses Fisher Exact tests to detect mutually exclusive, co-occurring, and altered genesets or pathways, then visualize the results.   
<br/>

```{r fig.height=6, fig.width=8,warning=FALSE,results='hold'}
#a mutually exclusive co-occurrence chart showing the top 20 genes for mutation rates
gene_fre<-apply(mut_matrix,1,function(x){length(which(x!=0))/length(x)})
genes<-names(sort(gene_fre,decreasing = TRUE))[1:20]
plotMutInteract(freq_matrix=mut_matrix, genes=genes,returnAll = FALSE)
#pathways' mutually exclusive co-occurrence chart
plotMutInteract(freq_matrix=PTMB_matrix,genes=final_character,
                nShiftSymbols =0.3,returnAll = FALSE)
```

(3) The function 'GenePathwayOncoplots' takes output generated by read.maf and draws a GenePathwayOncoplots.   
<br/>

```{r fig.height=6, fig.width=8,warning=FALSE,results='hold'}
#calculate the PTMB-related riskscore
riskscore<-plotKMcurves(t(PTMB_matrix[final_character,]),sur=sur,plots = FALSE)$risk_score
cut_off<-median(riskscore)
#draw an GenePathwayOncoplots
GenePathwayOncoplots(maffile=maf,freq_matrix =mut_matrix,risk_score=riskscore,    
cut_off=cut_off,final_character=final_character,gene_path = gene_path,removeNonMutated = FALSE)
```

(4) The function 'plotROC' uses to plot a ROC curve.   
<br/>

```{r fig.height=6, fig.width=8, warning=FALSE, results='hold'}
#get the path of samples' immunotherapy response data
res_path<- system.file("extdata","response.csv",package = "pathwayTMB")
response<-read.csv(res_path,header=TRUE,stringsAsFactors =FALSE,row.name=1)
plotROC(riskscore=riskscore,response=response,main="Objective Response",print.auc=TRUE,grid = TRUE)
```
