search

cortes-ciriano-lab / SComatic

A tool for detecting somatic variants in single cell data
Other
145 stars 20 forks source link

readme

SComatic

SComatic is a tool that provides functionalities to detect somatic single-nucleotide mutations in high-throughput single-cell genomics and transcriptomics data sets, such as single-cell RNA-seq and single-cell ATAC-seq.

If you use SComatic (see License at the bottom of this page), please cite our publication Muyas et al. 2023.

For further details on SComatic, its assumptions, limitations and applications, please see Muyas et al. 2023.

Algorithm

Created with BioRender.com

Installation and requirements

SComatic requires Python version >=3.7.0, R version >=3.6.0, samtools, bedtools and datamash (>=v1.1.0, only for building your own panel of normals).

Detection of somatic mutations in single-cell data sets using SComatic

We show below how to run SComatic for the detection of somatic mutations in scRNA-seq data. SComatic requires two data types as input:

SComatic consists of the following 4 steps, each of which is run using a different Python script as indicated below.

Step 1: Splitting alignment file into cell-type-specific bams

The first step consists of splitting the BAM file containing aligned sequencing reads for all cell types in a sample into cell-type-specific BAM files using precomputed cell type annotations. The BAM file must contain the cell type barcode information in the cell barcode tag “CB” (as reported by popular tools, such as Cell Ranger, 10x Genomics). In addition, we strongly suggest generating the input bam files to contain the "nM" tag (number of mismatches) and the "NH" tag (number of hits), which will be used to remove low-quality reads for downstream analysis (--max_nM and --max_NH parameters. However, we included the possibility of running SComatic without these last two tags (nM and NH), at the user's own risk.

Step 1 is executed using the script SplitBam/SplitBamCellTypes.py, which has the following parameters:

Split alignment file into cell type specific BAMs

optional arguments: -h, --help show this help message and exit --bam BAM BAM file to be analysed (Sorted by coordinate) --meta META Metadata file mapping cell barcodes to cell type information --id ID Sample ID --max_nM MAX_NM Maximum number of mismatches permitted to consider reads for analysis. By default, this filter is switched off, although we recommed using --max_nM 5. If applied, this filter requires having the nM tag in the bam file. [Default: Switched off] --max_NH MAX_NH Maximum number of alignment hits permitted to consider reads for analysis. By default, this filter is switched off, although we recommend using --max_NH 1. This filter requires having the NH tag in the bam file. [Default: Switched off] --min_MQ MIN_MQ Minimum mapping quality required to consider reads for analysis. Set this value to 0 to switch this filter off. --min_MQ 255 is recommended for RNA data, and --min_MQ 30 for DNA data. [Default: 255] --n_trim N_TRIM Number of bases trimmed by setting the base quality to 0 at the beginning and end of each read [Default: 0] --outdir OUTDIR Out directory


The precomputed cell type annotation file provided with the --meta parameter must contain at least the following two columns (Index for cell barcode ID and Cell_type for the precomputed cell type annotation) and must be a tab-separated file. Cell type annotations containing whitespaces or any of the following special characters (~ . ` ! @ # $ % ^ & * ( ) { | } / \ : ; " ' < > ? , = +) are not supported. Dashes and underscores are supported. Whitespace characters in the filenames are not supported. **Importanly, do not include doublets in this file or cell types with unkown cell type annotations.**

Index Cell_type AAACCTGCATGCTAGT Epithelial AAACCTGGTAGCCTAT Epithelial AAACCTGGTTGTCGCG Epithelial AAACCTGTCATGTGGT Epithelial AAACCTGTCCTTGGTC Epithelial AAACCTGTCGGATGTT T_cell AAACCTGTCGTACGGC T_cell AAACCTGTCTTGCAAG T_cell AAACGGGAGACGCACA T_cell


In addition to the cell-type specific bam files, this script creates a txt (\*.report.txt) file showing the number of reads filter, and importantly, why there were filtered out. 

- **Example:** check [here](/docs/SComaticExample.md) to see how to run this step with an example sample.  

## Step 2: Collecting base count information
Base count information for each cell type and for every position in the genome is recorded in a base count matrix indexed by cell types and genomic coordinates. 

The command line to run this step is: 

- List of parameters:

python scripts/BaseCellCounter/BaseCellCounter.py --help usage: BaseCellCounter.py [-h] --bam BAM --ref REF --chrom CHROM [--out_folder OUT_FOLDER] [--id ID] [--nprocs NPROCS] [--bin BIN] [--bed BED] [--bed_out BED_OUT] [--min_ac MIN_AC] [--min_af MIN_AF] [--min_dp MIN_DP] [--min_cc MIN_CC] [--min_bq MIN_BQ] [--min_mq MIN_MQ] [--tmp_dir TMP_DIR]

Script to obtain a list of base and cell counts in scRNA bam file

optional arguments: -h, --help Show this help message and exit --bam BAM Tumor bam file to be analysed --ref REF Reference genome. fai must be available in the same folder as reference --chrom CHROM Chromosome to be analysed. Use --chrom all to analyse all chromosomes --out_folder OUT_FOLDER Output folder --id ID Prefix of out file. If provided, please use the following format: .[cell_type] . Example: sample1.t_cell. If not provided, it is taken from the bam file --nprocs NPROCS Number of processes [Default = 1] --bin BIN Bin size for running the analysis [Default 50000] --bed BED Regions to focus the analysis on. Three-column bed file --bed_out BED_OUT Regions to ignore. Three-column bed file --min_ac MIN_AC Minimum alternative count to consider a genomic site for further analysis. Default = 0 --min_af MIN_AF Minimum alternative allele fraction to consider a genomic site for further analysis. Default = 0 --min_dp MIN_DP Minimum coverage to consider a genomic site for further analysis. Default = 5 --min_cc MIN_CC Minimum number of cells required to consider a genomic site for further analysis. Default = 5 --min_bq MIN_BQ Minimum base quality permited for the base counts. Default = 20 --min_mq MIN_MQ Minimum mapping quality required to analyse read. Default = 255 --tmp_dir TMP_DIR Temporary folder for tmp files


- **Example:** check [here](/docs/SComaticExample.md) to see how to run this step with an example sample.  

## Step 3: Merging  base count matrices
In Step 3, SComatic takes as input base count matrices computed in Step 2 for all cell types analysed to merge them into a single base count matrix, which is stored in tsv format. Individual base count matrices to be merged need to be stored in the same directory.

- List of parameters:
```python 
python scripts/MergeCounts/MergeBaseCellCounts.py --help
usage: MergeBaseCellCounts.py [-h] --tsv_folder TSV_FOLDER --outfile
                                       OUTFILE

Script to merge the cell/base counts tsv files per cell type in only one file

optional arguments:
  -h, --help            Show this help message and exit
  --tsv_folder TSV_FOLDER
                        Folder with cell/base count files (tsv) per cell type.
                        Avoid not desired tsv files in this folder
  --outfile OUTFILE   Out file name

Step 4: Detection of somatic mutations

The last step consists of running two scripts to call somatic mutations.

Step 4.1

SComatic applies a set of hard filters and Beta binomial tests to discount sites affected by recurrent technical artefacts as somatic mutations.

Script to perform the scRNA somatic variant calling

optional arguments: -h, --help Show this help message and exit --infile INFILE Input file with all samples merged in a single tsv --outfile OUTFILE Output file prefix --ref REF Reference fasta file (*fai must exist) --min_cov MIN_COV Minimum depth of coverage to consider a sample. [Default: 5] --min_cells MIN_CELLS Minimum number of cells with sequencing depth at a site to consider a genomic site for further analysis. [Default: 5] --min_ac_cells MIN_AC_CELLS Minimum number of cells supporting the alternative allele to call a mutation. [Default: 2] --min_ac_reads MIN_AC_READS Minimum number of reads supporting the alternative allele to call a mutation. [Default: 3] --max_cell_types MAX_CELL_TYPES Maximum number of cell types carrying a mutation to make a somatic call. [Default: 1] --min_cell_types MIN_CELL_TYPES Minimum number of cell types with enough coverage across enough cells to consider a site as callable [Default: 2] --fisher_cutoff FISHER_CUTOFF P value cutoff for the Fisher's exact test performed to detect strand bias. A float value is expected, if applied, we recommend 0.001. By default, this test is switched off with a value of 1 [Default: 1] --alpha1 ALPHA1 Alpha parameter for Beta distribution of read counts. [Default: 0.260288007167716] --beta1 BETA1 Beta parameter for Beta distribution of read counts. [Default: 173.94711910763732] --alpha2 ALPHA2 Alpha parameter for Beta distribution of cell counts. [Default: 0.08354121346569514] --beta2 BETA2 Beta parameter for Beta distribution of cell counts. [Default: 103.47683488327257]


To estimate new Beta binomial parameters whenever required by the user, SComatic provides the following [scripts](/docs/betabinomialestimation.md).

- **Example:** check [here](/docs/SComaticExample.md) to see how to run this step with an example sample.  

### Step 4.2 
Scomatic takes the output of the previous step (4.1) and applies additional filters based on external datasets (RNA editing and Panel of Normals), and flags clustered mutations. High quality mutations are marked with the label “PASS” in the FILTER column of the output file. When using the provided [RNA editing](RNAediting) file, please cite the following articles: [Tan et al. 2017](https://www.ncbi.nlm.nih.gov/pubmed/29022589), [Kiran et al. 2012](https://www.ncbi.nlm.nih.gov/pubmed/23074185) and [Picardi et al. 2017](https://www.ncbi.nlm.nih.gov/pubmed/27587585).

- List of parameters: 
```python 
python scripts/BaseCellCalling/BaseCellCalling.step2.py --help
usage: BaseCellCalling.step2.py [-h] --infile INFILE --outfile
                                         OUTFILE [--editing EDITING]
                                         [--pon PON]
                                         [--min_distance MIN_DISTANCE]

Script to perform the scRNA somatic variant calling

optional arguments:
  -h, --help            Show this help message and exit
  --infile INFILE       Input file with all samples merged in a single tsv
  --outfile OUTFILE     Output file prefix
  --editing EDITING     RNA editing file to be used to remove RNA-diting sites
  --pon PON             Panel of normals (PoN) file to be used to remove
                        germline polymorphisms and recurrent artefacts
  --min_distance MIN_DISTANCE
                        Minimum distance allowed between potential somatic
                        variants [Default: 5]

Estimating new beta-binomial parameters

SComatic models the background error rate of the technology used to generate the single-cell data (e.g., single-cell RNA-seq) using a Beta binomial distribution. Specifically, non-reference allele counts at homozygous reference sites are modelled using a binomial distribution with parameter P (error rate), which is a random variable that follows a Beta distribution with parameters α and β. Default values for the Beta binomial tests used in Step 4.1 are computed using the data sets described in the manuscript. However, we provide scripts to allow the user to reparameterize the Beta binomial using other data sets.

Generating a custom Panel of Normals

In Step 4.2, SComatic uses a Panel of Normals (PoN) to detect systematic errors and germline contamination in the somatic mutation callset. The PoN provided in this repository is computed using the data described in the manuscript (Hg38 reference genome). However, SComatic provides a script to build a custom PoN using other data sets if required.

Other SComatic functionalities

SComatic provides the following additional functionalities, which are described in detail here.

FAQs - Frequently asked questions

This section answers some of the users' most recurrent doubts when running SComatic.

  1. Are the SComatic parameters for scATAC-seq data the same as for scRNA-seq data?
  2. How can we perform the variant annotation with the SComatic output?
  3. Can SComatic work with other types of PoN files?
  4. Can we use the calls from other callers to genotype unique cells using SComatic?
  5. How do different cell type labels (e.g. different levels of granularity) affect the SComatic performance?
  6. What does it happen if CellRanger does not properly trim all non-genomic sequences (adapters) from the reads?
  7. How to interpret the SingleCellGenotype.py output?

Contact

If you have any comments or suggestions about SComatic please raise an issue or contact us:

Francesc Muyas: fmuyas@ebi.ac.uk

Isidro Cortes-Ciriano: icortes@ebi.ac.uk

License

SComatic is free for academic use only. If you are not a member of a public funded academic and/or education and/or research institution you must obtain a commercial license from EMBL Enterprise Management GmbH (EMBLEM); please email EMBLEM (info@embl-em.de).