search

KChen-lab / Monopogen

SNV calling from single cell sequencing
GNU General Public License v3.0
71 stars 17 forks source link

error in chr20.gl.vcf, broken VCF record, incorrect # of columns #5

Open chouellaine opened 1 year ago

chouellaine commented 1 year ago

Hello, I've tried running Monopogen with the below parameters and run into an error where there is a Incorrect number of FORMAT fields at chr20:2998993 .. DP, 0 vs 1 of the chr20.gl.vcf file that gets generated.

bcftools view chr20.gl.vcf.gz -r chr20:2998993 shows me that that position does not have the correct # of columns (see below) 

##FILTER=<ID=PASS,Description="All filters passed">
##samtoolsVersion=1.8+htslib-1.8
##samtoolsCommand=samtools mpileup -f /tmp/Monopogen/example/chr20_2Mb.hg38.fa -r chr20 -q 20 -Q 20 -t DP -d 10000000 -v /dbfs/FileStore/2023-h1-ancestry/data/monopogen_output/80_test/Bam/chr20.filter.bam
##reference=file:///tmp/Monopogen/example/chr20_2Mb.hg38.fa
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr10,length=133797422>
##contig=<ID=chr11,length=135086622>
##contig=<ID=chr12,length=133275309>
##contig=<ID=chr13,length=114364328>
##contig=<ID=chr14,length=107043718>
##contig=<ID=chr15,length=101991189>
##contig=<ID=chr16,length=90338345>
##contig=<ID=chr17,length=83257441>
##contig=<ID=chr18,length=80373285>
##contig=<ID=chr19,length=58617616>
##contig=<ID=chr2,length=242193529>
##contig=<ID=chr20,length=64444167>
##contig=<ID=chr21,length=46709983>
##contig=<ID=chr22,length=50818468>
##contig=<ID=chr3,length=198295559>
##contig=<ID=chr4,length=190214555>
##contig=<ID=chr5,length=181538259>
##contig=<ID=chr6,length=170805979>
##contig=<ID=chr7,length=159345973>
##contig=<ID=chr8,length=145138636>
##contig=<ID=chr9,length=138394717>
##contig=<ID=chrM,length=16569>
##contig=<ID=chrX,length=156040895>
##contig=<ID=chrY,length=57227415>
##contig=<ID=KI270728.1,length=1872759>
##contig=<ID=KI270727.1,length=448248>
##contig=<ID=KI270442.1,length=392061>
##contig=<ID=KI270729.1,length=280839>
##contig=<ID=GL000225.1,length=211173>
##contig=<ID=KI270743.1,length=210658>
##contig=<ID=GL000008.2,length=209709>
##contig=<ID=GL000009.2,length=201709>
##contig=<ID=KI270747.1,length=198735>
##contig=<ID=KI270722.1,length=194050>
##contig=<ID=GL000194.1,length=191469>
##contig=<ID=KI270742.1,length=186739>
##contig=<ID=GL000205.2,length=185591>
##contig=<ID=GL000195.1,length=182896>
##contig=<ID=KI270736.1,length=181920>
##contig=<ID=KI270733.1,length=179772>
##contig=<ID=GL000224.1,length=179693>
##contig=<ID=GL000219.1,length=179198>
##contig=<ID=KI270719.1,length=176845>
##contig=<ID=GL000216.2,length=176608>
##contig=<ID=KI270712.1,length=176043>
##contig=<ID=KI270706.1,length=175055>
##contig=<ID=KI270725.1,length=172810>
##contig=<ID=KI270744.1,length=168472>
##contig=<ID=KI270734.1,length=165050>
##contig=<ID=GL000213.1,length=164239>
##contig=<ID=GL000220.1,length=161802>
##contig=<ID=KI270715.1,length=161471>
##contig=<ID=GL000218.1,length=161147>
##contig=<ID=KI270749.1,length=158759>
##contig=<ID=KI270741.1,length=157432>
##contig=<ID=GL000221.1,length=155397>
##contig=<ID=KI270716.1,length=153799>
##contig=<ID=KI270731.1,length=150754>
##contig=<ID=KI270751.1,length=150742>
##contig=<ID=KI270750.1,length=148850>
##contig=<ID=KI270519.1,length=138126>
##contig=<ID=GL000214.1,length=137718>
##contig=<ID=KI270708.1,length=127682>
##contig=<ID=KI270730.1,length=112551>
##contig=<ID=KI270438.1,length=112505>
##contig=<ID=KI270737.1,length=103838>
##contig=<ID=KI270721.1,length=100316>
##contig=<ID=KI270738.1,length=99375>
##contig=<ID=KI270748.1,length=93321>
##contig=<ID=KI270435.1,length=92983>
##contig=<ID=GL000208.1,length=92689>
##contig=<ID=KI270538.1,length=91309>
##contig=<ID=KI270756.1,length=79590>
##contig=<ID=KI270739.1,length=73985>
##contig=<ID=KI270757.1,length=71251>
##contig=<ID=KI270709.1,length=66860>
##contig=<ID=KI270746.1,length=66486>
##contig=<ID=KI270753.1,length=62944>
##contig=<ID=KI270589.1,length=44474>
##contig=<ID=KI270726.1,length=43739>
##contig=<ID=KI270735.1,length=42811>
##contig=<ID=KI270711.1,length=42210>
##contig=<ID=KI270745.1,length=41891>
##contig=<ID=KI270714.1,length=41717>
##contig=<ID=KI270732.1,length=41543>
##contig=<ID=KI270713.1,length=40745>
##contig=<ID=KI270754.1,length=40191>
##contig=<ID=KI270710.1,length=40176>
##contig=<ID=KI270717.1,length=40062>
##contig=<ID=KI270724.1,length=39555>
##contig=<ID=KI270720.1,length=39050>
##contig=<ID=KI270723.1,length=38115>
##contig=<ID=KI270718.1,length=38054>
##contig=<ID=KI270317.1,length=37690>
##contig=<ID=KI270740.1,length=37240>
##contig=<ID=KI270755.1,length=36723>
##contig=<ID=KI270707.1,length=32032>
##contig=<ID=KI270579.1,length=31033>
##contig=<ID=KI270752.1,length=27745>
##contig=<ID=KI270512.1,length=22689>
##contig=<ID=KI270322.1,length=21476>
##contig=<ID=GL000226.1,length=15008>
##contig=<ID=KI270311.1,length=12399>
##contig=<ID=KI270366.1,length=8320>
##contig=<ID=KI270511.1,length=8127>
##contig=<ID=KI270448.1,length=7992>
##contig=<ID=KI270521.1,length=7642>
##contig=<ID=KI270581.1,length=7046>
##contig=<ID=KI270582.1,length=6504>
##contig=<ID=KI270515.1,length=6361>
##contig=<ID=KI270588.1,length=6158>
##contig=<ID=KI270591.1,length=5796>
##contig=<ID=KI270522.1,length=5674>
##contig=<ID=KI270507.1,length=5353>
##contig=<ID=KI270590.1,length=4685>
##contig=<ID=KI270584.1,length=4513>
##contig=<ID=KI270320.1,length=4416>
##contig=<ID=KI270382.1,length=4215>
##contig=<ID=KI270468.1,length=4055>
##contig=<ID=KI270467.1,length=3920>
##contig=<ID=KI270362.1,length=3530>
##contig=<ID=KI270517.1,length=3253>
##contig=<ID=KI270593.1,length=3041>
##contig=<ID=KI270528.1,length=2983>
##contig=<ID=KI270587.1,length=2969>
##contig=<ID=KI270364.1,length=2855>
##contig=<ID=KI270371.1,length=2805>
##contig=<ID=KI270333.1,length=2699>
##contig=<ID=KI270374.1,length=2656>
##contig=<ID=KI270411.1,length=2646>
##contig=<ID=KI270414.1,length=2489>
##contig=<ID=KI270510.1,length=2415>
##contig=<ID=KI270390.1,length=2387>
##contig=<ID=KI270375.1,length=2378>
##contig=<ID=KI270420.1,length=2321>
##contig=<ID=KI270509.1,length=2318>
##contig=<ID=KI270315.1,length=2276>
##contig=<ID=KI270302.1,length=2274>
##contig=<ID=KI270518.1,length=2186>
##contig=<ID=KI270530.1,length=2168>
##contig=<ID=KI270304.1,length=2165>
##contig=<ID=KI270418.1,length=2145>
##contig=<ID=KI270424.1,length=2140>
##contig=<ID=KI270417.1,length=2043>
##contig=<ID=KI270508.1,length=1951>
##contig=<ID=KI270303.1,length=1942>
##contig=<ID=KI270381.1,length=1930>
##contig=<ID=KI270529.1,length=1899>
##contig=<ID=KI270425.1,length=1884>
##contig=<ID=KI270396.1,length=1880>
##contig=<ID=KI270363.1,length=1803>
##contig=<ID=KI270386.1,length=1788>
##contig=<ID=KI270465.1,length=1774>
##contig=<ID=KI270383.1,length=1750>
##contig=<ID=KI270384.1,length=1658>
##contig=<ID=KI270330.1,length=1652>
##contig=<ID=KI270372.1,length=1650>
##contig=<ID=KI270548.1,length=1599>
##contig=<ID=KI270580.1,length=1553>
##contig=<ID=KI270387.1,length=1537>
##contig=<ID=KI270391.1,length=1484>
##contig=<ID=KI270305.1,length=1472>
##contig=<ID=KI270373.1,length=1451>
##contig=<ID=KI270422.1,length=1445>
##contig=<ID=KI270316.1,length=1444>
##contig=<ID=KI270340.1,length=1428>
##contig=<ID=KI270338.1,length=1428>
##contig=<ID=KI270583.1,length=1400>
##contig=<ID=KI270334.1,length=1368>
##contig=<ID=KI270429.1,length=1361>
##contig=<ID=KI270393.1,length=1308>
##contig=<ID=KI270516.1,length=1300>
##contig=<ID=KI270389.1,length=1298>
##contig=<ID=KI270466.1,length=1233>
##contig=<ID=KI270388.1,length=1216>
##contig=<ID=KI270544.1,length=1202>
##contig=<ID=KI270310.1,length=1201>
##contig=<ID=KI270412.1,length=1179>
##contig=<ID=KI270395.1,length=1143>
##contig=<ID=KI270376.1,length=1136>
##contig=<ID=KI270337.1,length=1121>
##contig=<ID=KI270335.1,length=1048>
##contig=<ID=KI270378.1,length=1048>
##contig=<ID=KI270379.1,length=1045>
##contig=<ID=KI270329.1,length=1040>
##contig=<ID=KI270419.1,length=1029>
##contig=<ID=KI270336.1,length=1026>
##contig=<ID=KI270312.1,length=998>
##contig=<ID=KI270539.1,length=993>
##contig=<ID=KI270385.1,length=990>
##contig=<ID=KI270423.1,length=981>
##contig=<ID=KI270392.1,length=971>
##contig=<ID=KI270394.1,length=970>
##ALT=<ID=*,Description="Represents allele(s) other than observed.">
##INFO=<ID=IDV,Number=1,Type=Integer,Description="Maximum number of reads supporting an indel">
##INFO=<ID=IMF,Number=1,Type=Float,Description="Maximum fraction of reads supporting an indel">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias for filtering splice-site artefacts in RNA-seq data (bigger is better)",Version="3">
##INFO=<ID=RPB,Number=1,Type=Float,Description="Mann-Whitney U test of Read Position Bias (bigger is better)">
##INFO=<ID=MQB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality Bias (bigger is better)">
##INFO=<ID=BQB,Number=1,Type=Float,Description="Mann-Whitney U test of Base Quality Bias (bigger is better)">
##INFO=<ID=MQSB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality vs Strand Bias (bigger is better)">
##INFO=<ID=SGB,Number=1,Type=Float,Description="Segregation based metric.">
##INFO=<ID=MQ0F,Number=1,Type=Float,Description="Fraction of MQ0 reads (smaller is better)">
##INFO=<ID=I16,Number=16,Type=Float,Description="Auxiliary tag used for calling, see description of bcf_callret1_t in bam2bcf.h">
##INFO=<ID=QS,Number=R,Type=Float,Description="Auxiliary tag used for calling">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Number of high-quality bases">
##bcftools_viewVersion=1.8+htslib-1.8
##bcftools_viewCommand=view; Date=Wed Mar 15 02:47:19 2023
##bcftools_normVersion=1.8+htslib-1.8
##bcftools_normCommand=norm -m-both -f /tmp/Monopogen/example/chr20_2Mb.hg38.fa; Date=Wed Mar 15 02:47:19 2023
##bcftools_viewCommand=view -r chr20:2998993 /dbfs/FileStore/2023-h1-ancestry/data/monopogen_output/80_test/SCvarCall/chr20.gl.vcf.gz; Date=Wed Mar 15 03:18:53 2023
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  10x3_Lobe_19_D006_NeuN_4
[E::bcf_write] Broken VCF record, the number of columns at chr20:2998993 does not match the number of samples (0 vs 1)```

**Here's the output log**

[2023-03-15 20:01:59,456] INFO Monopogen.py Preparing varint calling pipeline... [2023-03-15 20:01:59,457] INFO Monopogen.py Parameters in effect: [2023-03-15 20:01:59,457] INFO Monopogen.py --subcommand = [SCvarCall] [2023-03-15 20:01:59,457] INFO Monopogen.py --bamFile = [/dbfs/FileStore/2023-h1-ancestry/data/lattice/80.bam] [2023-03-15 20:01:59,457] INFO Monopogen.py --step = [germline] [2023-03-15 20:01:59,457] INFO Monopogen.py --mode = [single] [2023-03-15 20:01:59,457] INFO Monopogen.py --chr = [chr20] [2023-03-15 20:01:59,458] INFO Monopogen.py --out = [/dbfs/FileStore/2023-h1-ancestry/data/monopogen_output/80_test] [2023-03-15 20:01:59,458] INFO Monopogen.py --reference = [/tmp/Monopogen/example/chr20_2Mb.hg38.fa] [2023-03-15 20:01:59,458] INFO Monopogen.py --imputation_panel = [/tmp/Monopogen/example/CCDG_14151_B01_GRM_WGS_2020-08-05_chr20.filtered.shapeit2-duohmm-phased.vcf.gz] [2023-03-15 20:01:59,458] INFO Monopogen.py --depth_filter = [200] [2023-03-15 20:01:59,458] INFO Monopogen.py --depth_filter_monovar = [50] [2023-03-15 20:01:59,458] INFO Monopogen.py --alt_ratio = [0.1] [2023-03-15 20:01:59,458] INFO Monopogen.py --wildCluster = [2] [2023-03-15 20:01:59,458] INFO Monopogen.py --mutationCluster = [1] [2023-03-15 20:01:59,458] INFO Monopogen.py --max_mismatch = [3] [2023-03-15 20:01:59,458] INFO Monopogen.py --max_softClipped = [5] [2023-03-15 20:01:59,459] INFO Monopogen.py --app_path = [/tmp/Monopogen/apps] [2023-03-15 20:01:59,459] INFO Monopogen.py --cell_cluster = [/tmp/Monopogen/example/cell_cluster.csv] [2023-03-15 20:01:59,459] INFO Monopogen.py --logfile = [/dbfs/FileStore/2023-h1-ancestry/data/monopogen_output/80_test_chr20.log] [2023-03-15 20:01:59,459] INFO Monopogen.py --samtools = [/tmp/Monopogen/apps/samtools] [2023-03-15 20:01:59,459] INFO Monopogen.py --bcftools = [/tmp/Monopogen/apps/bcftools] [2023-03-15 20:01:59,459] INFO Monopogen.py --bgzip = [/tmp/Monopogen/apps/bgzip] [2023-03-15 20:01:59,459] INFO Monopogen.py --java = [java] [2023-03-15 20:01:59,459] INFO Monopogen.py --beagle = [/tmp/Monopogen/apps/beagle.27Jul16.86a.jar] [2023-03-15 20:01:59,459] INFO Monopogen.py Checking existence of essenstial resource files... [2023-03-15 20:01:59,460] INFO Monopogen.py Checking dependencies... [2023-03-15 20:01:59,566] INFO Monopogen.py Filtering bam files... [2023-03-15 20:03:06,093] INFO Monopogen.py Performing Variant Calling... [W::hts_idx_load2] The index file is older than the data file: /dbfs/FileStore/2023-h1-ancestry/data/monopogen_output/80_test/Bam/chr20.filter.bam.bai [mpileup] 1 samples in 1 input files (mpileup) Max depth is above 1M. Potential memory hog! Reference allele mismatch at chr20:3000001 .. REF_SEQ:'A' vs VCF:'N' Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid Incorrect number of FORMAT fields at chr20:2998993 .. DP, 0 vs 1 Exception in thread "main" java.lang.IllegalArgumentException: VCF header line has 10 fields, but data line has 8 fields File source:File source: /dbfs/FileStore/2023-h1-ancestry/data/monopogen_output/80_test/SCvarCall/chr20.gl.vcf.gz [chr20, 2998993, ., T, <*>, 0, ., DP=9;I16=6,3,0,0,333] at vcf.VcfRecord.fieldCountError(VcfRecord.java:221) at vcf.VcfRecord.delimiters(VcfRecord.java:203) at vcf.VcfRecord.(VcfRecord.java:87) at vcf.VcfRecord.fromGTGL(VcfRecord.java:193) at vcf.VcfIt.lambda$static$5(VcfIt.java:76) at vcf.VcfIt.lambda$new$8(VcfIt.java:192) at java.util.stream.ReferencePipeline$3$1.accept(ReferencePipeline.java:193) at java.util.Spliterators$ArraySpliterator.forEachRemaining(Spliterators.java:948) at java.util.stream.AbstractPipeline.copyInto(AbstractPipeline.java:482) at java.util.stream.AbstractPipeline.wrapAndCopyInto(AbstractPipeline.java:472) at java.util.stream.ReduceOps$ReduceTask.doLeaf(ReduceOps.java:747) at java.util.stream.ReduceOps$ReduceTask.doLeaf(ReduceOps.java:721) at java.util.stream.AbstractTask.compute(AbstractTask.java:327) at java.util.concurrent.CountedCompleter.exec(CountedCompleter.java:731) at java.util.concurrent.ForkJoinTask.doExec(ForkJoinTask.java:289) at java.util.concurrent.ForkJoinPool$WorkQueue.pollAndExecCC(ForkJoinPool.java:1190) at java.util.concurrent.ForkJoinPool.helpComplete(ForkJoinPool.java:1879) at java.util.concurrent.ForkJoinPool.externalHelpComplete(ForkJoinPool.java:2467) at java.util.concurrent.ForkJoinTask.externalAwaitDone(ForkJoinTask.java:324) at java.util.concurrent.ForkJoinTask.doInvoke(ForkJoinTask.java:405) at java.util.concurrent.ForkJoinTask.invoke(ForkJoinTask.java:734) at java.util.stream.ReduceOps$ReduceOp.evaluateParallel(ReduceOps.java:714) at java.util.stream.AbstractPipeline.evaluate(AbstractPipeline.java:233) at java.util.stream.ReferencePipeline.collect(ReferencePipeline.java:566) at vcf.VcfIt.fillEmissionBuffer(VcfIt.java:307) at vcf.VcfIt.next(VcfIt.java:363) at vcf.VcfIt.next(VcfIt.java:52) at vcf.IntervalVcfIt.readNextRecord(IntervalVcfIt.java:110) at vcf.IntervalVcfIt.next(IntervalVcfIt.java:92) at vcf.IntervalVcfIt.next(IntervalVcfIt.java:36) at main.Main.restrictToVcfMarkers(Main.java:343) at main.Main.allData(Main.java:313) at main.Main.main(Main.java:111) gzip: /dbfs/FileStore/2023-h1-ancestry/data/monopogen_output/80_test/SCvarCall/chr20.gp.vcf.gz: No such file or directory /dbfs/FileStore/2023-h1-ancestry/data/monopogen_output/80_test/SCvarCall/chr20.gp.vcf.gz: No such file or directory [2023-03-15 20:26:35,885] INFO Monopogen.py Generating variant statistical information ... [E::hts_open_format] Failed to open file "/dbfs/FileStore/2023-h1-ancestry/data/monopogen_output/80_test/SCvarCall/chr20.gp.vcf.gz" : No such file or directory Traceback (most recent call last): File "/tmp/Monopogen/src/Monopogen.py", line 519, in main() File "/tmp/Monopogen/src/Monopogen.py", line 514, in main args.func(args) File "/tmp/Monopogen/src/Monopogen.py", line 398, in SCvarCall getDPinfo(args) File "/tmp/Monopogen/src/Monopogen.py", line 152, in getDPinfo gp_vcf_in = VariantFile(out + "/SCvarCall/" + args.chr + ".gp.vcf.gz") File "pysam/libcbcf.pyx", line 4119, in pysam.libcbcf.VariantFile.init File "pysam/libcbcf.pyx", line 4344, in pysam.libcbcf.VariantFile.open FileNotFoundError: [Errno 2] could not open variant file b'/dbfs/FileStore/2023-h1-ancestry/data/monopogen_output/80_test/SCvarCall/chr20.gp.vcf.gz': No such file or directory