cellsnp-lite output not working with Vireo - reporting really high read depths

[dannyconrad]

I have a scATACseq dataset where I'm trying to use the mitochondrial reads to demultiplex two mouse strains. When I run cellSNP and cellsnp-lite in mode 2 using what I believe to be matching parameters and the same input files, both identify all of the same SNPs, but the read depths reported in cellSNP.base.vcf.gz and the sparse matrices are different; the depths reported by cellsnp-lite are ~100X greater. When I run vireo on these two different outputs, the cellSNP results seem to properly demultiplex my donors, whereas the cellsnp-lite results leave all cells unassigned. Importantly, I noticed in the donor_ids.tsv file for the cellsnp-lite vireo run, that all cells are reported to have the same results:

cell	donor_id	prob_max	prob_doublet	n_vars	best_singlet	best_doublet	doublet_logLikRatio
AAACGAAAGGTCGTTT-1	unassigned	4.76e-01	4.79e-02	370	donor0	donor0,donor1	0.000

I guess I'm not sure if this is a bug with cellsnp-lite or if I did something drastically different running the two programs without realizing it. Any ideas?

Here are the two batches of code I ran and the relevant console outputs:

cellsnp-lite

cellsnp-lite \
-s ./chrM_CB.bam \
-b ./barcodes_mt_revcomp.tsv \
-O ./cellsnp-lite_DeNovo \
-p 40 \
--minMAF 0.1 \
--UMItag None \
--chrom M
--exclFLAG UNMAP,SECONDARY,QCFAIL
[I::main] start time: 2021-11-24 11:49:27
[W::check_args] Max depth set to maximum value (2147483647)
[I::main] mode 2a: pileup 1 whole chromosomes in 4787 single cells.
[W::csp_pileup_core] Combined max depth is above 1M. Potential memory hog!
[I::csp_pileup_core][Thread-0] processing chrom M ...
[I::csp_pileup_core][Thread-0] has pileup-ed in total 370 SNPs for chrom M
[I::main] All Done!
[I::main] end time: 2021-11-24 13:12:30
[I::main] time spent: 4983 seconds.

less cellsnp-lite_DeNovo/cellSNP.base.vcf.gz | head 
##fileformat=VCFv4.2
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
M	55	.	T	G	.	PASS	AD=146211;DP=489229;OTH=247
M	165	.	C	A	.	PASS	AD=142350;DP=607292;OTH=268
M	348	.	T	C	.	PASS	AD=103804;DP=350149;OTH=251

vireo \
-c ./cellsnp-lite_DeNovo/ \
-N 2 \
-o ./vireo/cellsnp-lite_DeNovo \
-p 40
[vireo] Loading cell folder ...
[vireo] Demultiplex 4787 cells to 2 donors with 370 variants.
[vireo] lower bound ranges [-114928564.4, -114928564.4, -114928564.4]
[vireo] allelic rate mean and concentrations:
[[0.01  0.284 0.99 ]]
[[5.00000000e+01 2.39958576e+08 5.00000000e+01]]
[vireo] donor size before removing doublets:
donor0	donor1
2394	2394
[vireo] final donor size:
unassigned
4787
[vireo] All done: 0 min 16.0 sec

cellSNP

cellSNP \
-s ./chrM_CB.bam \
-b ./barcodes_mt_revcomp.tsv \
-O ./cellSNP_DeNovo \
-p 40 \
--minMAF=0.1 \
--UMItag=None \
--chrom=M \
--maxFLAG=4096
[cellSNP] mode 2: pileup 1 whole chromosomes in 4787 single cells.
[cellSNP] Whole genome pileupped, now merging all variants ...
[cellSNP] 408 lines in final vcf file
[cellSNP] All done: 64 min 43.5 sec

less cellSNP_DeNovo/cellSNP.base.vcf.gz | head 
##fileformat=VCFv4.2
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
chrM	55	.	T	G	.	PASS	AD=1924;DP=7534;OTH=2
chrM	165	.	C	A	.	PASS	AD=2100;DP=7740;OTH=0
chrM	348	.	T	C	.	PASS	AD=2749;DP=7402;OTH=2

vireo \
-c ./cellSNP_DeNovo/cellSNP.cells.vcf.gz \
-N 2 \
-o ./vireo/cellSNP_DeNovo \
-p 40
[vireo] Loading cell VCF file ...
[vireo] Demultiplex 4787 cells to 2 donors with 370 variants.
[vireo] lower bound ranges [-1585667.0, -1585667.0, -430338.7]
[vireo] allelic rate mean and concentrations:
[[0.038 0.907 0.938]]
[[1857003.8  450325.5  418458.7]]
[vireo] donor size before removing doublets:
donor0	donor1
1735	3052
[vireo] final donor size:
donor0	donor1	doublet	unassigned
1581	2891	286	29
[vireo] All done: 0 min 15.2 sec

[hxj5]

Hi, thanks for your feedback. Could you share the version of cellSNP and cellsnp-lite?

The reason why the read counts are different: cellSNP (actually the dependent pysam) has a default limitation that the max_depth (i.e., max pileup-ed read count) is 8000, so you may check that every DP<=8000 in the cellSNP output. But cellsnp-lite does not have this max_depth limitation, it will pileup as many reads as possible. (we are trying to add a --maxPileup option with similar function as max_depth, for the next release of cellsnp-lite)

The reason why cellsnp-lite+vireo does not work well: The large read counts in cellsnp-lite output makes it more likely for vireo to reach local optima so that the parameters of two donors become the same and hence vireo cannot assign the cells to certain donor.

Besides, vireo is designed for nuclear SNVs. For mito SNVs, you may want to try this tutorial, which was used by MQuad. Note that the duplicate reads should probably be removed beforehand, as there are not UMIs in your case.

Xianjie