vcfallelicprimitives: added shell output support for tests and show single example

Author: pjotrp

samples/10158243.vcf

@@ -0,0 +1,13 @@
+##fileformat=VCFv4.2
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##INFO=<ID=CONFLICT,Number=.,Type=String,Description="Sample names for which there are multiple paths in the graph with conflicting alleles">
+##INFO=<ID=AC,Number=A,Type=Integer,Description="Total number of alternate alleles in called genotypes">
+##INFO=<ID=AF,Number=A,Type=Float,Description="Estimated allele frequency in the range (0,1]">
+##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of samples with data">
+##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
+##INFO=<ID=LV,Number=1,Type=Integer,Description="Level in the snarl tree (0=top level)">
+##INFO=<ID=PS,Number=1,Type=String,Description="ID of variant corresponding to parent snarl">
+##INFO=<ID=AT,Number=R,Type=String,Description="Allele Traversal as path in graph">
+##contig=<ID=grch38#chr4,length=97395>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	HG00438	HG00621	HG00673	HG00733	HG00735	HG00741	HG01071	HG01106	HG01109	HG01123	HG01175	HG01243	HG01258	HG01358	HG01361	HG01891	HG01928	HG01952	HG01978	HG02055	HG02080	HG02109	HG02145	HG02148	HG02257	HG02486	HG02559	HG02572	HG02622	HG02630	HG02717	HG02723	HG02818	HG02886	HG03098	HG03453	HG03486	HG03492	HG03516	HG03540	HG03579	NA18906	NA20129	NA21309	chm13
+grch38#chr4	10158243	>3655>3662	ACCCCCACCCCCACC	ACC,AC,ACCCCCACCCCCAC,ACCCCCACC,ACA	60	.	AC=64,3,2,3,1;AF=0.719101,0.0337079,0.0224719,0.0337079,0.011236;AN=89;AT=>3655>3656>3657>3658>3659>3660>3662,>3655>3656>3660>3662,>3655>3660>3662,>3655>3656>3657>3658>3660>3662,>3655>3656>3657>3660>3662,>3655>3656>3661>3662;NS=45;LV=0	GT	0|0	1|1	1|1	1|0	5|1	0|4	0|1	0|1	1|1	1|1	1|1	1|1	1|1	1|1	1|1	4|3	1|1	1|1	1|1	1|0	1|0	1|0	1|0	1|1	1|1	1|4	1|1	1|1	3|0	1|0	1|1	0|1	1|1	1|1	2|1	1|2	1|1	1|1	0|1	1|1	1|1	1|0	1|2	1|1	0

…ta/regression/vcfallelicprimitives_2.vcf …ta/regression/vcfallelicprimitives_4.vcftest/data/regression/vcfallelicprimitives_2.vcf renamed to test/data/regression/vcfallelicprimitives_4.vcf test/data/regression/vcfallelicprimitives_2.vcf renamed to test/data/regression/vcfallelicprimitives_4.vcf

@@ -36,6 +36,18 @@ def head(cmd, lines=4):
     header = [l.replace("../build/", "") for l in header]
     print("\n".join(header))
 
+def sh(cmd):
+    cmds = ['bash','-c',cmd]
+    p = Popen(cmds, stdout=PIPE, stderr=PIPE, close_fds=True)
+    output = p.communicate()
+    out = output[0]
+    if len(out) == 0:
+        # if stdout is empty fetch stderr
+        out = output[1]
+    header = out.decode().split("\n")
+    header = out.decode().expandtabs(tabsize=8).split("\n")
+    print("\n".join(header))
+
 def run_stdout(cmd, ext = "vcf"):
     os.makedirs(tmpdir,exist_ok=True)
     curframe = inspect.currentframe()

…ta/regression/vcfallelicprimitives_3.vcf …ta/regression/vcfallelicprimitives_5.vcftest/data/regression/vcfallelicprimitives_3.vcf renamed to test/data/regression/vcfallelicprimitives_5.vcf test/data/regression/vcfallelicprimitives_3.vcf renamed to test/data/regression/vcfallelicprimitives_5.vcf

@@ -25,6 +25,8 @@ multiple SNPs unless the -m flag is provided.
 
 : shows help message and exits.
 
+See more below.
+
 # EXIT VALUES
 
 **0**
@@ -38,19 +40,54 @@ multiple SNPs unless the -m flag is provided.
 
 <!--
 
-    >>> from pytest.rtest import run_stdout, head, cat
+    >>> from pytest.rtest import run_stdout, head, cat, sh
 
 -->
 
 ```
 
->>> head("vcfallelicprimitives -h",1)
+>>> head("vcfallelicprimitives -h",20)
 usage: vcfallelicprimitives [options] [file]
+>
+If multiple allelic primitives (gaps or mismatches) are specified in
+a single VCF record, split the record into multiple lines, but drop all
+INFO fields.  Does not handle genotypes (yet).  MNPs are split into
+multiple SNPs unless the -m flag is provided.  Records generated by splits have the
+options:
+    -m, --use-mnps          Retain MNPs as separate events (default: false).
+    -t, --tag-parsed FLAG   Tag records which are split apart of a complex allele with this flag.
+    -L, --max-length LEN    Do not manipulate records in which either the ALT or
+                            REF is longer than LEN (default: 200).
+    -k, --keep-info         Maintain site and allele-level annotations when decomposing.
+                            Note that in many cases, such as multisample VCFs, these won't
+                            be valid post-decomposition.  For biallelic loci in single-sample
+                            VCFs, they should be usable with caution.
+    -g, --keep-geno         Maintain genotype-level annotations when decomposing.  Similar
+                            caution should be used for this as for --keep-info.
+>
+Type: transformation
+>
 
 ```
 
-vcfallelicprimitives picks complex regions and simplifies nested alignments
+vcfallelicprimitives picks complex regions and simplifies nested alignments. For example:
+
+```python
+
+>>> sh("grep 10158243 ../samples/10158243.vcf")
+grch38#chr4     10158243        >3655>3662      ACCCCCACCCCCACC ACC,AC,ACCCCCACCCCCAC,ACCCCCACC,ACA     60      .       AC=64,3,2,3,1;AF=0.719101,0.0337079,0.0224719,0.0337079,0.011236;AN=89;AT=>3655>3656>3657>3658>3659>3660>3662,>3655>3656>3660>3662,>3655>3660>3662,>3655>3656>3657>3658>3660>3662,>3655>3656>3657>3660>3662,>3655>3656>3661>3662;NS=45;LV=0     GT      0|0     1|1     1|1     1|0     5|1     0|4     0|1     0|1     1|1     1|1     1|1     1|1     1|1     1|1     1|1     4|3     1|1     1|1     1|1     1|0     1|0     1|0     1|0     1|1     1|1     1|4     1|1     1|1     3|0     1|0     1|1     0|1     1|1     1|1     2|1     1|2     1|1     1|1     0|1     1|1     1|1     1|0     1|2     1|1     0
+
+```
 
+After aligning it reduces into two records and adjusts the genotypes accordingly:
+
+```python
+
+>>> sh("../build/vcfallelicprimitives -m -L 1000 ../samples/10158243.vcf|grep -v ^\#")
+grch38#chr4     10158243        .       ACCCCCACCCCCAC  ACCCCCAC,ACAC,AC,A      60      .       AC=3,1,64,3;AF=0.0337079,0.011236,0.719101,0.0337079;LEN=6,10,12,13;TYPE=del,del,del,del        GT      0|0     3|3     3|3     3|0     2|3     0|1     0|3     0|3     3|3     3|3     3|3     3|3     3|3     3|3     3|3     1|0     3|3     3|3     3|3     3|0     3|0     3|0     3|0     3|3     3|3     3|1     3|3     3|3     0|0     3|0     3|3     0|3     3|3     3|3     4|3     3|4     3|3     3|3     0|3     3|3     3|3     3|0     3|4     3|3     0
+grch38#chr4     10158255        .       ACC     AC,A    60      .       AC=2,1;AF=0.0224719,0.011236;LEN=1,2;TYPE=del,del       GT      0|0     0|0     0|0     0|0     2|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|1     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     1|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0|0     0
+
+```
 
 ## Source code
 
@@ -60,10 +97,10 @@ vcfallelicprimitives picks complex regions and simplifies nested alignments
 
 ```python
 >>> run_stdout("vcfallelicprimitives -m -L 1000 ../samples/grch38#chr8_36353854-36453166.vcf", ext="vcf")
-output in <a href="../data/regression/vcfallelicprimitives_2.vcf">vcfallelicprimitives_2.vcf</a>
+output in <a href="../data/regression/vcfallelicprimitives_4.vcf">vcfallelicprimitives_4.vcf</a>
 
 >>> run_stdout("vcfallelicprimitives -m -L 1000 ../samples/grch38#chr4_10083863-10181258.vcf", ext="vcf")
-output in <a href="../data/regression/vcfallelicprimitives_3.vcf">vcfallelicprimitives_3.vcf</a>
+output in <a href="../data/regression/vcfallelicprimitives_5.vcf">vcfallelicprimitives_5.vcf</a>
 
 ```