Compute sample saturation curve by downsampling

TODOs:

• check if supplying the .bam file with unaligned reads makes a difference
• handle different .bam files - with mapped and unmapped reads - better reporting

This script subsamples the .bam file and fits the saturation curve based on the amount of input reads. For the estimation to be precise, the mapping rate and the number of cells should be supplied.
Example toy dataset: 1000 cells, 80% mapping rate Downsample the reads, extract the tags and compute saturation stats (use -c to use more cores):

python saturation_table.py --ncpu 5 -b test/sample.bam -n 1000 -r 0.8 -o output.tsv

CAVE: this creates a big tabular file

Fit the MM model, predict the number of input reads for --target saturation and plot:

python scripts/plot_curve.py  output.tsv saturation.png --target 0.7 

output.tsv - contains the sequencing saturation statistics for to 10 (-s) subsampling steps
saturation.png - contains the plot of the saturation curve

It's useful to examine residuals plot to see if the model tends to over or underestimate the coverage needeed.

How to compute saturation?

The saturation is defined as following:

Sequencing saturation is a measure of the fraction of library complexity that was sequenced in a given experiment. The inverse of one minus the sequencing saturation can be interpreted as the number of additional reads it would take to detect a new transcript.

How is it calculated?

Vis 10x reference for details.

The web_summary.html output from cellranger count includes a metric called "Sequencing Saturation". This metric quantifies the fraction of reads originating from an already-observed UMI. More specifically, this is the fraction of confidently mapped, valid cell-barcode, valid UMI reads that are non-unique (match an existing cell-barcode, UMI, gene combination).

The script attampts to mimic the cellranger's approach by extracting the relevant tags from the .bam file. Note that cellranger is computing saturation by considering confidently mapped reads which may inroduce some discrepancies with other methods.

The formula for calculating this metric is as follows:

$$ \text{Sequencing Saturation} = 1 - \frac{N_{\text{dedup}}}{N_{\text{dedup}} + N_{\text{dup}}}$$