Merge branch 'include_cnvs'

src/main/python/clinicalfilter/inheritance.py

@@ -148,12 +148,10 @@ def examine_variant(self, variant, inheritance):
         """
 
         if variant.is_cnv():
-            self.log_string = "skipping CNVs for now"
-            return "nothing"
-            # cnv_checker = CNVInheritance(variant, self.trio, self.known_genes, self.cnv_regions)
-            # check = cnv_checker.check_single_inheritance()
-            # self.log_string = cnv_checker.log_string
-            # return check
+            cnv_checker = CNVInheritance(variant, self.trio, self.known_genes, self.cnv_regions)
+            check = cnv_checker.check_single_inheritance()
+            self.log_string = cnv_checker.log_string
+            return check
 
         if not self.trio.has_parents():
             return self.check_variant_without_parents(inheritance)
@@ -443,7 +441,7 @@ def check_single_inheritance(self):
 
         # check that the inheritance status is consistent with the parental
         # affected status
-        inh = self.variant.child.format["INHERITANCE"]
+        inh = [self.variant.child.format["INHERITANCE"], self.variant.child.format["CIFER_INHERITANCE"]]
         if not self.inheritance_matches_parental_affected_status(inh):
             if self.check_compound_inheritance():
                 self.log_string = "possible compound het CNV"
@@ -498,8 +496,6 @@ def check_compound_inheritance(self):
         SNVs to see if their gene includes a CNV variant.
         """
 
-        # return True
-
         # we don't want CNVs that don't have copy number of 1 or 3, since
         # copy number = 1 or 3 are the only ones that could operate as compound
         # hets (other copy numbers such as 0 are implicitly dominant)
@@ -552,25 +548,42 @@ def check_variant_without_parents(self):
         return "nothing"
 
     def inheritance_matches_parental_affected_status(self, inh):
-        """ check that the inheritance matches the parental affected status
+        """ check that the inheritance matches the parental affected status.
+        
+        If the variant has been inherited from the mother (ie maternally), we
+        expect the mother to also be affected. For some variants we don't know
+        whether how the variant was transmitted (due to uncertainties in
+        classifying the transmission). In this case, we assume the inheritance
+        state would be correct for the parental affected states.
         
         Args:
-            inh: inheritace status of a CNV, eg maternal, deNovo etc
+            inh: list of inheritance statuses of a CNV. We have two inheritance
+                classifications, from VICAR (classified from parental likelihoods
+                from array CGH data) and CIFER (classified from exome based read
+                depths in populations). This gives lists such as:
+                [maternal, maternal_inh], [not_inherited, deNovo] etc
             
         Returns:
             True/False for whether the inheritance is consistent with the
                parental affected statuses
         """
 
-        # if the inheritance status indiates that the CNV was inherited, check
-        # that the pertinent parents actually are affected.
-        if inh not in ["paternal", "maternal", "biparental", "inheritedDuo"]:
-            return True
+        # figure out whether the inheritance classifications indicate whether
+        # the variant is paternally, maternally, or biparentally inherited
+        paternal = any(["paternal" in x for x in inh])
+        maternal = any(["maternal" in x for x in inh])
+        biparental = any([y in x for x in inh for y in ["biparental", "inheritedDuo"]])
 
-        elif (inh == "paternal" and self.trio.father.is_affected()) or \
-            (inh == "maternal" and self.trio.mother.is_affected()) or \
-            ((inh == "biparental" or inh == "inheritedDuo") and \
-            (self.trio.father.is_affected() or self.trio.mother.is_affected())):
+        if not (paternal or maternal or biparental):
+            # if the variant isn't inherited (or the inheritance isn't known),
+            # then the parental affected statuses are irrelevant.
+            return True
+        elif (paternal and self.trio.father.is_affected()) or \
+              (maternal and self.trio.mother.is_affected()) or \
+              ((biparental or inh == "inheritedDuo") and \
+              (self.trio.father.is_affected() or self.trio.mother.is_affected())):
+            # if the inheritance status indiates that the CNV was inherited,
+            # the pertinent parents need to be also affected.
             return True
 
         return False

src/main/python/clinicalfilter/variant/cnv.py

@@ -5,7 +5,6 @@
 from clinicalfilter.variant.variant import Variant
 from clinicalfilter.variant.cnv_acgh_filter import ACGH_CNV
 from clinicalfilter.variant.cnv_exome_filter import ExomeCNV
-from clinicalfilter.variant.cnv_breakdancer_filter import BreakdancerCNV
 
 class CNV(Variant, VariantInfo):
     """  class to take CNV data for an individual, and
@@ -130,22 +129,23 @@ def passes_filters(self):
         except ValueError:
             return False
 
+        # we rely on the CALLSOURCE field to inform us what the CNV has been
+        # called by. Raise an error if this is not present.
+        assert "CALLSOURCE" in self.info
+
         track_variant = False
         if self.get_chrom() == self.debug_chrom and self.get_position() == self.debug_pos:
             track_variant = True
 
         passes = True
-        if "CONVEX" in self.info and "CNSOLIDATE" not in self.info:
+        if "aCGH" in self.info["CALLSOURCE"]:
+            filt = ACGH_CNV(self)
+            passes = filt.filter_cnv(track_variant)
+        elif "EXOME" in self.info["CALLSOURCE"]:
             # currently return false for all exome-only CNVs, undergoing testing
             filt = ExomeCNV(self)
             # passes = filt.filter_cnv(track_variant)
             passes = False
-        elif "CNSOLIDATE" in self.info:
-            filt = ACGH_CNV(self)
-            passes = filt.filter_cnv(track_variant)
-        elif "BREAKDANCER" in self.info:
-            filt = BreakdancerCNV(self)
-            passes = filt.filter_cnv(track_variant)
         else:
             if track_variant:
                 print("CNV is not an aCGH or exome CNV")

src/main/python/clinicalfilter/variant/cnv_acgh_filter.py

@@ -40,15 +40,26 @@ def filter_cnv(self, track_variant):
             passes = False
             if track_variant:
                 print("failed no exons", self.cnv.info["NUMBEREXONS"])
+        elif self.fails_frequency():
+            passes = False
+            if track_variant:
+                print("failed frequency", self.cnv.info["ACGH_RC_FREQ50"])
+        elif self.fails_cifer_inh():
+            passes = False
+            if track_variant:
+                print("failed CIFER inheritance", self.cnv.format["CIFER_INHERITANCE"])
 
         return passes
 
     def fails_mad_ratio(self):
-        """ checks if the MAD ratio is too low
+        """ checks if the MAD ratio is too low.
+        
+        Note that this filter has been lowered, so that it will no longer
+        exclude any variants. This function could probably be removed.
         """
 
         try:
-            return abs(float(self.cnv.info["MEANLR2"])/float(self.cnv.info["MADL2R"])) < 15
+            return abs(float(self.cnv.info["MEANLR2"])/float(self.cnv.info["MADL2R"])) < 0
         except ZeroDivisionError:
             return True
 
@@ -86,3 +97,26 @@ def fails_no_exons(self):
         """
 
         return float(self.cnv.info["NUMBEREXONS"]) < 1
+
+    def fails_frequency(self):
+        """ checks that the CNV has a low population frequency.
+        
+        If the CNV has been observed in the unaffected controls (aka unaffected
+        parents), then we can determine the population frequency, which must be
+        sufficiently rare to pass.
+        
+        If the population frequency field is absent, assume the frequency is 0.
+        """
+
+        try:
+            return int(self.cnv.info["ACGH_RC_FREQ50"] > 0.01)
+        except KeyError:
+            # If the field isn't available, assume the frequency is 0.
+            return False
+
+    def fails_cifer_inh(self):
+        """ check that the CIFER inheritance classification isn't false_positive
+        """
+
+        return self.cnv.format["CIFER_INHERITANCE"] == "false_positive"
+

src/test/python/clinicalfilter/test_cnv_inheritance.py

@@ -33,7 +33,7 @@ def setUp(self):
 
         self.inh = CNVInheritance(self.variant, self.trio, self.known_genes, syndrome_regions)
 
-    def create_cnv(self, gender, inh, chrom, pos):
+    def create_cnv(self, gender, inh, cifer, chrom, pos):
         """ create a default variant
         """
 
@@ -46,8 +46,8 @@ def create_cnv(self, gender, inh, chrom, pos):
         var = CNV(chrom, pos, snp_id, ref, alt, filt)
 
         info = "HGNC=TEST;HGNC_ALL=TEST;END=16000000;SVLEN=5000"
-        format_keys = "INHERITANCE:DP"
-        sample_values = inh + ":50"
+        format_keys = "CIFER:INHERITANCE:DP"
+        sample_values = cifer + ":" + inh + ":50"
 
         var.add_info(info)
         var.add_format(format_keys, sample_values)
@@ -61,9 +61,9 @@ def create_variant(self, child_gender, chrom="1", position="15000000"):
         """
 
         # generate a test variant
-        child_var = self.create_cnv(child_gender, "unknown", chrom, position)
-        mom_var = self.create_cnv("F", "unknown", chrom, position)
-        dad_var = self.create_cnv("M", "unknown", chrom, position)
+        child_var = self.create_cnv(child_gender, "unknown", "uncertain", chrom, position)
+        mom_var = self.create_cnv("F", "unknown", "uncertain", chrom, position)
+        dad_var = self.create_cnv("M", "unknown", "uncertain", chrom, position)
 
         var = TrioGenotypes(child_var)
         var.add_mother_variant(mom_var)
@@ -89,20 +89,34 @@ def test_inheritance_matches_parental_affected_status(self):
 
         # check that paternally inherited CNVs that have affected fathers pass
         self.inh.trio.father.affected_status = "2"
-        inh = "paternal"
+        inh = ["paternal"]
+        self.assertTrue(self.inh.inheritance_matches_parental_affected_status(inh))
+
+        # check for a CIFER derived annotation
+        inh = ["paternal_inh"]
+        self.assertTrue(self.inh.inheritance_matches_parental_affected_status(inh))
+
+        # check for a list of the VICAR and CIFER inheritance annotations
+        inh = ["paternal", "paternal_inh"]
+        self.assertTrue(self.inh.inheritance_matches_parental_affected_status(inh))
+
+        # check when one annotation says inherited, but the other doesn't.
+        # Having one parentally inherited annotation is sufficient to classify
+        # the CNV as inherited.
+        inh = ["uknown", "paternal_inh"]
         self.assertTrue(self.inh.inheritance_matches_parental_affected_status(inh))
 
         # check that paternally inherited CNVs without an affected father fail
         self.inh.trio.father.affected_status = "1"
         self.assertFalse(self.inh.inheritance_matches_parental_affected_status(inh))
 
         # check that maternally inherited CNVs without an affected mother fail
-        inh = "maternal"
+        inh = ["maternal"]
         self.inh.trio.father.affected_status = "1"
         self.assertFalse(self.inh.inheritance_matches_parental_affected_status(inh))
 
         # check that biparentally inherited CNVs pass if either parent is affected
-        inh = "biparental"
+        inh = ["biparental"]
         self.assertFalse(self.inh.inheritance_matches_parental_affected_status(inh))
         self.inh.trio.father.affected_status = "2"
         self.assertTrue(self.inh.inheritance_matches_parental_affected_status(inh))
@@ -111,7 +125,7 @@ def test_inheritance_matches_parental_affected_status(self):
         self.assertTrue(self.inh.inheritance_matches_parental_affected_status(inh))
 
         # check that biparentally inherited CNVs pass if either parent is affected
-        inh = "inheritedDuo"
+        inh = ["inheritedDuo"]
         self.inh.trio.mother.affected_status = "1"
         self.assertFalse(self.inh.inheritance_matches_parental_affected_status(inh))
         self.inh.trio.father.affected_status = "2"
@@ -121,7 +135,7 @@ def test_inheritance_matches_parental_affected_status(self):
         self.assertTrue(self.inh.inheritance_matches_parental_affected_status(inh))
 
         # check that noninherited CNVs pass, regardless of parental affected status
-        inh = "deNovo"
+        inh = ["deNovo"]
         self.inh.trio.mother.affected_status = "1"
         self.assertTrue(self.inh.inheritance_matches_parental_affected_status(inh))
 

src/test/python/clinicalfilter/test_variant_cnv_acgh_filter.py

@@ -34,24 +34,19 @@ def setUp(self):
         self.var.cnv.set_gender("F")
 
     def test_fails_mad_ratio(self):
-        """ test that fails_mad_ratio() works correctly
+        """ test that fails_mad_ratio() works correctly.
         """
 
-        # check that var passes when MAD ratio > 15
+        # check that var passes when MAD ratio > 0
         self.var.cnv.info["MEANLR2"] = "0.5"
         self.var.cnv.info["MADL2R"] = "0.02"
         self.assertFalse(self.var.fails_mad_ratio())
 
-        # check that var passes when MAD ratio == 15
+        # check that var passes when MAD ratio == 0
         self.var.cnv.info["MEANLR2"] = "0.3"
-        self.var.cnv.info["MADL2R"] = "0.02"
+        self.var.cnv.info["MADL2R"] = float("inf")
         self.assertFalse(self.var.fails_mad_ratio())
 
-        # check that var fails when MAD ratio < 15
-        self.var.cnv.info["MEANLR2"] = "0.2"
-        self.var.cnv.info["MADL2R"] = "0.02"
-        self.assertTrue(self.var.fails_mad_ratio())
-
         # check that var fails when trying to divide by zero
         self.var.cnv.info["MEANLR2"] = "0.2"
         self.var.cnv.info["MADL2R"] = "0"