Merge pull request #178 from timothymillar/176-correct-AN

Version 0.9.3

.github/workflows/python-package.yml

@@ -15,7 +15,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["3.8", "3.9", "3.10"]
+        python-version: ["3.8", "3.9", "3.10", "3.11"]
 
     steps:
     - uses: actions/checkout@v2
@@ -34,6 +34,13 @@ jobs:
       run: |
         python setup.py sdist
         pip install dist/mchap-*.tar.gz
+    - name: Test with pytest (bounds checked)
+      env:
+        NUMBA_BOUNDSCHECK: 1
+      run: |
+        pytest -vv
     - name: Test with pytest
+      env:
+        NUMBA_BOUNDSCHECK: 0
       run: |
         pytest -vv

CHANGELOG.md

@@ -3,6 +3,12 @@
 ## Unreleased 
 
 
+## Beta v0.9.3
+
+Bug Fixes:
+- Correct usage of the `AN` field #176
+- Improved error messages for io #163 #177
+
 ## Beta v0.9.2
 
 Bug Fixes:

Dockerfile

@@ -9,7 +9,7 @@ RUN apt-get update \
 
 RUN git clone  https://github.com/PlantandFoodResearch/MCHap.git   \
     && cd MCHap \
-    && git checkout v0.9.2 \
+    && git checkout v0.9.3 \
     && pip install -r requirements.txt \
     && python3 setup.py sdist \
     && python3 -m pip install dist/mchap-*tar.gz

docs/assemble.rst

@@ -3,7 +3,7 @@ MCHap assemble
 
 De novo assembly of micro-haplotypes.
 
-*(Last updated for MCHap version 0.9.2)*
+*(Last updated for MCHap version 0.9.3)*
 
 Background
 ----------

docs/call.rst

@@ -3,7 +3,7 @@ MCHap call
 
 Calling genotypes from known haplotypes.
 
-*(Last updated for MCHap version 0.9.2)*
+*(Last updated for MCHap version 0.9.3)*
 
 Background
 ----------
@@ -341,7 +341,7 @@ at least ``1000`` steps should be kept to calculate posterior probabilities.
 Excluding input haplotypes
 ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-The speed of each MCMC step in ``mchap assemble`` is largely dependant on the
+The speed of each MCMC step in ``mchap call`` is largely dependant on the
 ploidy of an individual and the number of unique haplotypes in the input VCF file.
 Therefore, the speed of analysis can be improved by minimizing unnecessary
 haplotypes from the input VCF file.

mchap/application/baseclass.py

@@ -305,7 +305,7 @@ def sumarise_vcf_record(self, data):
                     allele_counts[a] += 1
         data.infodata["AC"] = allele_counts[1:]  # skip ref count
         # total number of alleles in called genotypes
-        data.infodata["AN"] = np.sum(allele_counts > 0)
+        data.infodata["AN"] = np.sum(allele_counts)
         # number of called samples
         data.infodata["NS"] = sum(
             np.any(a >= 0) for a in data.sampledata["alleles"].values()

mchap/io/bam.py

@@ -93,7 +93,7 @@ def extract_read_variants(
 
     """
 
-    assert id in {"ID", "SM"}
+    assert id in ID_TAGS
 
     # a single sample name may be given as a string
     if isinstance(samples, str):
@@ -182,7 +182,19 @@ def extract_read_variants(
                         idx = positions[ref_pos]
 
                         # references allele in bam should match reference allele in locus
-                        assert locus.alleles[idx][0].upper() == ref_char.upper()
+                        if locus.alleles[idx][0].upper() != ref_char.upper():
+                            path = alignment_file.filename.decode()
+                            locus_ref_char = locus.alleles[idx][0]
+                            locus_contig = locus.contig
+                            locus_name = locus.name
+                            vcf_pos = ref_pos + 1
+                            if locus_name:
+                                loc = f"'{locus_contig}:{vcf_pos}' in target '{locus_name}'"
+                            else:
+                                loc = f"'{locus_contig}:{vcf_pos}'"
+                            raise ValueError(
+                                f"Reference allele of variant '{locus_ref_char}' does not match alignment reference allele '{ref_char}' at position {loc} in '{path}'"
+                            )
 
                         char = read.seq[read_pos]
                         qual = util.qual_of_char(read.qual[read_pos])

mchap/io/loci.py

@@ -66,28 +66,82 @@ def set(self, **kwargs):
         data.update(kwargs)
         return type(self)(**data)
 
+    def validate_reference_alleles(self):
+        sequence_chars = list(self.sequence)
+        ref_alleles = (tup[0] for tup in self.alleles)
+        for pos, char in zip(self.positions, ref_alleles):
+            idx = pos - self.start
+            seq_char = sequence_chars[idx]
+            if seq_char != char:
+                contig = self.contig
+                name = self.name
+                vcf_pos = pos + 1
+                if name:
+                    loc = f"'{contig}:{vcf_pos}' in target '{name}'"
+                else:
+                    loc = f"'{contig}:{vcf_pos}'"
+                raise ValueError(
+                    f"Reference allele of variant '{char}' does not match reference sequence '{seq_char}' at {loc}"
+                )
+
     def set_sequence(self, fasta):
         with pysam.FastaFile(fasta) as f:
-            return _set_locus_sequence(self, f)
+            sequence = f.fetch(self.contig, self.start, self.stop).upper()
+            locus = self.set(sequence=sequence)
+            if locus.variants:
+                locus.validate_reference_alleles()
+            return locus
 
     def set_variants(self, vcf):
         with pysam.VariantFile(vcf) as f:
-            return _set_locus_variants(self, f)
+            variants = []
+            positions = set()
+            try:
+                records = f.fetch(self.contig, self.start, self.stop)
+            except ValueError as e:
+                if "fetch requires an index" in e.args:
+                    raise ValueError(
+                        f"Could not fetch variants from file '{vcf}', the file is not indexed"
+                    ) from e
+                else:
+                    raise e
+            for var in records:
+                alleles = (var.ref,) + var.alts
+
+                if (var.stop - var.start == 1) and all(len(a) == 1 for a in alleles):
+                    # is a SNP
+                    snp = SNP(
+                        contig=var.contig,
+                        start=var.start,
+                        stop=var.stop,
+                        name=var.id if var.id else ".",
+                        alleles=alleles,
+                    )
+                    if snp.start in positions:
+                        # attempt to merge duplicates
+                        variants = [
+                            _merge_snps(s, snp) if s.start == snp.start else s
+                            for s in variants
+                        ]
+                    else:
+                        variants.append(snp)
+                        positions.add(snp.start)
+                else:
+                    pass
+
+            variants = tuple(variants)
+            locus = self.set(variants=variants)
+            if locus.sequence:
+                locus.validate_reference_alleles()
+            return locus
 
     def _template_sequence(self):
         chars = list(self.sequence)
-        ref_alleles = (tup[0] for tup in self.alleles)
-        for pos, string in zip(self.positions, ref_alleles):
+        for pos in self.positions:
             idx = pos - self.start
-            for offset, char in enumerate(string):
-                if chars[idx + offset] != char:
-                    message = (
-                        "Reference allele does not match sequence at position {}:{}"
-                    )
-                    raise ValueError(message.format(self.contig, pos + offset))
 
-                # remove chars
-                chars[idx + offset] = ""
+            # remove chars
+            chars[idx] = ""
 
             # add template position
             chars[idx] = "{}"
@@ -309,12 +363,6 @@ def read_bed4(bed, region=None):
                     yield _parse_bed4_line(line.decode())
 
 
-def _set_locus_sequence(locus, fasta_file):
-    sequence = fasta_file.fetch(locus.contig, locus.start, locus.stop).upper()
-    locus = locus.set(sequence=sequence)
-    return locus
-
-
 def _merge_snps(x, y):
     match = [
         x.contig == y.contig,
@@ -340,35 +388,3 @@ def _merge_snps(x, y):
         name=x.name,
         alleles=alleles,
     )
-
-
-def _set_locus_variants(locus, variant_file):
-    variants = []
-    positions = set()
-
-    for var in variant_file.fetch(locus.contig, locus.start, locus.stop):
-        alleles = (var.ref,) + var.alts
-
-        if (var.stop - var.start == 1) and all(len(a) == 1 for a in alleles):
-            # is a SNP
-            snp = SNP(
-                contig=var.contig,
-                start=var.start,
-                stop=var.stop,
-                name=var.id if var.id else ".",
-                alleles=alleles,
-            )
-            if snp.start in positions:
-                # attempt to merge duplicates
-                variants = [
-                    _merge_snps(s, snp) if s.start == snp.start else s for s in variants
-                ]
-            else:
-                variants.append(snp)
-                positions.add(snp.start)
-        else:
-            pass
-
-    variants = tuple(variants)
-    locus = locus.set(variants=variants)
-    return locus

mchap/tests/test_io/data/mock.input.frequencies.vcf

@@ -1,6 +1,6 @@
 ##fileformat=VCFv4.3
 ##fileDate=20220406
-##source=mchap v0.9.2
+##source=mchap v0.9.3
 ##phasing=None
 ##commandline="mockup"
 ##randomseed=11

mchap/tests/test_io/data/simple.output.assemble.vcf

@@ -1,6 +1,6 @@
 ##fileformat=VCFv4.3
 ##fileDate=20230710
-##source=mchap v0.9.2
+##source=mchap v0.9.3
 ##phasing=None
 ##commandline="mchap assemble --bam simple.sample1.bam simple.sample2.bam simple.sample3.bam --ploidy 4 --targets simple.bed.gz --variants simple.vcf.gz --reference simple.fasta --mcmc-steps 500 --mcmc-burn 100 --mcmc-seed 11"
 ##randomseed=11
@@ -31,7 +31,7 @@
 ##FORMAT=<ID=PHPM,Number=1,Type=Float,Description="Phenotype posterior mode probability">
 ##FORMAT=<ID=MCI,Number=1,Type=Integer,Description="Replicate Markov-chain incongruence, 0 = none, 1 = incongruence, 2 = putative CNV">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	SAMPLE1	SAMPLE2	SAMPLE3
-CHR1	6	CHR1_05_25	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAAGAAAAAATAA,ACAAAAAAAAGAAAAAACAA	.	PASS	AN=3;AC=3,2;NS=3;DP=39;RCOUNT=60;END=25;NVAR=3;SNVPOS=2,11,18	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/1/2:7:8:13:20:40:0:0:0.785:0.835:0	0/0/1/1:4:5:13:20:40:0:0:0.636:0.678:0	0/0/0/2:4:5:13:20:40:0:0:0.619:0.698:0
-CHR1	31	CHR1_30_50	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=1;AC=.;NS=3;DP=.;RCOUNT=72;END=50;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:24:0:0:.:1:1:0	0/0/0/0:60:60:.:24:0:0:.:1:1:0	0/0/0/0:60:60:.:24:0:0:.:1:1:0
-CHR2	11	CHR2_10_30	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAGAAAAAAAAAA,AAAAAAAAATAAAAAAAAAA,AAAATAAAAGAAAAAAAAAA	.	PASS	AN=4;AC=3,2,1;NS=3;DP=42;RCOUNT=72;END=30;NVAR=2;SNVPOS=5,10	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/2:6:8:14:24:28:0:0:0.72:0.84:0	0/0/1/2:7:11:14:24:28:0:0:0.815:0.929:0	0/1/1/3:7:13:14:24:28:0:0:0.801:0.949:0
-CHR3	21	CHR3_20_40	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=1;AC=.;NS=3;DP=.;RCOUNT=0;END=40;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:0:0:0:.:1:1:0	0/0/0/0:60:60:.:0:0:0:.:1:1:0	0/0/0/0:60:60:.:0:0:0:.:1:1:0
+CHR1	6	CHR1_05_25	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAAGAAAAAATAA,ACAAAAAAAAGAAAAAACAA	.	PASS	AN=12;AC=3,2;NS=3;DP=39;RCOUNT=60;END=25;NVAR=3;SNVPOS=2,11,18	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/1/2:7:8:13:20:40:0:0:0.785:0.835:0	0/0/1/1:4:5:13:20:40:0:0:0.636:0.678:0	0/0/0/2:4:5:13:20:40:0:0:0.619:0.698:0
+CHR1	31	CHR1_30_50	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=12;AC=.;NS=3;DP=.;RCOUNT=72;END=50;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:24:0:0:.:1:1:0	0/0/0/0:60:60:.:24:0:0:.:1:1:0	0/0/0/0:60:60:.:24:0:0:.:1:1:0
+CHR2	11	CHR2_10_30	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAGAAAAAAAAAA,AAAAAAAAATAAAAAAAAAA,AAAATAAAAGAAAAAAAAAA	.	PASS	AN=12;AC=3,2,1;NS=3;DP=42;RCOUNT=72;END=30;NVAR=2;SNVPOS=5,10	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/2:6:8:14:24:28:0:0:0.72:0.84:0	0/0/1/2:7:11:14:24:28:0:0:0.815:0.929:0	0/1/1/3:7:13:14:24:28:0:0:0.801:0.949:0
+CHR3	21	CHR3_20_40	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=12;AC=.;NS=3;DP=.;RCOUNT=0;END=40;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:0:0:0:.:1:1:0	0/0/0/0:60:60:.:0:0:0:.:1:1:0	0/0/0/0:60:60:.:0:0:0:.:1:1:0

mchap/tests/test_io/data/simple.output.basis.minad2.vcf

@@ -1,6 +1,6 @@
 ##fileformat=VCFv4.3
 ##fileDate=20230710
-##source=mchap v0.9.2
+##source=mchap v0.9.3
 ##commandline="mchap find-snvs --targets simple.bed --reference simple.fasta --bam simple.sample1.bam simple.sample2.bam simple.sample3.bam --ind-mad 2"
 ##reference=file:simple.fasta
 ##contig=<ID=CHR1,length=60>

mchap/tests/test_io/data/simple.output.basis.minaf0.3.vcf

@@ -1,6 +1,6 @@
 ##fileformat=VCFv4.3
 ##fileDate=20230710
-##source=mchap v0.9.2
+##source=mchap v0.9.3
 ##commandline="mchap find-snvs --targets simple.bed --reference simple.fasta --bam simple.sample1.bam simple.sample2.bam simple.sample3.bam --ind-maf 0.3"
 ##reference=file:simple.fasta
 ##contig=<ID=CHR1,length=60>

mchap/tests/test_io/data/simple.output.basis.minaf0.minad0.vcf

@@ -1,6 +1,6 @@
 ##fileformat=VCFv4.3
 ##fileDate=20230710
-##source=mchap v0.9.2
+##source=mchap v0.9.3
 ##commandline="mchap find-snvs --targets simple.bed --reference simple.fasta --bam simple.sample1.bam simple.sample2.bam simple.sample3.bam --ind-maf 0.0 --ind-mad 0"
 ##reference=file:simple.fasta
 ##contig=<ID=CHR1,length=60>

mchap/tests/test_io/data/simple.output.basis.mixed_depth.mad10.vcf

@@ -1,6 +1,6 @@
 ##fileformat=VCFv4.3
 ##fileDate=20230710
-##source=mchap v0.9.2
+##source=mchap v0.9.3
 ##commandline="mchap find-snvs --targets simple.bed --reference simple.fasta --bam simple.sample1.bam simple.sample2.deep.bam simple.sample3.bam --ind-maf 0 --ind-mad 0 --mad 10"
 ##reference=file:simple.fasta
 ##contig=<ID=CHR1,length=60>

mchap/tests/test_io/data/simple.output.basis.mixed_depth.maf0.1.vcf

@@ -1,6 +1,6 @@
 ##fileformat=VCFv4.3
 ##fileDate=20230710
-##source=mchap v0.9.2
+##source=mchap v0.9.3
 ##commandline="mchap find-snvs --targets simple.bed --reference simple.fasta --bam simple.sample1.bam simple.sample2.deep.bam simple.sample3.bam --ind-maf 0 --ind-mad 0 --maf 0.1"
 ##reference=file:simple.fasta
 ##contig=<ID=CHR1,length=60>

mchap/tests/test_io/data/simple.output.basis.mixed_depth.vcf

@@ -1,6 +1,6 @@
 ##fileformat=VCFv4.3
 ##fileDate=20230710
-##source=mchap v0.9.2
+##source=mchap v0.9.3
 ##commandline="mchap find-snvs --targets simple.bed --reference simple.fasta --bam simple.sample1.bam simple.sample2.deep.bam simple.sample3.bam"
 ##reference=file:simple.fasta
 ##contig=<ID=CHR1,length=60>

mchap/tests/test_io/data/simple.output.basis.vcf

@@ -1,6 +1,6 @@
 ##fileformat=VCFv4.3
 ##fileDate=20230710
-##source=mchap v0.9.2
+##source=mchap v0.9.3
 ##commandline="mchap find-snvs --targets simple.bed --reference simple.fasta --bam simple.sample1.bam simple.sample2.bam simple.sample3.bam"
 ##reference=file:simple.fasta
 ##contig=<ID=CHR1,length=60>

mchap/tests/test_io/data/simple.output.call-exact.vcf

@@ -1,6 +1,6 @@
 ##fileformat=VCFv4.3
 ##fileDate=20230710
-##source=mchap v0.9.2
+##source=mchap v0.9.3
 ##phasing=None
 ##commandline="mchap call-exact --bam simple.sample1.bam simple.sample2.bam simple.sample3.bam --ploidy 4 --haplotypes simple.output.assemble.vcf"
 ##randomseed=None
@@ -31,7 +31,7 @@
 ##FORMAT=<ID=PHPM,Number=1,Type=Float,Description="Phenotype posterior mode probability">
 ##FORMAT=<ID=MCI,Number=1,Type=Integer,Description="Replicate Markov-chain incongruence, 0 = none, 1 = incongruence, 2 = putative CNV">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	SAMPLE1	SAMPLE2	SAMPLE3
-CHR1	6	CHR1_05_25	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAAGAAAAAATAA,ACAAAAAAAAGAAAAAACAA	.	PASS	AN=3;AC=3,2;NS=3;DP=39;RCOUNT=60;END=25;NVAR=3;SNVPOS=2,11,18	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/1/2:12:60:13:20:40:0:0:0.941:1:.	0/0/1/1:11:24:13:20:40:0:0:0.926:0.996:.	0/0/0/2:10:22:13:20:40:0:0:0.896:0.994:.
-CHR1	31	CHR1_30_50	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=1;AC=.;NS=3;DP=.;RCOUNT=72;END=50;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:24:0:0:.:1:1:.	0/0/0/0:60:60:.:24:0:0:.:1:1:.	0/0/0/0:60:60:.:24:0:0:.:1:1:.
-CHR2	11	CHR2_10_30	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAGAAAAAAAAAA,AAAAAAAAATAAAAAAAAAA,AAAATAAAAGAAAAAAAAAA	.	PASS	AN=4;AC=3,2,1;NS=3;DP=42;RCOUNT=72;END=30;NVAR=2;SNVPOS=5,10	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/2:7:14:14:24:28:0:0:0.812:0.962:.	0/0/1/2:9:19:14:24:28:0:0:0.876:0.986:.	0/1/1/3:8:21:14:24:28:0:0:0.827:0.992:.
-CHR3	21	CHR3_20_40	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=1;AC=.;NS=3;DP=.;RCOUNT=0;END=40;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:0:0:0:.:1:1:.	0/0/0/0:60:60:.:0:0:0:.:1:1:.	0/0/0/0:60:60:.:0:0:0:.:1:1:.
+CHR1	6	CHR1_05_25	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAAGAAAAAATAA,ACAAAAAAAAGAAAAAACAA	.	PASS	AN=12;AC=3,2;NS=3;DP=39;RCOUNT=60;END=25;NVAR=3;SNVPOS=2,11,18	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/1/2:12:60:13:20:40:0:0:0.941:1:.	0/0/1/1:11:24:13:20:40:0:0:0.926:0.996:.	0/0/0/2:10:22:13:20:40:0:0:0.896:0.994:.
+CHR1	31	CHR1_30_50	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=12;AC=.;NS=3;DP=.;RCOUNT=72;END=50;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:24:0:0:.:1:1:.	0/0/0/0:60:60:.:24:0:0:.:1:1:.	0/0/0/0:60:60:.:24:0:0:.:1:1:.
+CHR2	11	CHR2_10_30	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAGAAAAAAAAAA,AAAAAAAAATAAAAAAAAAA,AAAATAAAAGAAAAAAAAAA	.	PASS	AN=12;AC=3,2,1;NS=3;DP=42;RCOUNT=72;END=30;NVAR=2;SNVPOS=5,10	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/2:7:14:14:24:28:0:0:0.812:0.962:.	0/0/1/2:9:19:14:24:28:0:0:0.876:0.986:.	0/1/1/3:8:21:14:24:28:0:0:0.827:0.992:.
+CHR3	21	CHR3_20_40	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=12;AC=.;NS=3;DP=.;RCOUNT=0;END=40;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:0:0:0:.:1:1:.	0/0/0/0:60:60:.:0:0:0:.:1:1:.	0/0/0/0:60:60:.:0:0:0:.:1:1:.

mchap/tests/test_io/data/simple.output.call.vcf

@@ -1,6 +1,6 @@
 ##fileformat=VCFv4.3
 ##fileDate=20230710
-##source=mchap v0.9.2
+##source=mchap v0.9.3
 ##phasing=None
 ##commandline="mchap call --bam simple.sample1.bam simple.sample2.bam simple.sample3.bam --ploidy 4 --haplotypes simple.output.assemble.vcf --mcmc-steps 500 --mcmc-burn 100 --mcmc-seed 11"
 ##randomseed=11
@@ -31,7 +31,7 @@
 ##FORMAT=<ID=PHPM,Number=1,Type=Float,Description="Phenotype posterior mode probability">
 ##FORMAT=<ID=MCI,Number=1,Type=Integer,Description="Replicate Markov-chain incongruence, 0 = none, 1 = incongruence, 2 = putative CNV">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	SAMPLE1	SAMPLE2	SAMPLE3
-CHR1	6	CHR1_05_25	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAAGAAAAAATAA,ACAAAAAAAAGAAAAAACAA	.	PASS	AN=3;AC=3,2;NS=3;DP=39;RCOUNT=60;END=25;NVAR=3;SNVPOS=2,11,18	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/1/2:13:60:13:20:40:0:0:0.952:1:0	0/0/1/1:11:26:13:20:40:0:0:0.928:0.998:0	0/0/0/2:10:19:13:20:40:0:0:0.894:0.989:0
-CHR1	31	CHR1_30_50	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=1;AC=.;NS=3;DP=.;RCOUNT=72;END=50;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:24:0:0:.:1:1:0	0/0/0/0:60:60:.:24:0:0:.:1:1:0	0/0/0/0:60:60:.:24:0:0:.:1:1:0
-CHR2	11	CHR2_10_30	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAGAAAAAAAAAA,AAAAAAAAATAAAAAAAAAA,AAAATAAAAGAAAAAAAAAA	.	PASS	AN=4;AC=3,2,1;NS=3;DP=42;RCOUNT=72;END=30;NVAR=2;SNVPOS=5,10	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/2:7:15:14:24:28:0:0:0.821:0.969:0	0/0/1/2:10:21:14:24:28:0:0:0.9:0.991:0	0/1/1/3:7:21:14:24:28:0:0:0.806:0.992:0
-CHR3	21	CHR3_20_40	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=1;AC=.;NS=3;DP=.;RCOUNT=0;END=40;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:0:0:0:.:1:1:0	0/0/0/0:60:60:.:0:0:0:.:1:1:0	0/0/0/0:60:60:.:0:0:0:.:1:1:0
+CHR1	6	CHR1_05_25	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAAGAAAAAATAA,ACAAAAAAAAGAAAAAACAA	.	PASS	AN=12;AC=3,2;NS=3;DP=39;RCOUNT=60;END=25;NVAR=3;SNVPOS=2,11,18	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/1/2:13:60:13:20:40:0:0:0.952:1:0	0/0/1/1:11:26:13:20:40:0:0:0.928:0.998:0	0/0/0/2:10:19:13:20:40:0:0:0.894:0.989:0
+CHR1	31	CHR1_30_50	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=12;AC=.;NS=3;DP=.;RCOUNT=72;END=50;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:24:0:0:.:1:1:0	0/0/0/0:60:60:.:24:0:0:.:1:1:0	0/0/0/0:60:60:.:24:0:0:.:1:1:0
+CHR2	11	CHR2_10_30	AAAAAAAAAAAAAAAAAAAA	AAAAAAAAAGAAAAAAAAAA,AAAAAAAAATAAAAAAAAAA,AAAATAAAAGAAAAAAAAAA	.	PASS	AN=12;AC=3,2,1;NS=3;DP=42;RCOUNT=72;END=30;NVAR=2;SNVPOS=5,10	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/2:7:15:14:24:28:0:0:0.821:0.969:0	0/0/1/2:10:21:14:24:28:0:0:0.9:0.991:0	0/1/1/3:7:21:14:24:28:0:0:0.806:0.992:0
+CHR3	21	CHR3_20_40	AAAAAAAAAAAAAAAAAAAA	.	.	PASS	AN=12;AC=.;NS=3;DP=.;RCOUNT=0;END=40;NVAR=0;SNVPOS=.	GT:GQ:PHQ:DP:RCOUNT:RCALLS:MEC:MECP:GPM:PHPM:MCI	0/0/0/0:60:60:.:0:0:0:.:1:1:0	0/0/0/0:60:60:.:0:0:0:.:1:1:0	0/0/0/0:60:60:.:0:0:0:.:1:1:0