dataAnalysis.Rmd

---
title: "## INSERT PROJECT TITLE ## - data analysis"
author: "Lorenzo Fanchi"
date: "`r format(Sys.Date(), '%d-%b-%Y')`"
output: 
  html_document: 
    theme: cerulean
    toc: yes
    toc_depth: 5
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
options(knitr.table.format = 'html')

setwd('~/projects/## PROJECT DIR ##')

############# KEEP IN MIND: correct data analysis depends on the addition of gs identifiers to germline variants #############
############# these should be added to the ID field (this has been done, if ngs-tools is used)                   #############

source('runConfig.R')
source('helperFunctions.R')

theme_set(theme_gray())
```

<link href="style.css" rel="stylesheet" type="text/css" />

<br>

#### Introduction

This document contains in-depth analysis of variants & predictions for the ## PROJECT TITLE ## project.

<br>

#### Sample information

```{r variant_snpEff_info}
snpeff_output = list.files(path = '4_snpEff',
                           pattern = 'genes[.]txt$',
                           full.names = T)
snpeff_output = snpeff_output[!grepl(pattern = 'snps', x = snpeff_output, fixed = TRUE)]

sample_info = fread('sample_info.tsv')
sample_info = sample_info[naturalorder(sample_info$patient_id)]

effects_table = lapply(snpeff_output,
                                 function(path) {
                                   data = fread(path, skip = 1)
                                   
                                   sample_prefix = gsub(pattern = regexPatterns$seqdata_prefix,
                                                        replacement = '\\1',
                                                        x = basename(path),
                                                        perl = T)
                                   
                                   if (sample_prefix %nin% sample_info$dna_data_prefix) {
                                     message(sample_prefix, ' not found in sample info file; please check DNA data prefix')
                                     return(data.table())
                                   }
                                   
                                   data.table(patient_id = sample_info[dna_data_prefix == sample_prefix,
                                                                       patient_id],
                                              syn_coding_variants = data %>%
                                                .[, grepl('synonymous_variant', names(.)), with = F] %>% 
                                                { if (length(.) < 1) { 0 } else { sum(.) } },
                                              non_syn_coding_variants = data %>%
                                                .[, grepl('missense_variant', names(.)), with = F] %>%
                                                { if (length(.) < 1) { 0 } else { sum(.) } },
                                              stop_gained_variants = data %>%
                                                .[, grepl('stop_gained', names(.)), with = F] %>%
                                                { if (length(.) < 1) { 0 } else { sum(.) } },
                                              stop_lost_variants = data %>%
                                                .[, grepl('stop_lost', names(.)), with = F] %>%
                                                { if (length(.) < 1) { 0 } else { sum(.) } },
                                              inframe_indel_variants = data %>%
                                                .[, grepl('inframe_(inser|dele)tion', names(.)), with = F] %>%
                                                { if (length(.) < 1) { 0 } else { sum(.) } },
                                              frameshift_variants = data %>%
                                                .[, grepl('frameshift_variant', names(.)), with = F] %>%
                                                { if (length(.) < 1) { 0 } else { sum(.) } })
                                 }) %>% rbindlist

effects_table = effects_table[naturalorder(effects_table$patient_id)]

# pander(sample_info, split.table = Inf, missing = '')

knitr::kable(effects_table) %>%
  kable_styling(bootstrap_options = c('striped', 'hover'))
```

```{r load_prediction_data, message = FALSE}
# parse data
project_predictions = parseEpitopePredictions(path = '3_neolution/predictions_output',
                                              pattern = '_epitopes[.]csv$')

# take subsets based on cutoffs specified in runConfig (same as used in predictions, but can be altered, if required)
project_predictions_subset = applyCutoffs(predictions = project_predictions,
                                          model = runOptions$neolution$model_cutoff,
                                          expression = runOptions$neolution$expression_cutoff,
                                          selfsim = runOptions$neolution$selfsim_filter)

# exclude alleles
project_predictions_subset = project_predictions_subset[!grepl(pattern = regexPatterns$allele_exclusion,
                                                               x = names(project_predictions_subset))]

# split by xmer
project_predictions_by_xmer = prepareEpitopeLists(list_of_predictions = project_predictions_subset)
```

```{r prepare_order_and_setup_lists}
# bind all together for order list
all_predictions = rbindlist(project_predictions_by_xmer)
all_predictions = all_predictions[naturalorder(patient_id)]
setkey(x = all_predictions, xmer, tumor_peptide_resin)

# take subsets for individual patients for combi-coding setup lists
predictions_by_patient = setNames(object = lapply(sample_info$dna_data_prefix,
                                                  function(prefix) {
                                                    data_subset = all_predictions[sample_prefix == prefix]
                                                    
                                                    data_unique = unique(x = data_subset,
                                                                         by = c('hla_allele', 'tumor_peptide'))
                                                    
                                                    setkey(x = data_unique, xmer, hla_allele)
                                                    return(data_unique)
                                                  }),
                                  nm = sample_info$patient_id)

# prepare unfiltered predictions to generate some stats
all_predictions_unfiltered = rbindlist(project_predictions)
setnames(x = all_predictions_unfiltered,
         old = grep('affinity|percentile_rank', names(all_predictions_unfiltered), value = T),
         new = c('tumor_affinity_nM', 'tumor_affinity_rank', 'normal_affinity_nM', 'normal_affinity_rank'))
all_predictions_unfiltered_by_patient = setNames(object = lapply(sample_info$patient_id,
                                                                 function(id) {
                                                                   data_subset = all_predictions_unfiltered[patient_id == id]
                                                                   
                                                                   return(data_subset)
                                                                 }),
                                                 nm = sample_info$patient_id)
```

<br>

#### Predicted epitopes summary
_Yield rates are defined as unique epitopes per non-synonymous coding variant_

```{r generate_prediction_summary, results = 'hide'}
prediction_table = merge(x = data.table(patient_id = names(predictions_by_patient),
                                        total_alleles = sapply(names(all_predictions_unfiltered_by_patient),
                                                         function(id) {
                                                           length(all_predictions_unfiltered_by_patient[[id]][, unique(hla_allele)])
                                                         }),
                                        `total_epitopes_\npredicted` = sapply(all_predictions_unfiltered_by_patient,
                                                                              function(dt) dt[!duplicated(tumor_peptide), .N]),
                                        selected_alleles = sapply(names(predictions_by_patient),
                                                         function(id) {
                                                           length(predictions_by_patient[[id]][, unique(hla_allele)])
                                                         }),
                                        `selected_epitopes_\npredicted` = sapply(predictions_by_patient,
                                                                                 function(dt) dt[!duplicated(tumor_peptide), .N])),
                         y = effects_table[, -c('syn_coding_variants'), with = F],
                         by = 'patient_id')
prediction_table = prediction_table[naturalorder(patient_id)]
prediction_table[, `total_epitope_\nyield_rate_per_allele` := round(x = `total_epitopes_\npredicted` / non_syn_coding_variants / total_alleles,
                                                                    digits = 1)]
prediction_table[, `selected_epitope_\nyield_rate_per_allele` := round(x = `selected_epitopes_\npredicted` / non_syn_coding_variants / selected_alleles,
                                                                       digits = 2)]

order_first = c('patient_id',
                'total_alleles', 'total_epitopes_\npredicted', 'total_epitope_\nyield_rate_per_allele',
                'selected_alleles', 'selected_epitopes_\npredicted', 'selected_epitope_\nyield_rate_per_allele')
setcolorder(x = prediction_table,
            neworder = c(order_first, 
                         setdiff(x = names(prediction_table),
                                 y = order_first)))
```

```{r knit_prediction_summary}
knitr::kable(prediction_table[, -grep('variants', names(prediction_table), value = T), with = F]) %>%
  kable_styling(bootstrap_options = c('striped', 'hover'))
```

<br>

__Total number of unique selected epitopes: `r nrow(unique(all_predictions,by = 'tumor_peptide'))`__  

<br>

##### Filter effect analysis

Following table shows percentages of predicted epitopes passing various filtering steps and medians for tumor_affinity_rank, tumor_processing_score & rna_expression  

```{r calculate_filter_consequences}
# annotate consequences of various filtering steps per variant per patient
all_predictions_unfiltered_by_patient = lapply(all_predictions_unfiltered_by_patient,
                                               function(dt) {
                                                 dt[model_prediction >= runOptions$neolution$model_cutoff,
                                                    passes_model_cutoff := TRUE]
                                                 dt[model_prediction < runOptions$neolution$model_cutoff,
                                                    passes_model_cutoff := FALSE]
                                                 
                                                 dt[rna_expression > runOptions$neolution$expression_cutoff
                                                    | is.na(rna_expression),
                                                    passes_expression_cutoff := TRUE]
                                                 dt[rna_expression <= runOptions$neolution$expression_cutoff,
                                                    passes_expression_cutoff := FALSE]
                                               })

filter_consequences = rbindlist(lapply(all_predictions_unfiltered_by_patient,
                                       function(dt) {
                                         data.table(patient_id = dt[, unique(patient_id)],
                                                    passes_model_cutoff = round(dt[passes_model_cutoff == TRUE
                                                                                   & !duplicated(tumor_peptide), .N] / dt[!duplicated(tumor_peptide), .N] * 100,
                                                                                digits =  2),
                                                    passes_expression_cutoff = round(dt[passes_expression_cutoff == TRUE
                                                                                        & !duplicated(tumor_peptide), .N] / dt[!duplicated(tumor_peptide), .N] * 100,
                                                                                     digits =  1),
                                                    passes_both_cutoffs = round(dt[passes_model_cutoff == TRUE
                                                                                   & passes_expression_cutoff == TRUE
                                                                                   & !duplicated(tumor_peptide), .N] / dt[!duplicated(tumor_peptide), .N] * 100,
                                                                                digits =  2),
                                                    median_tumor_affinity_rank = median(dt[!duplicated(tumor_peptide),
                                                                                           tumor_affinity_rank]),
                                                    median_tumor_processing_score = round(median(dt[!duplicated(tumor_peptide),
                                                                                                    tumor_processing_score]),
                                                                                          digits =  2),
                                                    median_rna_expression = round(median(dt[!duplicated(tumor_peptide), rna_expression],
                                                                                         na.rm = TRUE),
                                                                                  digits =  2))
                                       }))

knitr::kable(filter_consequences) %>%
  kable_styling(bootstrap_options = c('striped', 'hover'))
```

<br>

##### DNA VAF analysis

To determine whether differences in epitope yield rate (# of epitopes / coding snv / allele) are due to differences in tumor purity or sequencing depth, the VAFs of coding variants are compared between patients. Lower than average VAF could indicate lower sample purity and/or lower sequencing depth.

```{r dna_vaf_check}
# get all unique coding variants from varcontext output
varcontext_output = setNames(object = lapply(list.files(path = '2_varcontext', pattern = 'varcontext.tsv', full.names = TRUE),
                                             function(path) {
                                               data = fread(path, colClasses = c(chromosome = 'character'),
                                                            select = c('variant_id','chromosome','start_position', 'end_position',
                                                                       'ref_allele', 'alt_allele', 'dna_vaf', 'variant_classification'))
                                               data = unique(x = data, by = names(data))
                                               data[, patient_id := sample_info[dna_data_prefix == gsub(pattern = regexPatterns$seqdata_prefix,
                                                                                                        replacement = '\\1',
                                                                                                        x = basename(path),
                                                                                                        perl = T),
                                                                                patient_id]]
                                             }),
                             nm = list.files(path = '2_varcontext', pattern = 'varcontext.tsv'))

all_variants = rbindlist(varcontext_output)
all_variants$patient_id = factor(x = all_variants$patient_id, levels = unique(naturalsort(all_variants$patient_id)))

dna_vaf_summary = data.table(patient_id = sapply(varcontext_output, function(dt) dt[, unique(patient_id)]),
                             selected_epitope_yield_rate_per_allele = prediction_table[sapply(sapply(varcontext_output,
                                                                                                     function(dt) dt[, unique(patient_id)]),
                                                                                              function(id) grep(id, prediction_table$patient_id)),
                                                                                       `selected_epitope_\nyield_rate_per_allele`],
                             non_syn_coding_variants = prediction_table[sapply(sapply(varcontext_output,
                                                                                      function(dt) dt[, unique(patient_id)]),
                                                                               function(id) grep(id, prediction_table$patient_id)),
                                                                        non_syn_coding_variants],
                             median_dna_vaf = sapply(varcontext_output,
                                                     function(dt) round(x = median(dt[!grepl(pattern = regexPatterns$gs_identifier,
                                                                                             x =  variant_id),
                                                                                      dna_vaf],
                                                                                   na.rm = TRUE),
                                                                        digits = 2)))
dna_vaf_summary = dna_vaf_summary[naturalorder(patient_id)]
dna_vaf_summary$patient_id = factor(x = dna_vaf_summary$patient_id,
                                    levels = naturalsort(dna_vaf_summary$patient_id))

# pander(dna_vaf_summary, split.table = 200)
```

```{r plot_dna_vaf_scatter, warning = F, message = F, fig.align="center"}
ggplot(data = dna_vaf_summary, 
       mapping = aes(x = selected_epitope_yield_rate_per_allele,
                     y = median_dna_vaf,
                     color = patient_id)) +
  geom_point(aes(size = non_syn_coding_variants)) +
  scale_size(limits = c(1, 2000),
             range = c(1, 5)) +
  scale_x_continuous(limits = c(0, max(dna_vaf_summary$selected_epitope_yield_rate_per_allele))) +
  scale_y_continuous(limits = c(0, max(dna_vaf_summary$median_dna_vaf)))
```

```{r plot_dna_vaf_histograms, fig.width = 14, fig.height = 7, warning = F, message = F, fig.align="center"}
ggplot(data = all_variants[!grepl(regexPatterns$gs_identifier, variant_id)], 
       mapping = aes(x = dna_vaf, fill = patient_id)) +
  geom_histogram() +
  facet_wrap(~ patient_id, ncol = 3)
```

<br>

##### RNA expression median

To determine whether differences in selected epitope yield rate are due to differences in rna sequencing data quality, the median TPM values of a random draw of protein_coding genes are compared between patients.

```{r rna_median_check, cache = TRUE}
gtf_data = as.data.table(readGFF(runOptions$general$gtf_annotation))

rna_data = setNames(object = lapply(list.files(path = '1b_rnaseq_data/processed_salmon', pattern = '\\.tsv$', full.names = T), fread),
                    nm = sample_info[unlist(sapply(gsub(pattern = regexPatterns$seqdata_prefix,
                                                        replacement = '\\1',
                                                        x = basename(list.files('1b_rnaseq_data/processed_salmon', '\\.tsv$')),
                                                        perl = T),
                                                   function(prefix) grep(prefix, rna_data_prefix))),
                                     patient_id])

gene_sample = sample(gtf_data[gene_biotype == 'protein_coding', gene_id], size = 2500)

protein_coding_medians = data.table(patient_id = names(rna_data),
                                    median_rna_expression = sapply(rna_data,
                                                                   function(dt) round(median(x = dt[gene_id %in% 
                                                                                                      gene_sample, tpm]),
                                                                                      digits = 1)))

protein_coding_medians = merge(x = prediction_table[, .(patient_id, non_syn_coding_variants,
                                                        selected_epitope_yield_rate_per_allele = `selected_epitope_\nyield_rate_per_allele`)],
                               y = protein_coding_medians,
                               by = 'patient_id')
protein_coding_medians = protein_coding_medians[naturalorder(patient_id)]
protein_coding_medians$patient_id = factor(x = protein_coding_medians$patient_id,
                                           levels = naturalsort(protein_coding_medians$patient_id))

# pander(protein_coding_medians, split.table = 200)
```

```{r plot_rna_median_scatter, warning = F, message = F, fig.align="center"}
ggplot(data = protein_coding_medians, 
       mapping = aes(x = selected_epitope_yield_rate_per_allele,
                     y = median_rna_expression,
                     color = patient_id)) +
  geom_point(aes(size = non_syn_coding_variants)) +
  scale_size(limits = c(1, 2000),
             range = c(1, 5)) +
  scale_x_continuous(limits = c(0, max(protein_coding_medians$selected_epitope_yield_rate_per_allele))) + 
  scale_y_continuous(limits = c(0, max(protein_coding_medians$median_rna_expression)))
```

<br>

##### Driver gene analysis

List of driver genes was obtained from [COSMIC: Cancer Gene census](http://cancer.sanger.ac.uk/census) and cross-checked with hugo_symbols observed in varcontext output

```{r check_oncogenes_selected_patients}
# load oncogene db list, obtained from http://cancer.sanger.ac.uk/census
oncogene_census = list.files(path = '.', pattern = 'oncogene-census\\.tsv$', full.names = TRUE)

if (length(oncogene_census) < 1) {
  stop('Please download recent oncogene census file from COSMIC and save as "oncogene_census.tsv" in working directory')
} 

oncogene_db = fread(oncogene_census)
oncogene_db = oncogene_db[, list(Synonyms = unlist(strsplit(Synonyms, ','))), by = `Gene Symbol`]

# load varcontext data
varcontext_full_data = setNames(object = lapply(list.files(path = '2_varcontext', pattern = 'varcontext\\.tsv', full.names = TRUE),
                                                fread, colClasses = c(chromosome = 'character')),
                                nm = sapply(gsub(pattern = regexPatterns$seqdata_prefix,
                                                 replacement =  '\\1',
                                                 x =  list.files(path = '2_varcontext', pattern = 'varcontext\\.tsv', full.names = FALSE),
                                                 perl = T),
                                            function(prefix) sample_info[dna_data_prefix == prefix, patient_id]))

varcontext_full_data = varcontext_full_data[naturalorder(names(varcontext_full_data))]

driver_table = data.table(patient_id = names(varcontext_full_data),
                          non_syn_coding_variants = sapply(names(varcontext_full_data),
                                                           function(id) effects_table[patient_id == id, non_syn_coding_variants]),
                          driver_genes = sapply(varcontext_full_data,
                                                function(dt) {
                                                  dt_subset = dt[variant_classification != 'silent_mutation'
                                                                 & transcript_remark != 'identical'
                                                                 & !duplicated(hugo_symbol)]
                                                  drivers_symbol = sort(dt_subset$hugo_symbol)[which(sort(dt_subset$hugo_symbol) %in% oncogene_db[, `Gene Symbol`])]
                                                  drivers_synonym = sort(dt_subset$hugo_symbol)[which(sort(dt_subset$hugo_symbol) %in% oncogene_db[, Synonyms])]
                                                  
                                                  return(paste0(unique(c(drivers_symbol, drivers_synonym)), collapse = ', '))
                                                }))

knitr::kable(driver_table) %>%
  kable_styling(bootstrap_options = c('striped', 'hover'))
```

<br>

##### Mutational signature analysis

Following analysis shows distribution of mutational signatures of all SNVs (excluding germline SNPs) for all patients. Signatures from Alexandov et al. Nature 2013 have been used in the analysis.  

The top panel is the tumor mutational profile displaying the fraction of mutations found in each trinucleotide context, the middle panel is the reconstructed mutational profile created by multiplying the calculated weights by the signatures, and the bottom panel is the error between the tumor mutational profile and reconstructed mutational profile, with SSE (sum-squared error) shown.  

For mutational signature information, click <a href="http://cancer.sanger.ac.uk/cosmic/signatures" target="_blank">here</a>  

```{r mutational_signature_analysis, warning = FALSE, message = FALSE, cache = TRUE, fig.width = 10, fig.heigth = 7}
mutational_signatures = mutationalSignatureAnalysis(table = rbindlist(varcontext_output)[!grepl(regexPatterns$gs_identifier, variant_id)
                                                                                         & (nchar(ref_allele) == 1 & nchar(alt_allele) == 1)],
                                                    genome_build = )
mutational_signatures = mutational_signatures[naturalorder(names(mutational_signatures))]

analysed_mutation_counts = sapply(names(mutational_signatures), 
                                  function(id) nrow(rbindlist(varcontext_output)[patient_id == id
                                                                                 & !grepl(regexPatterns$gs_identifier, variant_id)
                                                                                 & (nchar(ref_allele) == 1 & nchar(alt_allele) == 1)]))

for (i in 1:length(mutational_signatures)) {
  plotSignaturesSmallLabels(mutational_signatures[[i]], sub = paste('Analysed mutations:', analysed_mutation_counts[names(mutational_signatures)[i]]))
}
```

<br>

##### Immune response pathway mutation analysis


```{r pathway_mutation_analysis, result = 'asis'}
ifng_symbols = c('IFN', 'JAK', 'STAT', 'IRF', 'SOCS', 'DAPK', 'SLFN')
ag_pres_symbols = c('B2M', 'HLA', 'PSM', 'TAP', 'ERAP')

vcfs_parsed = parseVcfs(vcf_path = '4_snpEff', extract_info = T, extract_fields = F, write = F)

vcfs_muts = lapply(vcfs_parsed,
                   function(vcf_dt) {
                     vcf_ifn = vcf_dt[grep(pattern = paste0(ifng_symbols, collapse = '|'), x = vcf_dt$info)]
                     vcf_ag_pres = vcf_dt[grep(pattern = paste0(ag_pres_symbols, collapse = '|'), x = vcf_dt$info)]
                     
                     vcf_muts = rbindlist(list(vcf_ifn, vcf_ag_pres))
                     vcf_muts = vcf_muts[grep(pattern = 'MODERATE|HIGH', info), .(chromosome, start_position, info)]
                     
                     return(vcf_muts)
                   })

for (i in 1:length(vcfs_muts)) {
  knitr::kable(vcfs_muts[[i]]) %>%
    kable_styling(bootstrap_options = c('striped', 'hover')) %>%
    print()
}
```