Skip to content

tools.md

Latest commit

 

History

History
232 lines (138 loc) · 12.4 KB

tools.md

File metadata and controls

232 lines (138 loc) · 12.4 KB
Raw

bbglab/deepCSA: Tools explanation

Table of contents

Introduction

Here, you can find an explanation of the different computations, tools or metrics implemented in deepCSA.

Interpreting outputs (sanity checks and key metrics)

Sanity checks / QC

Use these outputs to assess overall data quality before interpreting biological signals:

  • Depth summaries (depthssummary/): verify consistent coverage across samples and genes.
  • Mutation density vs depth (qc/metrics_vs_depth/): check that mutation density does not collapse in low-depth samples.
  • Omega QC (qc/metrics_vs_depth/ + qc/annotated_omegas): highlights genes/samples with unstable omega estimates.
  • Mutational profile stability (computeprofile/*.profile_stability.tsv): higher deviations indicate unstable mutational profiles (see below).

Omega vs omegagloballoc

  • omega/ uses per-sample mutational profiles and per-sample synonymous rates to estimate selection.
  • omegagloballoc/ uses a global cohort mutational profile and global synonymous rates (shared across samples), which stabilizes estimates in low-burden samples and facilitates cohort-level comparisons.

Use omega for sample-specific selection signals and omegagloballoc for conservative cohort-level estimates.

Site selection values

Outputs in sitecomparison/ and sitecomparisongloballoc/ compare observed vs expected mutations per site or residue:

  • OBSERVED_MUTS: number of observed mutations.
  • EXPECTED_MUTS: expected mutations from mutability models.
  • OBS/EXP: selection enrichment ratio.
  • p_value: Poisson p-value for observing at least OBSERVED_MUTS given EXPECTED_MUTS.

The resolution is controlled by site_comparison_grouping (site, aminoacid, or aminoacid_change).

Mutational signatures

  • sigprofilerassignment/: assignments of known COSMIC signatures; includes activity tables and plots.
  • signatures_hdp/: extracted signatures using a hierarchical Dirichlet process.
  • sigprobs/ / muts2sigs/: per-mutation signature probabilities (useful for downstream stratification).

Interpret signature results alongside mutation counts and profile stability to avoid over-interpreting low-burden samples.

Mutational profile stability

The file *.profile_stability.tsv is generated by adding a single mutation to each of the 96 SBS channels and measuring the L1 deviation from the original profile. Reported statistics include:

  • mean_deviation, min_deviation, max_deviation, std_deviation

Lower deviations indicate a more stable (less noisy) profile.

Publications with detailed explanation

We are in the process of completing the documentation, but in the meantime you can check the recently published paper and its supplementary material for more details.

Adjusted mutation density

This can be found in the mutdensityadj folder of the deepCSA output.

Goal

The goal is to define a mean mutation density estimate for a specific subset of mutation sites. For example, we might be interested in measuring what is the mutation density at the set of missense mutation sites -- i.e., mutation sites in the CDS that induce missense mutations.

Motivation

In the context of bulk-ultradeep sequencing, there are several factors that can influence the occurrence of mutations at a given site.

  • trinucleotide context of the site

  • neutral mutagenesis, defined as:

    • normalized profile: vector of relative mutabilities for each trinucleotide context

    • exposure associated with the profile, expressed as mutation burden

  • sequencing depth of the site

  • selection of the site

We would like to define a way to compute the mutation density that corrects for potential confounders like depth and triplet content, thus rendering the estimates more comparable. For example, the missense sites in gene A may imply that more missense mutations per site are observed than in gene B, even if both genes are subject to the same exact mutational process.

Given a collection $\mathcal{C}=\{i\}$ of mutation sites (e.g. missense mutations) and a number of mutations $m$ observed in this context, we would like to come up with an expression of the type $m/L$, where $L$ denotes the effective number of sites.

Assuming even depth $D$ for simplicity, the current flat method would simply take :

$$L=\frac{1}{3}\cdot|\mathcal{C}|\cdot D$$

Assumptions

Throughout we make the following assumptions:

  • the only way mutations come about is by neutral mutagenesis

  • discrepancies between neutral mutagenesis and the observed mutations are due to selection

Definitions/Notation

Mutation sites are defined as specific single base nucleotide changes. Each mutation site corresponds to a tri-nucleotide context of the form $F_5RF_3>A$ where $F_5, F_3$ are the 5' and 3' flanking nucleotides, respectively, R is the reference nucleotide and A is the alternate allele. Therefore, there are 3 mutation sites per genomic position.

Let $\{P_c\}$ be the relative mutabilities per trinucleotide context $c$: given a mutation site with tri-nucleotide context $c_1$, we expect to observe $P_{c_1}/P_{c_2}$ more mutations than in a mutation site with triplet context $c_2$. For a more univocal definition, we may want to impose the following condition:

$$\sum_c P_c = 1$$

Note that in order to render absolute mutabilities, the $P_c$ should be multiplied by a suitable scaling factor $\alpha$.

Let $D_i$ be the sequencing depth at mutation site $i$.

Let $m$ denote the total number of mutations observed in the target class -- i.e., mapping to the prescribed collection of mutation sites where we are interested to compute the mutation density.

For each mutation site $i$, we will denote $c(i)$ its triplet context.

Method description

Given a collection $\mathcal{C}=\{i\}$ of mutation sites (e.g. missense mutations) and a number of mutations $m$ observed in this context, we want to compute an effective length $L$ that represents the entire collection of sites

$$L=\sum_{i\in\mathcal{C}} P_{c(i)}\cdot D_i$$

Note that this definition of length has an arbitrary interpretation, because of the multiple possible choices that we can make for the relative mutational profile. To render a universal length, we must set a scale to represent the mutational profile with a concrete meaning. What would be a good practical reference?

One possibility would be to choose a scaling factor $\alpha$ that represents an absolute mutability equivalent to an expected 1 mutation per megabase in one mappable copy of the genome.

If $n_c$ is the number of mutation sites with trinucleotide context $c$ comprised in the mappable genome and $G$ is the total length of the mappable genome (expressed in nucleotides), $\alpha$ must be such that:

$$\alpha\sum_c n_cP_c = G / 10^6$$

Because $G=\frac{1}{3}\sum_c n_c$, we can compute $\alpha$:

$$\hat\alpha = \frac{G}{10^6\cdot\sum_c n_cP_c} = \frac{1}{3\cdot 10^6} \frac{\sum_c n_c}{\sum_c n_cP_c}$$

Applying this scaling, the effective length would take the following form:

$$L=\hat\alpha\sum_i P_{c(i)}\cdot D_i = \hat\alpha\sum_c \left( P_c\cdot \sum_{i\in c} D_i \right)$$

Units

Defining the mutation density as $m/L$, according to the preceding explanation, the units in which the readout is expressed are the following: mutations / Mb sequenced.

Omega

For more explanations on omega go to the corresponding repo.

deepCSA applies Benjamini-Hochberg multiple-testing correction separately for each of the following comparison sets: all-samples (cohort), sample groups, and per-sample results, and it does this independently for gene-level and subgenic regions. The corrected values are reported in the pvalue_adj column. P-values equal to 0 are set to 1.17e-38 (minimum non-zero float32) before correction to avoid underflow.

Site comparison

The site comparison step takes advantage of the computation of mutabilities in omega, and then compares these mutabilities either by residue, residue change or nucleotide change.

dNdScv

deepCSA wraps the dNdScv R package and runs it with a dynamically built RefCDS reference instead of relying on a pre-baked .rda transcripts file.

The dnds subworkflow performs three steps for every run:

  1. ADAPT_PANEL_REFCDS (dNdScv_panel_prep.py) — filter the biomart export referenced by params.dnds_biomart_ref to the transcripts overlapping the panel BED.
  2. BUILD_REFCDS — call dndscv::buildref using params.fasta to produce a fresh RefCDS_custom.rda.
  3. DNDSRUN (dNdS_run.R) — run dNdScv on the cohort, producing *.cv.tsv, *.globaldnds.tsv and *.loc.tsv under selection/dndscv/{cv,persample,local}/.

Instructions for regenerating the biomart TSV are in assets/build_datasets/dndscv/instructions.txt. The previously required dnds_ref_transcripts parameter has been removed.

dN/dS proxy

When mutation density and the all-regions profile are computed, deepCSA also generates a quick dN/dS proxy per gene by taking the ratio of non-synonymous vs synonymous adjusted mutation densities. The implementation is in mut_density_adjusted_dnds.py and the results are published to selection/dndsproxy/ as *.gene_mutdensities_n_dnds.tsv.

This metric is intended as a fast sanity check and is independent of the R-based dNdScv run and of omega, both of which provide dN/dS estimates with significance testing. It is gated by run_mutdensity (which itself is enabled by either mutationdensity or omega) combined with profileall.

Depth-vs-metric QC

The qc/metrics_vs_depth/ directory is produced by PLOT_METRICS_VS_DEPTH_QC (in the plotting_qc subworkflow) and is always generated. It joins per-gene/sample average depth (from the PLOTDEPTHSEXONSCONS step) with:

  • raw mutation densities (mutdensity/)
  • adjusted mutation densities (mutdensity_adjusted/)
  • omega-globalloc estimates (selection/omegagloballoc/)

Each combination yields a scatter PDF and a status TSV under *.metrics_depth_qc/, used to flag samples/genes whose metric values may be confounded by sequencing depth.

Mutational signatures

We provide two different strategies for signature analysis.

  • Using SigProfilerAssignment with a set of known SBS signatures

  • Using a Hierarchical Dirichlet Process algorithm developed by Nicola Robets and compacted by the McGranahan lab into a wrapped version.

    • The outputs of the signature extraction process are then further processed downstream using SigProfilerAssignment to decompose the de novo signature and reassign mutational processes to samples.

  • Additionally we also output mutation count matrices that are ready to be run through MSA which is another method for mutational signature attribution.

Additionally one could run SigProfilerExtractor on the data but this needs to be done externally.

Containers and reproducibility

deepCSA defines container images directly in module files and conf/modules.config. For bbglab-maintained images (bbglab/*), Dockerfile recipes are tracked in the lab repository: https://github.com/bbglab/containers-recipes. External images (e.g., ferriolcalvet/*, rblancomi/*, biocontainers/*) should be mirrored locally if strict reproducibility is required.

Key images used by the pipeline:

Component Image
Core utilities docker.io/bbglab/deepcsa-core:0.1.0
Panel BED tools docker.io/bbglab/deepcsa_bed:latest
Omega docker.io/bbglab/omega:0.2.1
Oncodrive3D docker.io/bbglab/oncodrive3d:1.0.5
Oncodrive3D (ChimeraX plots) docker.io/spellegrini87/oncodrive3d_chimerax:latest
OncodriveFML docker.io/ferriolcalvet/oncodrivefml:latest
OncodriveCLUSTL docker.io/ferriolcalvet/oncodriveclustl:latest
SigProfilerAssignment docker.io/ferriolcalvet/sigprofiler_assignment:1.1.3
SigProfilerMatrixGenerator docker.io/ferriolcalvet/sigprofilermatrixgenerator:1.3.5
mSigHdp (HDP) docker.io/ferriolcalvet/msighdp:latest
bbgregressions docker.io/rblancomi/bbgregressions:dev
Ensembl VEP biocontainers/ensembl-vep:111.0--pl5321h2a3209d_0 (version depends on vep_cache_version)
SAMtools biocontainers/samtools:1.18--h50ea8bc_1

To override any image, set process.container or the relevant module label in your nextflow.config.