Repository: brentp/cyvcf2 (language: python)

Issue:
Use Python's logging module instead of writing directly to stderr
There are a few occasions where cyvcf2.pyx writes to stderr. It would be nice to be able to capture and/or filter these in some cases. Would it be ok if these were replaced with calls to the [`logging`](https://docs.python.org/3/library/logging.html) module? There seem to be some caveats with using the logger inside cython, but I think it should be ok if logger attributes related to stack frames aren't used.

The specific message I want to silence right now is "no intervals found for ...", because I'm deliberately doing `next(vcf(contig_id))` for each contig id found in the header to see if any variants exist for that contig (I then create a mapping from contig id to vcf file, which is handy for vcfs split by chromosome).
https://github.com/brentp/cyvcf2/blob/2addcaf097bba1b6fa219a260e6ec39a2b2e936c/cyvcf2/cyvcf2.pyx#L403
https://github.com/brentp/cyvcf2/blob/2addcaf097bba1b6fa219a260e6ec39a2b2e936c/cyvcf2/cyvcf2.pyx#L457


Discussion / hints:
> Would it be ok if these were replaced with calls to the logging module

Yes, that would be great!

Fix patch:
diff --git a/cyvcf2/cyvcf2.pyx b/cyvcf2/cyvcf2.pyx
index 2fc7399..5875d77 100644
--- a/cyvcf2/cyvcf2.pyx
+++ b/cyvcf2/cyvcf2.pyx
@@ -12,6 +12,7 @@ import numpy as np
 from array import array
 import math
 import ctypes
+import warnings
 try:
   from pathlib import Path
 except ImportError:
@@ -350,12 +351,12 @@ cdef class VCF(HTSFile):
 
         ret = bcf_hdr_set_samples(self.hdr, <const char *>samples, 0)
         assert ret >= 0, ("error setting samples", ret)
+        self.n_samples = bcf_hdr_nsamples(self.hdr)
         if ret != 0 and samples != "-":
             s = from_bytes(samples).split(",")
-            if ret < len(s):
-                sys.stderr.write("warning: not all requested samples found in VCF\n")
+            if self.n_samples < len(s):
+                warnings.warn("not all requested samples found in VCF")
 
-        self.n_samples = bcf_hdr_nsamples(self.hdr)
 
     def update(self, id, type, number, description):
         """Update the header with an INFO field of the given parameters.
@@ -400,7 +401,7 @@ cdef class VCF(HTSFile):
 
         itr = bcf_itr_querys(self.hidx, self.hdr, to_bytes(region))
         if itr == NULL:
-            sys.stderr.write("no intervals found for %s at %s\n" % (self.fname, region))
+            warnings.warn("no intervals found for %s at %s" % (self.fname, region))
             raise StopIteration
         try:
             while True:
@@ -409,10 +410,8 @@ cdef class VCF(HTSFile):
                 if ret < 0:
                     bcf_destroy(b)
                     break
-                if bcf_subset_format(self.hdr, b) != 0:
-                    sys.stderr.write("could not subset variant")
-                    bcf_destroy(b)
-                    break
+                ret = bcf_subset_format(self.hdr, b)
+                assert ret == 0, ("could not subset variant", self.fname, region)
                 yield newVariant(b, self)
         finally:
             if itr != NULL:
@@ -454,7 +453,7 @@ cdef class VCF(HTSFile):
             itr = tbx_itr_querys(self.idx, cregion)
 
         if itr == NULL:
-            sys.stderr.write("no intervals found for %s at %s\n" % (self.fname, region))
+            warnings.warn("no intervals found for %s at %s" % (self.fname, region))
             raise StopIteration
 
         try:
@@ -905,7 +904,7 @@ cdef class VCF(HTSFile):
             nv += 1
             if nv == n_variants:
                 break
-        sys.stderr.write("tested: %d variants out of %d\n" % (nv, nvt))
+        warnings.warn("tested: %d variants out of %d" % (nv, nvt))
         return self._relatedness_finish(ibs, n, hets)
 
     cdef dict _relatedness_finish(self,
@@ -929,7 +928,7 @@ cdef class VCF(HTSFile):
         for sj in range(ns):
             sample_j = samples[sj]
             if _hets[sj] == 0:
-                print("peddy: no hets found for sample %s\n" % sample_j, file=sys.stderr)
+                warnings.warn("peddy: no hets found for sample %s" % sample_j)
             for sk in range(sj, ns):
                 if sj == sk: continue
                 sample_k = samples[sk]
@@ -1898,9 +1897,7 @@ cdef class Variant(object):
                 new_rid = bcf_hdr_id2int(self.vcf.hdr, BCF_DT_CTG, new_chrom.encode())
                 if new_rid < 0:
                     raise ValueError("Unable to add {} to CHROM".format(new_chrom))
-                sys.stderr.write(
-                    "[cyvcf2]: added new contig {} to header".format(new_chrom)
-                )
+                warnings.warn("added new contig {} to header".format(new_chrom))
             self.b.rid = new_rid
 
     property var_type:


Test patch:
diff --git a/cyvcf2/tests/test_reader.py b/cyvcf2/tests/test_reader.py
index 0aca3f9..4549085 100644
--- a/cyvcf2/tests/test_reader.py
+++ b/cyvcf2/tests/test_reader.py
@@ -11,6 +11,7 @@
   from pathlib import Path
 except ImportError:
   from pathlib2 import Path  # python 2 backport
+import warnings
 
 
 HERE = os.path.dirname(__file__)
@@ -87,11 +88,17 @@ def test_format_str():
 
 def test_missing_samples():
     samples = ['101976-101976', 'sample_not_in_vcf']
-    vcf = VCF(VCF_PATH, gts012=True, samples=samples)
+    with warnings.catch_warnings(record=True) as w:
+        vcf = VCF(VCF_PATH, gts012=True, samples=samples)
+        assert len(w) == 1
+        assert "not all requested samples found" in str(w[-1].message)
     assert len(vcf.samples) == 1
     vcf.close()
     samples = '101976-101976,sample_not_in_vcf'
-    vcf = VCF(VCF_PATH, gts012=True, samples=samples)
+    with warnings.catch_warnings(record=True) as w:
+        vcf = VCF(VCF_PATH, gts012=True, samples=samples)
+        assert len(w) == 1
+        assert "not all requested samples found" in str(w[-1].message)
     assert len(vcf.samples) == 1
 
 def test_ibd():
@@ -104,7 +111,11 @@ def test_ibd():
 
 def test_relatedness():
     vcf = VCF(VCF_PATH, gts012=True)
-    df = vcf.relatedness(gap=0, linkage_max=2)
+    with warnings.catch_warnings(record=True) as w:
+        warnings.simplefilter("ignore", category=DeprecationWarning)
+        df = vcf.relatedness(gap=0, linkage_max=2)
+        assert len(w) == 1
+        assert "tested:" in str(w[-1].message)
     assert "ibs0" in df, df
     assert "rel" in df
     #vcf = VCF(VCF_PATH, gts012=True)
@@ -643,13 +654,24 @@ def test_bcf():
     viter = vcf("1:69260-69438")
     sys.stderr.write("\nOK\n")
     sys.stderr.flush()
-    l = list(viter)
+    with warnings.catch_warnings(record=True) as w:
+        l = list(viter)
+        assert len(w) == 1
+        assert "no intervals found" in str(w[-1].message)
     assert len(l) == 0, len(l)
 
     iter = vcf("chr1:69260-69438")
     l = list(iter)
     assert len(l) == 2, len(l)
 
+def test_vcf_no_intervals():
+    vcf = VCF('{}/test.vcf.gz'.format(HERE))
+    viter = vcf("not_a_chrom")
+    with warnings.catch_warnings(record=True) as w:
+        l = list(viter)
+        assert len(w) == 1
+        assert "no intervals found" in str(w[-1].message)
+    assert len(l) == 0, len(l)
 
 def test_issue12():
     fields = "ADP_ALL ADPD ADPO ADP_PASS ADPR AFR AMBIG BMF_PASS BMF_QUANT AF_FAILED FA_FAILED FM_FAILED FP_FAILED FR_FAILED MD_FAILED IMPROPER MQ_FAILED OVERLAP PV_FAILED QSS".split()
@@ -1013,7 +1035,10 @@ def test_set_chrom_when_contig_not_in_header():
     assert new_chrom not in original_seqnames
     v = next(vcf)
 
-    v.CHROM = new_chrom
+    with warnings.catch_warnings(record=True) as w:
+        v.CHROM = new_chrom
+        assert len(w) == 1
+        assert "added new contig" in str(w[-1].message)
     assert v.CHROM == new_chrom
     expected_seqnames = sorted(original_seqnames + [new_chrom])
     assert vcf.seqnames == expected_seqnames