Repository: bihealth/vcfpy

Issue:
add FORMAT and calls to vcf?
vcfpy version: 0.13.2
Python version: 3.6
Operating System: Centos 7

### Question

I would like to modify a vcf (from Lofreq) that does not contain FORMAT or call columns in order for it to be accepted by Picard within the bcbio-nextgen framework. 

I have tried like this but it does not work. Nothing comes out:

```
reader = vcfpy.Reader.from_path("/path/to/sample.vcf.gz")
reader.header.samples = vcfpy.SamplesInfos(names=['D20-XXX'], name_to_idx={'D20-XXX': 0})

writer = vcfpy.Writer.from_path('/dev/stdout', reader.header)

for rec in reader:
    data = {'GT': '1/1'}
    rec.calls.append(vcfpy.Call(sample=sample, data=data))
    rec.FORMAT = []
    rec.add_format('GT')
    writer.write_record(rec)
```

Is there some inherent limitation that prevents adding this info? I noticed a caveat that the FORMAT field cannot be altered after modifying the SamplesInfos object.

Incidently, I have also noticed that data objects in the call instances are dicts and not ordered dicts, as the documentation says. Perhaps a remnant from earlier times?




Discussion:

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index f3109ae..c20e5f3 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -5,6 +5,7 @@ History
 HEAD
 ----
 
+- Adding ``Record.update_calls``.
 - Making ``Record.{format,calls}`` use list when empty
 
 v0.13.2 (2020-08-20)
diff --git a/vcfpy/record.py b/vcfpy/record.py
index 4c89104..4dad26f 100644
--- a/vcfpy/record.py
+++ b/vcfpy/record.py
@@ -86,11 +86,16 @@ class Record:
         self.FORMAT = FORMAT or []
         #: A list of genotype :py:class:`Call` objects.  Optional, must be given if
         #: and only if ``FORMAT`` is also given.
-        self.calls = list(calls or [])
-        for call in self.calls:
-            call.site = self
+        self.calls = calls or []
         #: A mapping from sample name to entry in self.calls.
-        self.call_for_sample = {call.sample: call for call in self.calls}
+        self.call_for_sample = {}
+        self.update_calls(self.calls)
+
+    def update_calls(self, calls):
+        """Update ``self.calls`` and other fields as necessary."""
+        for call in calls:
+            call.site = self
+        self.call_for_sample = {call.sample: call for call in calls}
 
     def is_snv(self):
         """Return ``True`` if it is a SNV"""


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 429a80e..b346875 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -67,3 +67,30 @@ def multisample_vcf_file(tmpdir, multisample_vcf):
     with open(str(p), "wt") as f:
         f.write(multisample_vcf)
     return str(p)
+
+
+@pytest.fixture
+def nosample_vcf():
+    """Return string VCF that has no samples and no FORMAT"""
+    return textwrap.dedent(
+        """
+        ##fileformat=VCFv4.3
+        ##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>
+        ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+        #CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO
+        20\t14370\t.\tG\tA\t29\t.\t.
+        20\t17330\t.\tT\tA\t3\t.\t.
+        20\t1110696\t.\tA\tG,T\t67\t.\t.
+        20\t1230237\t.\tT\t.\t47\t.\t.
+        20\t1234567\t.\tGTC\tG,GTCT\t50\t.\t.
+        """
+    ).lstrip()
+
+
+@pytest.fixture
+def nosample_vcf_file(tmpdir, nosample_vcf):
+    """Return path to file without samples/FORMAT"""
+    p = tmpdir.mkdir("input").join("input.vcf")
+    with open(str(p), "wt") as f:
+        f.write(nosample_vcf)
+    return str(p)
diff --git a/tests/test_reader_nosample.py b/tests/test_reader_nosample.py
new file mode 100644
index 0000000..e982b93
--- /dev/null
+++ b/tests/test_reader_nosample.py
@@ -0,0 +1,50 @@
+# -*- coding: utf-8 -*-
+"""Tests for reading with file without samples"""
+
+from io import StringIO
+import textwrap
+
+import vcfpy
+
+
+def test_reading_parse_nosample(tmpdir, nosample_vcf_file):
+    """Test reading of VCF file without any samples."""
+    # Perform record-wise copying, saving results in records
+    records = []
+    with vcfpy.Reader.from_path(nosample_vcf_file) as reader:
+        for record in reader:
+            records.append(record)
+    # Check resulting records, checking the first and last records is enough
+    assert len(records) == 5
+
+
+def test_reading_parse_nosample(tmpdir, nosample_vcf_file):
+    """Read VCF file without samples, write file with samples."""
+    # Perform record-wise copying, saving results in records
+    path_out = tmpdir.mkdir("output").join("output.vcf")
+    with vcfpy.Reader.from_path(nosample_vcf_file) as reader:
+        header = reader.header.copy()
+        header.samples = vcfpy.SamplesInfos(["NA00001", "NA00002", "NA00003"])
+        with vcfpy.Writer.from_path(str(path_out), header) as writer:
+            for record in reader:
+                record.update_calls(
+                    [vcfpy.Call(sample, {}) for sample in ("NA00001", "NA00002", "NA00003")]
+                )
+                record.add_format("GT", "./.")
+                writer.write_record(record)
+
+    expected = textwrap.dedent(
+        """
+    ##fileformat=VCFv4.3
+    ##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>
+    ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+    #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NA00001	NA00002	NA00003
+    20	14370	.	G	A	29	.	.	GT	.	.	.
+    20	17330	.	T	A	3	.	.	GT	.	.	.
+    20	1110696	.	A	G,T	67	.	.	GT	.	.	.
+    20	1230237	.	T	.	47	.	.	GT	.	.	.
+    20	1234567	.	GTC	G,GTCT	50	.	.	GT	.	.	.
+    """
+    ).lstrip()
+
+    assert path_out.open("rt").read() == expected