Sanity

Sampling Noise based Inference of Transcription ActivitY : Filtering of Poisson noise on a single-cell RNA-seq UMI count matrix

Single-cell RNA sequencing normalization algorithm presented in the publication Bayesian inference of gene expression states from single-cell RNA-seq data - J Breda, M Zavolan, E van Nimwegen - Nature Biotechnology, 2021.

Sanity infers the log expression levels x_gc of gene g in cell c by filtering out the Poisson noise on the UMI count matrix n_gc of gene g in cell c.

Reproducibility

The raw UMI count and normalized datasets mentioned in benchmarking in the associated publication are available on . Files are named [dataset name]_UMI_counts.txt.gz and [dataset name]_[tool name]_normalization.txt.gz.

The scripts used for running the benchmarked normalization methods and for making the figures of the preprint are in the reproducibility folder.

Input

• -f: UMI count matrix: (N_g x N_c) matrix with N_g the number of genes and N_c the number of cells. ('path/to/count_matrix.tsv')

The count matrix must be tab-separated (TSV). The tab character is the only field separator Sanity accepts: comma-separated (CSV) and space-separated files are not supported and will be misread. Convert them to TSV before running Sanity, for example with pandas.DataFrame.to_csv(path, sep='\t') or awk 'BEGIN{FS=",";OFS="\t"}{$1=$1;print}' in.csv > out.tsv. Because tabs alone delimit fields, gene and cell names may contain spaces and commas.

The first line is the header with the cell names. It may either start with a label for the gene-ID column (the same number of fields as a data row) or start straight away with the first cell name (one field fewer, as written by pandas.DataFrame.to_csv). Every data row must have exactly one gene ID followed by one value per cell; Sanity stops with an error if a row has too few or too many fields.

The count matrix can be compressed with gzip (with .tsv.gz extension).

• -f: (Alternatively) Matrix Market File Format: Sparse matrix of UMI counts. Automatically recognized by .mtx extension of the input file. Example: matrix.mtx by cellranger 2.1.0 and 3.1.0 (10x Genomics). ('path/to/text_file.mtx') Important: MTX file should be sorted by row (gene) indices! If not then you can use sort_mtx_by_row.py script in the Sanity/scripts folder to sort it. See Sorting MTX file by row (gene) indices. The MTX file can be compressed with gzip (with .gz extension).

  • -mtx_genes: (optional) Gene ID file: text file with one gene ID per line. The order of gene IDs should match the order of genes in the count matrix. Examples: genes.tsv by cellranger 2.1.0 and features.tsv by cellranger 3.1.0 (10x Genomics). ('path/to/text_file')
  • -mtx_cells: (optional) Cell ID file: text file with one cell ID per line. The order of cell IDs should match the order of cells in the count matrix. Examples: barcodes.tsv by cellranger 2.1.0 and 3.1.0 (10x Genomics). ('path/to/text_file')

• -d: (optional) Destination folder ('path/to/output/folder', default: cwd)

• -n: (optional) Number of threads (integer, default: 4). The same setting controls input indexing and inference.

  Uncompressed, row-sorted MTX files with nonnegative integral counts are indexed in parallel using bounded input buffers; no index cache is used. Counts written as 2.0 or 2e0 also qualify.

  Fractional counts, gzip input, and counts or totals outside the exact integer range of double precision use the existing serial MTX reader. If an unsupported count is encountered during parallel indexing, the input scan restarts with that reader. The indexing log reports the reader, actual thread count and elapsed time. TSV input retains its existing reader. During inference, each worker reuses its count buffer. Plain MTX rows are read through a bounded buffer with direct numeric parsing; fractional counts are supported. TSV and gzip retain their existing numeric conversion.

• --bonsai-output-mode: (optional flag, no value) Write only the numerical files and metadata required by Bonsai (NPY, MAP).

• -e: (optional) Print extended output (Boolean, 'true', 'false', '1' or '0', default: false)

• -v_m: (optional, expert-user-only) Choose the method to estimate gene-variances $v_g$. In the MAP, EAP, MLE-options, one value for $v_g$ is fixed, and the corresponding gene expression estimates are returned, in the MARG-option, the gene expression estimates are obtained by marginalizing $v_g$. The options are:

  • MAP (default): Use the maximum a posteriori estimate for $v_g$.
  • EAP: expected value of $v_g$ over the posterior.
  • MLE: maximum likelihood estimate for $v_g$.
  • MARG: (original method) The reported $v_g$ is the EAP estimate, but unlike the EAP option above, this value is not fixed when computing the gene expression estimates — instead, the gene expression estimates are obtained by marginalizing over the posterior of $v_g$.

• -vmin/-vmax: (optional, expert-user-only) Minimal and maximal considered values of the variance in log transcription quotients (double, default: v_min=0.001 v_max=50)

• -nbin: (optional, expert-user-only) Number of bins for the variance in log transcription quotients (integer, default: 160)

• -no_norm: (optional, expert-user-only) Option to skip cell size normalization (Boolean, 'true', 'false', '1' or '0', default: false)

For MAP, MLE and EAP, Sanity searches the complete variance grid first and calculates cell uncertainties only for the selected output bin. It retains compact fitting results for each bin and reconstructs the selected bin’s cell fractions and expression estimates, instead of storing cell-sized arrays for every bin. MARG still calculates uncertainties at every bin because it averages over the grid. The variance grid, bin selection and uncertainty formulas are unchanged. Library-size constants are computed once after the normalization setting is applied, and each gene’s digamma and trigamma values are computed once. Each worker also reuses its fitting workspace across bins and genes. These changes avoid repeated calculations and allocations without changing the fitting algorithm or its convergence tolerance.

Text output uses fixed six-decimal formatting. Each inference worker formats one gene row, and complete rows are written in gene order. Formatting uses the same -n worker count as inference; there is no separate writer thread setting. The number of formatted rows held in memory is bounded by the worker count.

For zero-count cells, uncertainty is calculated from the positive displacement e at which e*(2*delta+e)/(2*v) + n*log(1+f*(exp(e)-1)) = 0.5. Here n is the total count for the gene and f its fitted fraction for the cell. Sanity evaluates the full logarithm with log1p and expm1, with an overflow-safe form for large displacements. The likelihood-drop tolerance remains 1e-7. This changes zero-count uncertainty relative to the previous logarithm linearization; it retains the existing Gaussian uncertainty model and does not change the fitted means, likelihood grid or selected gene variance.

Output

• log_transcription_quotients.txt: This file contains the estimated values of the log-transcription quotients (LTQs) for each gene in each cell. The LTQ x_gc of gene g in cell c corresponds to the estimated logarithm of the fraction of mRNAs in cell c that belong to gene g. The LTQs are thus normalized such that Σ_g exp(x_gc) = 1 for each cell c. In order to get an estimate of the number of mRNAs for gene g in cell c one would thus need to multiply exp(x_gc) by the estimated total number of mRNAs M in the cell.

  GeneID	Cell 1	Cell 2	Cell 3	...
  Gene 1	-13.7227	-13.722	-13.729	...
  Gene 2	-9.96744	-10.2522	-10.1453	...
  ...	...	...	...	...

• ltq_error_bars.txt: Table with the error-bars on the estimates of the LTQs x_gc for each gene g in each cell c.

  GeneID	Cell 1	Cell 2	Cell 3	...
  Gene 1	0.630111	0.630198	0.624802	...
  Gene 2	0.315551	0.325912	0.301861	...
  ...	...	...	...	...

Extended output (optional)

• mu.txt: Estimated average LTQ μ_g of each gene g averaged over all cells.One value per line. Order corresponds to the order of genes in the geneID.txt file.

• d_mu.txt: Error bars on the inferred mean LTQs μ_g. One value per line. Order corresponds to the order of genes in the geneID.txt file.

• variance.txt: Estimated variance of the LTQs x_gc across cells c for each gene g. Note that these variances are different, and generally larger, than what one would obtain when directly calculating the variance of the estimates of x_gc from the file log_transcription_quotients.txt. This is because the estimates in this file take into account the uncertainty on the estimates of the x_gc. Thus, when estimates of true gene expression variability are needed, you are strongly advised to use the results in this file. Order corresponds to the order of genes in the geneID.txt file.

• delta.txt: Matrix of inferred log-fold changes δ_gc = x_gc-μ_g for each gene g in each cell c. TSV file, columns correspond to cell IDs from cellID.txt file, rows correspond to gene IDs from geneID.txt file.

• d_delta.txt: Matrix of error-bars for the inferred log fold-changes δ_gc. TSV file, columns correspond to cell IDs from cellID.txt file, rows correspond to gene IDs from geneID.txt file.

• likelihood.txt: This file encodes the posterior distribution of each gene's true variance in log-expression. For the numerical calculation of this distribution, the variance is a priori assumed to lie in the range [v_min,v_max], and the distribution is evaluated on a grid of N_b values spaced uniformly on a logarithmic scale. The file contains a matrix with likelihoods for each gene g and each grid point b. The first line lists the grid values v_b = v_min exp(b · Δv) for b = 0,...,N_b-1 with Δv = log(v_max/v_min)/(N_b-1), so the first and last values are v_min and v_max themselves; the following lines correspond to the likelihood values of genes.

  Variance	0.001000	0.001070	0.001146	...
  Gene 1	0.018	0.019	0.020	...
  Gene 2	0.0006	0.0051	0.0031	...
  ...	...	...	...	...

Bonsai output mode (optional)

Use --bonsai-output-mode to write only:

• delta.npy and d_delta.npy: float32 arrays of shape (genes, cells), containing deviations and their standard deviations (not variances).
• mu.npy and variance.npy: float64 arrays of shape (genes, 1), containing one baseline mean and one prior variance per gene.
• geneID.txt, cellID.txt and sanity_command.txt: identifiers and version, method, output format and command metadata.

./bin/Sanity -f counts.tsv -d bonsai_output -n 12 --bonsai-output-mode

This flag requires MAP inference. It selects these seven files regardless of -e; it cannot be combined with gzip or the experimental NPY output flag. Use a folder without previous numerical text output; Sanity reports an error if conflicting files are present. Ordinary and extended output remain text.

sanity_command.txt records # Output mode: bonsai so consumers can select the NPY reader without another command-line option.

Arrays use NPY version 1.0, little-endian values and C order. Row order matches geneID.txt; matrix columns match cellID.txt. Inference and uncertainty calculations remain double precision. Workers prepare float32 matrix rows, then write them in gene order with memory bounded by the thread count. Gene vectors retain double precision. Unlike six-decimal text, NPY values are not rounded to a fixed number of decimal places.

The consuming Bonsai reader must support this format; text-only readers cannot use these outputs. In Python, arrays can be opened without loading the whole matrix using numpy.load(path, mmap_mode="r", allow_pickle=False).

Usage

  ./Sanity <option(s)>
  Options:
    -h, --help                          Show this help message
    -v, --version                       Show the current version
    -f, --file                          Specify the input transcript count text file (.mtx for Matrix Market File Format)
    -mtx_genes, --mtx_gene_name_file    Specify the gene name text file (only needed if .mtx input file)
    -mtx_cells, --mtx_cell_name_file    Specify the cell name text file (only needed if .mtx input file)
    -d, --destination                   Specify the destination path (default: pwd)
    -n, --n_threads                     Specify the number of threads to be used (default: 4)
    --bonsai-output-mode               Write only the four numerical files and IDs/metadata required by Bonsai (NPY, MAP)
    -e, --extended_output               Option to print extended output (default: false, choice: false,0,true,1)
    -v_m, --v_method                    Option to specify the method for variance estimation (default: MAP, choice: MAP, EAP, MLE, MARG)
    -vmin, --variance_min               Minimal value of variance in log transcription quotient (default: 0.001)
    -vmax, --variance_max               Maximal value of variance in log transcription quotient (default: 50)
    -nbin, --number_of_variance_bins    Number of bins for the variance in log transcription quotient  (default: 160)
    -no_norm, --no_cell_size_normalization  Option to skip cell size normalization (default: false, choice: false,0,true,1)

Installation

Text formatting requires floating-point std::to_chars support in the C++17 standard library.

• Clone the GitHub repository

git clone https://github.com/jmbreda/Sanity.git

• Install OpenMP library
  • On Linux
    If not already installed (Check with ldconfig -p | grep libgomp, no output if not installed), do

   sudo apt-get update
   sudo apt-get install libgomp1
  

  • On mac OS
    Apple's own g++/clang++ does not support OpenMP, so you need a real GCC.

    • Using brew (recommended)

     brew install gcc
    

    The Makefile detects this automatically and picks the right g++ binary, no further changes needed.

    • Using macports

     port install gcc13
    

           Tell the Makefile to use it by compiling with make CC=g++-mp-13 (adjust the version number to the package you installed).

• Move to the source code directory and compile.

cd Sanity/src
make

• The binary file is located in

Sanity/bin/Sanity

• Alternatively, pre-built binaries are published on the repository's Releases page.

Sanity_distance

Compute cell-cell distances from Sanity output files. Needs extended outputs of Sanity (-e 1 option).

Input

• -f: The output folder of the Sanity run, specified with the -d option in Sanity ('path/to/folder')
• -s2n: (optional) The gene signal to noise ratio used as gene cut-off (double, default: 1.0)
• -err: (optional) Compute cell-cell distance with or without error bars (boolean, default: 1 or true)
• -n: (optional) Number of threads (integer, default: 4)

Output

• Cell-cell distance: (N_c(N_c-1)/2) vector of cell to cell distances dist(cell_i,cell_j), i=1,...,N_c-1, j=i+1,...,N_c, with N_c the number of cells.

  dist(cell1,cell2)
  dist(cell1,cell3)
  dist(cell1,cell4)
  ...
  dist(cellNc-2,cellNc-1)
  dist(cellNc-2,cellNc)
  dist(cellNc-1,cellNc)

located in the Sanity output folder (specified with -f option), named cell_cell_distance_[...].txt, depending on the -err and -s2n options.

Usage

./Sanity_distance <option(s)>
Options:
	-h,--help		Show this help message
	-v,--version		Show the current version
	-f,--folder		Specify the input folder with extended output from Sanity
	-s2n,--signal_to_noise_cutoff	Minimal signal/noise of genes to include in the distance calculation (default: 1.0)
	-err,--with_error_bars	Compute cell-cell distance taking the errorbar epsilon into account (default: true)
	-n,--n_threads		Specify the number of threads to be used (default: 4)

Installation

Same dependencies as Sanity (see above).

• Move to the source code directory and compile.

cd Sanity/src
make Sanity_distance

• The binary file is located in

Sanity/bin/Sanity_distance

Sorting MTX file by row (gene) indices

Use script sort_mtx_by_row.py in the Sanity/scripts folder to sort the MTX file by row (gene) indices. The script takes as input the MTX file. The output is a sorted MTX file by default saved as sorted_[original filename].gz in the same folder as the input file. User can specify the output file name and path with the -o option. The saved MTX file is compressed with gzip.

Requirements: Python3, sort

Example usage:

python3 sort_mtx_by_row.py -i input.mtx

Tests

A simple test comparing Sanity's output across the -v_m methods (MAP, EAP, MLE, MARG) on a small count matrix is available in the tests directory. See tests/README.md for details on running it.

Help

For any questions or assistance regarding Sanity, please post your question in the issues section.