Push updates to main

.gitignore

@@ -15,3 +15,4 @@ scratch/
 scratchhhh/
 tests/test_data/all_samples.somatic.mutations.maf
 tests/test_data/all_samples_indv.depths.tsv.gz
+assets/useful_scripts/*.ipynb

CITATIONS.md

@@ -10,6 +10,10 @@
 
 ## Sources of data and tools
 
+- **Ensembl VEP**
+
+  > McLaren W, Gil L, Hunt SE, et al. The Ensembl Variant Effect Predictor. Genome Biol. 2016;17:122. doi: 10.1186/s13059-016-0974-4.
+
 - Nanoseq masks
 
   > Abascal, F., Harvey, L.M.R., Mitchell, E. et al. Somatic mutation landscapes at single-molecule resolution. Nature 593, 405–410 (2021). https://doi.org/10.1038/s41586-021-03477-4
@@ -22,7 +26,7 @@
 
   > https://cancer.sanger.ac.uk/signatures/sbs
 
-- **dNdScv covariates**
+- **dNdScv method + covariates**
 
   > Martincorena I, et al. (2017) Universal Patterns of Selection in Cancer and Somatic Tissues. Cell. http://www.cell.com/cell/fulltext/S0092-8674(17)31136-4
 
@@ -38,11 +42,55 @@
 
   > Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016 Oct 1;32(19):3047-8. doi: 10.1093/bioinformatics/btw354. Epub 2016 Jun 16. PubMed PMID: 27312411; PubMed Central PMCID: PMC5039924.
 
-- Python
-- SigProfilerAssignment, MatrixGenerator
-- HDP
-- OncodriveFML
-- OncodriveCLUSTL
+- **bgreference / bgdata**
+
+  > Repository: https://github.com/bbglab/bgreference  
+  > Repository: https://github.com/bbglab/bgdata
+
+- **SAMtools**
+
+  > Li H, Handsaker B, Wysoker A, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25(16):2078-9. doi: 10.1093/bioinformatics/btp352.
+
+- **BEDTools**
+
+  > Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841-842. doi: 10.1093/bioinformatics/btq033.
+
+- **HTSlib / Tabix**
+
+  > Li H. Tabix: fast retrieval of sequence features from generic TAB-delimited files. Bioinformatics. 2011;27(5):718-719. doi: 10.1093/bioinformatics/btq671.
+
+- **Python**
+
+  > Python Software Foundation. Python Language Reference, version 3.x. https://www.python.org/
+
+- **SigProfilerAssignment, SigProfilerMatrixGenerator**
+
+  > Alexandrov LB, et al. The repertoire of mutational signatures in human cancer. Nature 578, 94–101 (2020). doi:10.1038/s41586-020-1943-3.
+  > SigProfiler tools (SigProfilerMatrixGenerator, SigProfilerAssignment). Alexandrov Lab. https://github.com/AlexandrovLab
+  > Repository: https://github.com/AlexandrovLab/SigProfilerAssignment  
+  > Repository: https://github.com/AlexandrovLab/SigProfilerMatrixGenerator
+
+
+- **HDP (Hierarchical Dirichlet Processes)**
+
+  > https://github.com/nicolaroberts/hdp
+  >   Roberts, N. D. (2018). Patterns of somatic genome rearrangement in human cancer. https://doi.org/10.17863/CAM.22674
+
+- **OncodriveFML**
+
+  > Repository: https://github.com/bbglab/oncodrivefml (see repository for citation details)
+  > Mularoni L, Sabarinathan R, Deu-Pons J, Gonzalez-Perez A, Lopez-Bigas N. OncodriveFML: a general framework to identify coding and non-coding regions with cancer driver mutations. Genome Biology. 2016;17:128. doi:10.1186/s13059-016-0994-0. https://github.com/bbglab/oncodrivefml
+
+- **OncodriveCLUSTL**
+
+  > Repository: https://github.com/bbglab/oncodriveclustl (see repository for citation details)
+  > Claudia Arnedo-Pac, Loris Mularoni, Ferran Muiños, Abel Gonzalez-Perez, Nuria Lopez-Bigas, OncodriveCLUSTL: a sequence-based clustering method to identify cancer drivers, Bioinformatics, Volume 35, Issue 22, November 2019, Pages 4788–4790, https://doi.org/10.1093/bioinformatics/btz501; https://github.com/bbglab/oncodriveclustl
+
+- **Omega (dN/dS)**
+
+  > Repository: https://github.com/bbglab/omega (see repository for citation details)
+
+
 
 ## Software packaging/containerisation tools
 

README.md

@@ -6,6 +6,34 @@
 
 ![deepCSA workflow overview](docs/images/deepCSA.png)
 
+## Documentation
+
+Find the documentation ([link to docs](https://github.com/bbglab/deepCSA/tree/main/docs)).
+
+We are working to provide the biggest possible detail on the [usage](docs/usage.md) and explanation of the rationale and [tools](docs/tools.md).
+
+You can also find an explanation with examples on the [output](docs/output.md).
+
+For more examples and description of the entire process, you can check the publications listed below.
+
+## Publications
+
+> **Sex and smoking bias in the selection of somatic mutations in human bladder**
+>
+> Ferriol Calvet*, Raquel Blanco Martinez-Illescas*, Ferran Muiños, Maria Tretiakova, Elena S. Latorre-Esteves, Jeanne Fredrickson, Maria Andrianova, Stefano Pellegrini, Axel Rosendahl Huber, Joan Enric Ramis-Zaldivar, Shuyi Charlotte An, Elana Thieme, Brendan F. Kohrn, Miguel L. Grau, Abel Gonzalez-Perez, Nuria Lopez-Bigas & Rosa Ana Risques
+>
+>_Nature_ (2025) doi:[10.1038/s41586-025-09521-x](https://doi.org/10.1038/s41586-025-09521-x)
+>
+> *these authors contributed equally and the order was decided randomly
+
+&
+
+> **DeepClone, an end-to-end protocol to study somatic mutagenesis and selection at high resolution**
+>
+> Ferriol Calvet, Morena Pinheiro-Santin, Erika Lopez, Raquel Blanco Martinez-Illescas, Núria Samper, Miguel L. Grau, Ferran Muiños, Rocío Chamorro González, Maria Andrianova, Federica Brando, Stefano Pellegrini, Marta Huertas, Elisabet Figuerola-Bou, Coohleen Coombes, Brendan F. Kohrn, Jeanne Fredrickson, Rosa Ana Risques, Nuria Lopez-Bigas, Abel Gonzalez-Perez
+>
+> _protocols.io_ (2026) doi:[10.17504/protocols.io.dm6gp1jodgzp/v2](https://dx.doi.org/10.17504/protocols.io.dm6gp1jodgzp/v2)
+
 ## Usage
 
 You can find a detailed documentation in the [docs section](docs/README.md), but here there is a minimal summary on how to prepare the inputs. Still for your first runs if you need to make the complete set up you have to check the deeper documentation.
@@ -22,6 +50,8 @@ sample2,sample2.filtered.vcf,sample2.sorted.bam
 
 Each row represents a single sample with a single-sample VCF containing the mutations called in that sample and the BAM file that was used for getting those variant calls. The mutations will be obtained from the VCF and the BAM file will be used for computing the sequencing depth at each position and using this for the downstream analysis.
 
+Two alternative input modes are also supported: a samplesheet with only `sample,vcf` columns combined with a precomputed depths table, or a single cohort-level MAF file passed via `--input_maf` together with a precomputed depths table. See [Input scenarios](docs/input_scenarios.md) for details.
+
 **Make sure that you do not use any '.' in your sample names, and also use text-like names for the samples, try to avoid having only numbers.** This second case should be handled properly but using string-like names will ensure consistency.
 
 **There are specific datasets that need to be prepared before running deepCSA. You can find a list of those, and instructions for downloading them in [the documentation section of the repo](docs/usage.md#mandatory-parameter-configuration).**
@@ -57,31 +87,3 @@ We thank the following people for their extensive assistance in the development
 An extensive list of references for the tools used by the pipeline can be found in the [`CITATIONS.md`](CITATIONS.md) file.
 
 This pipeline uses code and infrastructure developed and maintained by the [nf-core](https://nf-co.re) community, reused here under the [MIT license](https://github.com/nf-core/tools/blob/master/LICENSE).
-
-> **The nf-core framework for community-curated bioinformatics pipelines.**
->
-> Philip Ewels, Alexander Peltzer, Sven Fillinger, Harshil Patel, Johannes Alneberg, Andreas Wilm, Maxime Ulysse Garcia, Paolo Di Tommaso & Sven Nahnsen.
->
-> _Nat Biotechnol._ 2020 Feb 13. doi: [10.1038/s41587-020-0439-x](https://dx.doi.org/10.1038/s41587-020-0439-x).
-
-## Documentation
-
-Find the documentation ([link to docs](https://github.com/bbglab/deepCSA/tree/main/docs)).
-
-We are working to provide the biggest possible detail on the [usage](docs/usage.md) and explanation of the rationale and [tools](docs/tools.md), but this is still in progress.
-
-## Publications
-
-> **Sex and smoking bias in the selection of somatic mutations in human bladder**
->
-> Ferriol Calvet*, Raquel Blanco Martinez-Illescas*, Ferran Muiños, Maria Tretiakova, Elena S. Latorre-Esteves, Jeanne Fredrickson, Maria Andrianova, Stefano Pellegrini, Axel Rosendahl Huber, Joan Enric Ramis-Zaldivar, Shuyi Charlotte An, Elana Thieme, Brendan F. Kohrn, Miguel L. Grau, Abel Gonzalez-Perez, Nuria Lopez-Bigas & Rosa Ana Risques
->
->_Nature_ (2025) doi:[10.1038/s41586-025-09521-x](https://doi.org/10.1038/s41586-025-09521-x)
->
-> *these authors contributed equally and the order was decided randomly
-
-> **DeepClone, an end-to-end protocol to study somatic mutagenesis and selection at high resolution**
->
-> Ferriol Calvet, Morena Pinheiro-Santin, Erika Lopez, Raquel Blanco Martinez-Illescas, Núria Samper, Miguel L. Grau, Ferran Muiños, Rocío Chamorro González, Maria Andrianova, Federica Brando, Stefano Pellegrini, Marta Huertas, Elisabet Figuerola-Bou, Coohleen Coombes, Brendan F. Kohrn, Jeanne Fredrickson, Rosa Ana Risques, Nuria Lopez-Bigas, Abel Gonzalez-Perez
->
-> _protocols.io_ (2026) doi:[10.17504/protocols.io.dm6gp1jodgzp/v2](https://dx.doi.org/10.17504/protocols.io.dm6gp1jodgzp/v2)

assets/assess_panel/assess_panel.md

@@ -0,0 +1,44 @@
+# Experimental design principles
+
+## Definition of relationship between number of mutations detected and sequencing depth
+
+$$
+\text{Expected number of mutations} = (\text{ average depth * region of interest (bp) }) \times (\text{expected mutation density (muts/bp)})
+$$
+
+## Use cases
+
+### Predict average depth required
+
+Given a fixed panel size based on which are the regions of interest, a desired number of mutations (enough to get good selection/mutagenesis metrics) and a known estimation of the expected mutation density of the sample, you can estimate which is the sequencing depth required using this formula:
+
+$$
+\text{average depth} = \frac{\text{minimum desired mutations}}{(\text{region of interest (bp)}) \times (\text{expected mutation density (muts/bp)})}
+$$
+
+To get an estimation of the expected mutation density you can use this formula even if it is with some approximate values:
+
+$$
+\text{expected mutation density (muts/bp)} = \text{mutation rate (mutations} \times \text{bp} \times \text{year)} \times \text{time (year)}
+$$
+
+For the mutation rate this could be a possible formula; but you may also rely on previously published data.
+
+$$
+\text{mutation rate (} \frac{\text{mutations}}{\text{bp * year}} \text{)}  = \frac{\text{observed mutations}}{ \text{genome size} \times \text{age (year)} }
+$$
+
+In this example, we are assuming that the units of the mutation rate are mutations x genome x year, but this could be in any other units
+Particularly the time units could be different from year.
+
+In case that the sequencing depth is restricted by the availability of DNA or any other reason, you can decide to adjust the panel so that by increasing the depth you can still reach a stable measurement of the mutation density and other variables.
+
+An additional scenario to consider is the trinucleotide mutation probabilities. In a given sample, the different sites may have a very different mutation probablility depending on which are the active mutational processes. If these are known, you can factor in this information to adjust the requirements of the experimental design.
+
+### Expected number of mutations given an average depth
+
+The average depth represents the number of times that you sequence your region of interest, then it can be seen as:
+
+$$
+\text{Expected number of mutations} = (\text{ average depth * region of interest (bp) }) \times (\text{expected mutation density (muts/bp)})}
+$$

assets/assess_panel/assess_panel.py

@@ -0,0 +1,101 @@
+#!/usr/bin/env python
+
+
+import sys
+import pandas as pd
+import numpy as np
+import seaborn as sns
+import matplotlib.pyplot as plt
+import click
+
+
+sys.path.append("../../bin")
+
+from utils_context import triplet_context_iterator
+
+# we need a stable background mutation density estimation,
+# for this we need at least 3 mutations
+
+def mutation_based_assessment(panel_tsv, expected_mutation_density_in_mb, samples, depth, mutational_profile = None):
+    if mutational_profile is None:
+        mutational_profile = pd.DataFrame(np.ones(96) / 96)
+        mutational_profile.index = triplet_context_iterator()
+    mutational_profile = mutational_profile / mutational_profile.mean()
+    mutational_profile.columns = ["Probability"]
+
+    size_per_consequence = panel_tsv.groupby(by = ["GENE", "IMPACT", "CONTEXT_MUT"]).size().reset_index(name = "DEPTH")
+    size_per_consequence["DEPTH"] *= depth/3
+
+    size_per_consequence = size_per_consequence.merge(mutational_profile,
+                                                      left_on = "CONTEXT_MUT",
+                                                      right_index = True,
+                                                      how = "left")
+    size_per_consequence["ADJUSTED_DEPTH"] = size_per_consequence["DEPTH"] * size_per_consequence["Probability"]
+
+    total_adjusted_size = size_per_consequence.groupby(by = ["GENE", "IMPACT"])["ADJUSTED_DEPTH"].sum() / 3
+    expected_mutations_per_cnsq = (expected_mutation_density_in_mb * total_adjusted_size * samples) / 1e6
+
+    return expected_mutations_per_cnsq
+
+
+
+@click.command()
+@click.option('--panel', '-p', 'panel_path', required=True, type=click.Path(exists=True, dir_okay=False),
+              help='Path to panel TSV file')
+@click.option('--expected-mutation-density', '-e', 'expected_mutation_density', default=1.2, show_default=True,
+              type=float, help='Expected mutation density (mutations per Mb)')
+@click.option('--samples', '-s', default=100, show_default=True, type=int,
+              help='Number of samples to assume')
+@click.option('--depth', '-d', default=5000, show_default=True, type=int,
+              help='Sequencing depth (used in size adjustment)')
+@click.option('--mutational-profile', '-m', 'mutational_profile_path', default=None, type=click.Path(exists=True, dir_okay=False),
+              help='Optional mutational profile TSV with CONTEXT_MUT and Probability (or column named like all_samples.all)')
+@click.option('--output-prefix', '-o', default=None, help='Optional prefix to write TSV outputs (will create <prefix>*.tsv)')
+def main(panel_path, expected_mutation_density, samples, depth, mutational_profile_path, output_prefix):
+    """CLI entrypoint for assess_panel.
+
+    Loads the panel TSV and optional mutational profile, runs the mutation-based assessment
+    (with and without the provided profile) and optionally writes results to TSV files.
+    """
+    panel = pd.read_table(panel_path)
+
+    mutational_prof = None
+    if mutational_profile_path:
+        mutational_prof = pd.read_table(mutational_profile_path)
+        # Try to normalize common formats: expect CONTEXT_MUT as index and a column named Probability
+        if 'CONTEXT_MUT' in mutational_prof.columns:
+            # If the profile already has a Probability column, keep it; otherwise try to rename
+            if 'Probability' in mutational_prof.columns:
+                mutational_prof = mutational_prof.set_index('CONTEXT_MUT')[['Probability']]
+            elif 'all_samples.all' in mutational_prof.columns:
+                mutational_prof = mutational_prof.set_index('CONTEXT_MUT').rename(columns={'all_samples.all': 'Probability'})[['Probability']]
+            else:
+                # Fallback: use all other numeric column(s) and take the first as Probability
+                numeric_cols = mutational_prof.select_dtypes(include=[np.number]).columns.tolist()
+                if numeric_cols:
+                    mutational_prof = mutational_prof.set_index('CONTEXT_MUT')[[numeric_cols[0]]].rename(columns={numeric_cols[0]: 'Probability'})
+                else:
+                    # If no numeric columns, create uniform profile (will be normalized in function)
+                    mutational_prof = None
+
+    # Run assessments
+    mutations_per_gene_consequence = mutation_based_assessment(panel, expected_mutation_density, samples, depth)
+    mutations_per_gene_consequence_mut_prof = mutation_based_assessment(panel, expected_mutation_density, samples, depth, mutational_prof)
+
+
+    signatures_based_assessment(panel_tsv, expected_mutation_density_in_mb, samples, depth, mutational_profile)
+
+
+    # Print a short summary to stdout
+    click.echo("Mutation-based assessment (uniform profile) - top rows:")
+    click.echo(mutations_per_gene_consequence.head().to_string())
+    click.echo("\nMutation-based assessment (provided mutational profile) - top rows:")
+    click.echo(mutations_per_gene_consequence_mut_prof.head().to_string())
+
+    # Optionally write outputs
+    if output_prefix:
+        out1 = f"{output_prefix}.mutations_per_gene_consequence.tsv"
+        out2 = f"{output_prefix}.mutations_per_gene_consequence_mutprof.tsv"
+        mutations_per_gene_consequence.to_csv(out1, sep='\t')
+        mutations_per_gene_consequence_mut_prof.to_csv(out2, sep='\t')
+        click.echo(f"Wrote results to: {out1} and {out2}")