search

single-cell-genetics / cellsnp-lite

Efficient genotyping bi-allelic SNPs on single cells
https://cellsnp-lite.readthedocs.io
Apache License 2.0
124 stars 11 forks source link

0 variants found #129

Open Sun-storm opened 1 week ago

Sun-storm commented 1 week ago

Hi!

I'm having an issue running cellSNP. When I do it, it runs smoothly but when I try to run VireoSNP with the results I get this error message:

$ vireo -c /a/b/c/demultiplexing/cellsnp_results/second_try -N 8 -o /a/b/c/demultiplexing/vireoSNP/first_try [vireo] Loading cell folder ... [vireo] Demultiplex 25508 cells to 8 donors with 0 variants. Traceback (most recent call last): File "/services/tools/anaconda3/2023.09-0/bin/vireo", line 8, in <module> sys.exit(main()) ^^^^^^ File "/services/tools/anaconda3/2023.09-0/lib/python3.11/site-packages/vireoSNP/vireo.py", line 204, in main res_vireo = vireo_wrap(cell_dat['AD'], cell_dat['DP'], n_donor=n_donor, ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/services/tools/anaconda3/2023.09-0/lib/python3.11/site-packages/vireoSNP/utils/vireo_wrap.py", line 86, in vireo_wrap _models_all[im].fit(AD, DP, min_iter=5, max_iter=max_iter_init, File "/services/tools/anaconda3/2023.09-0/lib/python3.11/site-packages/vireoSNP/utils/vireo_model.py", line 313, in fit ELBO += np.sum(get_binom_coeff(AD, DP)) ^^^^^^^^^^^^^^^^^^^^^^^ File "/services/tools/anaconda3/2023.09-0/lib/python3.11/site-packages/vireoSNP/utils/vireo_base.py", line 18, in get_binom_coeff binom_coeff = np.log(binom(_DP, _AD)) ^^^^^^^^^^^^^^^ TypeError: ufunc 'binom' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

I looked at the DP and AD files and they show this

DP file

%%MatrixMarket matrix coordinate integer general % 0 25508 0

AD file

%%MatrixMarket matrix coordinate integer general % 0 25508 0

I'm thus guessing that CellSNP is not finding any variants at all, which doesn't make sense to me as I'm applying the suggested common human snps file (hg38) to human samples. Here is the code I use to run CellSNP

cellsnp-lite -s /a/b/c/demultiplexing/base_data/four_lanes_merge.bam /a/b/c/demultiplexing/base_data/four_lanes_merge.bam.bai -b /a/b/c/demultiplexing/base_data/merged_barcodes.tsv -O /a/b/c/demultiplexing/cellsnp_results/second_try -R /a/b/c/demultiplexing/base_data/common_human_snps/genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf.gz -p 10 --minMAF 0.1 --minCOUNT 20 —gzip

Any suggestions for what might be the problem?

Thank you!! :)

hxj5 commented 1 week ago

Hi, the "-s" option only requires the BAM file(s), the ".bai" index file should not be specified. You may re-run cellsnp-lite without the ".bai" file. If there are still no SNPs outputed, you may checkout the manual FAQ: Why the output files are empty, no SNPs genotyped?.

Sun-storm commented 1 week ago

Thank you Xianjie. I first tried without the .bai and it gave me an error at the time, while adding it solved it. Now though I see it works without it as well. Not sure why that is. I'm also adding adding the tag --UMItag None, as it is scATAC data.

I will post if this solves the problem :)

Sun-storm commented 1 day ago

Hi again :)

After multiple attempts, I'm still unable to get any hits.

I removed the .bai, attached the --UMItag None tag and tried everything again with the big SNP file, genome1K.phase3.SNP_AF5e4.chr1toX.hg38.vcf.gz. I also changed the name from the .bam file from chrM to chrMT, to match the VCF file. I'm still getting this:

%%MatrixMarket matrix coordinate integer general % 0 25508 0

The only thing I've failed at is checking for unmatched cell barcodes. I don't know how to do that. Are you aware of any reasorces to learn to do this?

Thank you for your help ^^

hxj5 commented 19 hours ago

Hi, you may share more details here for debugging, e.g.,

Sun-storm commented 10 hours ago

Wonderful.

The command line is: $ cellsnp-lite -s a/b/c/four_lanes_merge_reheadered.bam -b a/b/c/merged_barcodes.tsv -O a/b/e/eigth_try_2 -R a/b/c/d/genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf.gz -p 10 --minMAF 0.1 —gzip --UMItag None

The logging file is:

[I::main] start time: 2024-06-27 12:02:37 [W::check_args] Max depth set to maximum value (2147483647) [W::check_args] Max pileup set to maximum value (2147483647) [I::main] global settings after checking: num of input files = 1 out_dir = a/b/e/eigth_try_2 is_out_zip = 0, is_genotype = 0 is_target = 0, num_of_pos = 0 num_of_barcodes = 25508, num_of_samples = 0 refseq file = (null) 0 chroms: cell-tag = CB, umi-tag = (null) nthreads = 10, tp_max_open = 32768 mthreads = 10, tp_errno = 0, tp_ntry = 0 min_count = 20, min_maf = 0.10, double_gl = 0 min_len = 30, min_mapq = 20 rflag_filter = 1796, rflag_require = 0 max_depth = 2147483647, max_pileup = 2147483647, no_orphan = 1 [I::main] loading the VCF file for given SNPs ... [I::main] fetching 7352497 candidate variants ... [I::main] mode 1a: fetch given SNPs in 25508 single cells. [I::csp_fetch_core][Thread-1] 2.00% SNPs processed. [I::csp_fetch_core][Thread-2] 2.00% SNPs processed. [I::csp_fetch_core][Thread-4] 2.00% SNPs processed. [I::csp_fetch_core][Thread-8] 2.00% SNPs processed. [I::csp_fetch_core][Thread-5] 2.00% SNPs processed. [I::csp_fetch_core][Thread-9] 2.00% SNPs processed. [I::csp_fetch_core][Thread-7] 2.00% SNPs processed. [I::csp_fetch_core][Thread-6] 2.00% SNPs processed. [I::csp_fetch_core][Thread-3] 2.00% SNPs processed. [I::csp_fetch_core][Thread-0] 2.00% SNPs processed. [I::csp_fetch_core][Thread-1] 4.00% SNPs processed. [I::csp_fetch_core][Thread-8] 4.00% SNPs processed. ... [I::csp_fetch_core][Thread-2] 98.00% SNPs processed. [I::csp_fetch_core][Thread-5] 96.00% SNPs processed. [I::csp_fetch_core][Thread-1] 98.00% SNPs processed. [I::csp_fetch_core][Thread-8] 96.00% SNPs processed. [I::csp_fetch_core][Thread-0] 98.00% SNPs processed. [I::csp_fetch_core][Thread-3] 98.00% SNPs processed. [I::csp_fetch_core][Thread-7] 98.00% SNPs processed. [I::csp_fetch_core][Thread-6] 98.00% SNPs processed. [I::csp_fetch_core][Thread-4] 98.00% SNPs processed. [I::csp_fetch_core][Thread-5] 98.00% SNPs processed. [I::csp_fetch_core][Thread-8] 98.00% SNPs processed. [I::main] All Done! [I::main] Version: 1.2.3 (htslib 1.20) [I::main] CMD: cellsnp-lite -s a/b/c/four_lanes_merge_reheadered.bam -b a/b/c/merged_barcodes.tsv -O /a/b/e/eigth_try_2 -R a/b/c/d/genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf.gz -p 10 --minMAF 0.1 --UMItag None —gzip [I::main] end time: 2024-06-27 13:04:46 [I::main] time spent: 3729 seconds.

The head of the .bam is:

A00642:443:HJFT7DRX3:2:2272:28348:4022 121 chr1 9990 0 58M 9990 0 TGTCTGACTGATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACC FFFFFFFFFFFFFFFF,FFF,:FFFF:FFFFFFF:FFFFFFFFFFFFFFFF,FFFFF: NM:i:2 MD:Z:6T1C49 AS:i:49 XS:i:48 XA:Z:chr22,+50808422,48M10S,0; CR:Z:ATTTGTCAGTGAGTGC CY:Z:F,FFF,FFFFF,F,FF CB:Z:ATTTGTCAGTGAGTGC-1 BC:Z:GTGGAGCA QT:Z:FFFF::FF RG:Z:CtrlMix_D1:MissingLibrary:1:HJFT7DRX3:2 A00642:443:HJFT7DRX3:2:2272:28348:4022 181 chr1 9990 0 9990 0 TCAGTAGGAAAACAGGTCATGCCATTCAATTTAATTATTTTCGTCAGATC FFFFFFFFFFFFFFFF,FF:F,FFFFFFFFFFFFF,FFF:FF,FFFFFFF AS:i:0 XS:i:0 CR:Z:ATTTGTCAGTGAGTGC CY:Z:F,FFF,FFFFF,F,FF CB:Z:ATTTGTCAGTGAGTGC-1 BC:Z:GTGGAGCA QT:Z:FFFF::FF RG:Z:CtrlMix_D1:MissingLibrary:1:HJFT7DRX3:2 A00642:443:HJFT7DRX3:2:2169:5945:20713 113 chr1 9995 10 58M chr3 69774748 0 TTCCGATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC FFF,:FF:FF:FFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFF NM:i:1 MD:Z:0G57 AS:i:57 XS:i:53 XA:Z:chr22,+50808417,53M5S,0; CR:Z:GTGTCCTGTAGAAAGG CY:Z:F:F:FFFFFFFFFFFF CB:Z:GTGTCCTGTAGAAAGG-1 BC:Z:GATGCAGT QT:Z:FFFFFFFF RG:Z:CtrlMix_C1:MissingLibrary:1:HJFT7DRX3:2 A00642:443:HJFT7DRX3:2:2107:5186:8093 129 chr1 9996 0 3M1D37Mchr13 40154541 0 TCCATAACCCTAACCCTAACCCTAACCCTAACCCTAACCC FFF,FFFFFFFFFFFFFFFFFFFF:FFF::FFFFFFFFFF NM:i:1 MD:Z:3^G37 AS:i:37 XS:i:37CR:Z:CGTTGCACAAAGCATA CY:Z:F,FFFFFF:FF:FF,F CB:Z:CCTTGCACAAAGCATA-1 BC:Z:GTGGAGCA QT:Z:,FFFFF,F TR:Z:CTGTCTCTTA TQ:Z:FFFFF:F,FF RG:Z:CtrlMix_D1:MissingLibrary:1:HJFT7DRX3:2 A00642:443:HJFT7DRX3:1:2269:18810:16376 81 chr1 9996 0 58M chr10 42284224 0 TCCTATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACC FFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:1 MD:Z:3G54 AS:i:54 XS:i:54 CR:Z:TTACCCGTCCCACTAC CY:Z:FFFFFFFFFFFFFFFF CB:Z:TTACCCGTCCCACTAC-1 BC:Z:CCGAGGCA QT:Z:FFFFFFFF RG:Z:CtrlMix_C1:MissingLibrary:1:HJFT7DRX3:1 A00642:443:HJFT7DRX3:1:2269:13449:27727 81 chr1 9996 0 3M1D55Mchr9 38147954 0 TCCATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCC FFFFF:,F:FFF,:FFF:FFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFF:F,FFF NM:i:1 MD:Z:3^G55 AS:i:55 XS:i:55 CR:Z:CGGACTGTCCACGAGC CY:Z:FFFFFFFFFFFFFFFF CB:Z:CGGACTGTCCACGAGC-1 BC:Z:GTTTGCCT QT:Z:FF:FFFFF RG:Z:CtrlMix_E1:MissingLibrary:1:HJFT7DRX3:1 A00642:443:HJFT7DRX3:2:2205:4933:20055 121 chr1 9997 9 44M 9997 0 GCGATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC ::,FFFFFFFFFFFFFFFFF:FFFFF:FFFFFFFFFF:FFFFFF NM:i:1 MD:Z:0C43 AS:i:43 XS:i:41 XA:Z:chr22,+50808429,41M3S,0; CR:Z:GAGGATGGTCATGGTA CY:Z:F,F,:,F,F::F,FFF CB:Z:GAGGATGGTCATCGTA-1 BC:Z:AGACTTTC QT:Z:FF,,FFF: TR:Z:CTGTCTCTTATACA TQ:Z:FFF:,,F:F,FFFF RG:Z:CtrlMix_C1:MissingLibrary:1:HJFT7DRX3:2 A00642:443:HJFT7DRX3:2:2205:4933:20055 181 chr1 9997 0 9997 0 ACACAGCTGCATTGAGTTTATGCTCTAGCCATTGTATATGGACTTTGGTC F,:FFFFFF,:FFFF,,FF,:F:F,,,,,,FFFFFF:,F,:F,::F,,,, AS:i:0 XS:i:0 CR:Z:GAGGATGGTCATGGTA CY:Z:F,F,:,F,F::F,FFF CB:Z:GAGGATGGTCATCGTA-1 BC:Z:AGACTTTC QT:Z:FF,,FFF: RG:Z:CtrlMix_C1:MissingLibrary:1:HJFT7DRX3:2 A00642:443:HJFT7DRX3:2:2105:6189:24706 69 chr1 9998 0 * 9998 0 GAGAAGAAATCCAAAAAAATAAATACATATTATACAATACCAGGCATGTCTCTTATAA F,FFFFF:,FFFF:,F,,:F,FF,FF::,FFF,F:,:F:F::::,,F,FF,,F:FF,, AS:i:0 XS:i:0 CR:Z:TCAGTTGGTATCAGTG CY:Z:,::F,,,F:,:,:,,: BC:Z:GATTGGTG QT:Z::F,,,,:F RG:Z:CtrlMix_B1:MissingLibrary:1:HJFT7DRX3:2 A00642:443:HJFT7DRX3:2:2105:6189:24706 137 chr1 9998 0 45M5S 9998 0 CGATAACCCTAACACTAACCCTAACCCTAACCCTAACCCTAACCCATGTC FF,,FFFFF:,FF,F:FFFF,,F,,FFFFFFFFFFFF,FFF,,F:FFFFF NM:i:1 MD:Z:13C31 AS:i:40XS:i:38 CR:Z:TCAGTTGGTATCAGTG CY:Z:,::F,,,F:,:,:,,: BC:Z:GATTGGTG QT:Z::F,,,,:F RG:Z:CtrlMix_B1:MissingLibrary:1:HJFT7DRX3:2

And the barcode file:

"AAACGAAAGGCGTCCT-1" "AAACGAAAGTACTCTG-1" "AAACGAAAGTTAGCAA-1" "AAACGAACAAGCAATA-1" "AAACGAACACATAAAG-1" "AAACGAACACGAACGA-1" "AAACGAACAGTGCTCG-1" "AAACGAAGTAAAGCTA-1" "AAACGAAGTACGGTTT-1" "AAACGAATCACACGTA-1"

Thank you for your continued help :)