Introduction: PCAWG signatures

In May 2018 a new and extended set of mutational signatures was published as a preprint by the Pan-Cancer Analysis of Whole Genomes (PCAWG) consortium (Alexandrov et al. 2020). In the following, we refer to that set of signatures as the PCAWG signatures. The underlying data is the largest available cohort of sequenced (Whole Genome Sequencing (WGS) and Whole Exome Sequencing (WES)) cancer samples to date and represents a substantial increase in statistical power as compared to (Alexandrov, Nik-Zainal, Wedge, Aparicio, et al. 2013) and (Alexandrov, Nik-Zainal, Wedge, Campbell, et al. 2013). On this data set, de novo mutational signature discovery analysis was performed with two different algorithms of the Non-negative Matrix Factorization (NMF) family of algorithms, SigPorfiler and SignatureAnalyzer (Bayesian NMF). In total 65 SNV mutational signatures (abbreciated SBS for single base substitution) were found, 49 of which were called validated, either because of good consensus between the two calling algorithms or because of presence in an orthogonal sequencing technology. In addition to the SNV or SBS mutational signatures, that analysis for the first time had sufficient statistical power to also extract 17 Indel (ID) signatures, based on a classification of Indels into 83 features, and Doublet Base Signatures (DBS) based on somatic di-nucleotide exchanges. In this vignette, we present an example analysis with the new PCAWG Indel signatures. The setup and line of arguments is analogous to the first vignette. However, YAPSA doesn’t provide functions for DBS signatures as mutation counts are often not sufficient for a deconvolution analysis.

As for the COSMIC SNV mutational signatures, the patterns for both the SNV and Indel PCAWG mutational signatures are stored in the software package and can be loaded into the workspace as follows:

data(sigs_pcawg)

For a precise description of all the data loaded by the above command, please refer to 2. Signature-specific cutoffs

Classification of Indels

A prerequisite for an analysis of Indel mutational signatures is a consistent classification system of Indels. This is non-trivial and critically depends on the choice of features, which have to include whether a given Indel is actually an insertion or a deletion, the size of the Indel,as well as the motif context in order to assess whether micro-homologies are present. In the classification proposed by the PCAWG consortium there are 83 features in total, which form five major groups:

1. Deletions of 1 bp C/(G) or T/(A) in a repetitive context. Taking into account the motif context of the deletion, the mutational event is further classified into categories of 1, 2, 3, 4, 5 or larger or equal to 6 identical nucleotide(s).
2. Insertions of 1 bp C/(G) or T/(A) in a repetitive context. Taking into account the motif context of the deletion, the mutational event is further classified into categories of 0, 1, 2, 3, 4, or larger or equal to 5 identical nucleotide(s).
3. Deletions of 2 bp, 3 bp, 4 bp or more or equal to 5 bp in a repetitive context. Each deletion is classified in a context of 1, 2, 3, 4, 5 or larger or equal to 6 times the same motif.
4. Insertions of 2 bp, 3 bp, 4 bp or more or equal to 5 bps in a repetitive context. Each deletion is classified in a context of 0, 1, 2, 3, 4 or larger or equal to 5 times the same motif.
5. Deletions of 2 bp, 3 bp, 4 bp or more or equal to 5 bp in a partially repetitive context with microhomology at the breakpoints. This partially repetitive context is defined by motif length of minus 1 bp, 2 bp, 3 bp, 4 bps or more or equal to 5 bp, located before and after the break-points of the deletion.

Indel classifications in the range of 1 - 5 bp are assigned to individual categories, while all variants larger than 5bp are considered as one category. Of note, microhomology is defined by partial sequence similarity between the motif of the respective Indel and the immediate neighbouring sequence context. In the following, we illustrate the microhomology classification with an example. We take a deletion of the motif ATGCGA, being more than 5 bp in length with microhomology, i.e., partially repetitive context of 4 bp at the breakpoint junction with the motifs ATGC upstream and GCGA downstream of the deletion. The annotated feature category is then MH_5+_4bp.

In YAPSA, the function plotExchangeSpectra_indel() can be used to plot the nucleotide exchange spectra of samples and/or signatures. The representation is analogous to the function plotExchangeSpectra() for SNV mutational signatures, while the colors and feature annotation is taken from (Alexandrov et al. 2020).

plotExchangeSpectra_indel(PCAWG_SP_ID_sigs_df[,c(3,6)])

: Nucleotide exchange spectra of the Indel signatures ID3, associated with tobacco smoking, and ID6, related to deficiencies in homologous recombination repair.

PCAWG Indel signatures

The publication (Alexandrov et al. 2020) reports signatures decomposed with an NMF approach (SigProfiler) and a Bayesian NMF approach (SignatureAnalyser). It is worth noting that YAPSA only has mutational signatures decomposed with SigProfiler integrated so far, of which the validated signatures are also found with SignatureAnalyser. Signatures only decomposed with SignatureAnalyzer are not part of YAPSA. Signature ID15 was excluded as none of the samples had a contribution to this signature. For seven of the Indel signatures, underlying mutational processes were asserted. Indel Signatures ID1 and ID2 can be found across all entities, while others such as ID13 is associated with UV light exposure and is found exclusively in melanoma (cf. Alexandrov et al. (2020)).

current_caption <- paste0("Information on Indel mutational signatures.")
if(!exists("repress_tables"))
  kable(PCAWG_SP_ID_sigInd_df, row.names=FALSE, caption=current_caption)

Example data: Genome of the Netherlands

In order to demonstrate the functionalities of the Indel signature analysis with YAPSA, a publicly available data set was chosen, the Genome of the Netherlands (https://www.nlgenome.nl/menu/main/app-go-nl/?page_id=9) release5. The data set is not related to cancer and instead contains germline variation in the genomes of the Dutch population. Hence no biological interpretation or conclusion of this analysis with respect to mutational processes should be made, it is of sole technical and demonstrative purpose.

Supervised Indel signature analysis

data(GenomeOfNl_raw)
GenomeOfNl_raw <- GenomeOfNl_raw[, c(1,2,4,5)]

load_data_new <- FALSE
if(load_data_new){
  data <- data.frame(matrix(ncol = 8, nrow = 0))
  for(index in seq_along(1:22)){
              print(index)
              temp <- tempfile()
              file_path <- paste0("https://molgenis26.gcc.rug.nl/
                                  downloads/gonl_public/variants/release5/
                                  gonl.chr",
                                  index, ".snps_indels.r5.vcf.gz")
              download.file(file_path, temp)
       
              data <- rbind(data, read.table(gzfile(temp, paste0("gonl.chr",
                                                    index,
                                                    ".snps_indels.r5.vcf")),
                                             header=FALSE, sep="\t", 
                                             stringsAsFactors = FALSE))
              data <- data[grep("INDEL", data$V8),]
              unlink(temp)
            
  }
  colnames(data) <- c("CHROM", "POS", "ID", "REF", "ALT", "QUAL", 
                      "FILTER","INFO")

  GenomeOfNl_raw <- data[, c(1,2,4,5)]
}

seed = 2
set.seed(seed)
number_of_indels <- sample(c(30:70), 15,  replace = TRUE)

index=0
seed=3
set.seed(seed)
vcf_like_indel_lists <- lapply(number_of_indels, function(size){
  df_per_PID <- GenomeOfNl_raw[sample(nrow(GenomeOfNl_raw), size, 
                                      replace = FALSE), ]
  index <<- index+1
  df_per_PID$PID <- rep(paste0("PID_", index), length(size))
  df_PIDs <- df_per_PID[order(df_per_PID$CHROM),]
  return(df_PIDs)
  })

vcf_like_indel_df <- do.call(rbind.data.frame, vcf_like_indel_lists)
kable(head(vcf_like_indel_df), caption="Head of the vcf_like_df
      containing the subsampled GoNL Indel data")

vcf_like_indel_trans_df <- translate_to_hg19(vcf_like_indel_df,"CHROM")
mutational_cataloge_indel_df <- create_indel_mutation_catalogue_from_df(
  in_dat = vcf_like_indel_trans_df,
  in_signature_df = PCAWG_SP_ID_sigs_df)

## [1] "Indel sequence context attribution of total  670  indels. This could take a while..."
## [1] "INDEL classification of total  670  INDELs This could take a while..."

kable(head(mutational_cataloge_indel_df[,1:5]))

data(cutoffs_pcawg)

current_catalogue_df <- mutational_cataloge_indel_df 
current_sig_df <- PCAWG_SP_ID_sigs_df
current_cutoff_pid_vector <- cutoffPCAWG_ID_WGS_Pid_df[3,]
current_sigInd_df <- PCAWG_SP_ID_sigInd_df

current_LCDlistsList <- LCD_complex_cutoff_combined(
  current_catalogue_df,
  current_sig_df,
  in_cutoff_vector = current_cutoff_pid_vector,
  in_filename = NULL,
  in_method = "abs",
  in_sig_ind_df = current_sigInd_df)

current_consensus_LCDlist <- current_LCDlistsList$consensus
if(!exists("repress_tables")) 
  as.character(current_consensus_LCDlist$out_sig_ind_df$sig)

##  [1] "ID1"  "ID2"  "ID3"  "ID4"  "ID5"  "ID6"  "ID8"  "ID9"  "ID12" "ID14"
## [11] "ID17"

Visualization

The exposures, i.e., the contributions of the different signatures to the mutational catalog of every sample, can be visualized with the function exposures_barplot(). All the functionalities which are available for SNV mutational signatures, such as the annotation of subgroups (cf.  1. Usage of YAPSA) are also available for the Indel mutational signatures.

exposures_barplot(current_LCDlistsList$perPID$exposures,
                  current_LCDlistsList$perPID$out_sig_ind_df)

:Exposures to Indel mutational signatures in the artificial data created by sampling GoNL variants. Exposures were obtained from a decomposition with PCAWG Indel signatures as well as their signature specific-cutoffs (cutoffPCAWG_ID_WGS_Pid_df).

Confidence Intervals

In order to assess trustworthiness of the computed exposures, YAPSA provides the calculation of confidence intervals (CIs). Analogously to CIs for SNV mutational signatures, the CIs for Indel mutational signatures are computed using the concept of profile likelihood.

confidence_intervals_ID <- confidence_indel_only_calulation(
  in_current_indel_df = current_catalogue_df)

## [1] "PCAWGValidID_abs"

plot(confidence_intervals_ID$p_complete_PCAWG_ID)

Confidence interval calculation for exposures to Indel mutational signatures

Note: With confidence_indel_only_calulation() only for INDELs the confidence intervals are calculated. To calculate the confidence intervals for SNVs and INDELs please use the function confidence_indel_calulation() and supply the mutational catalogs of both SNV and INDEL counts.