Mercurial > repos > ulfschaefer > filter_vcf

comparison phe/variant_filters/DP4Filter.py @ 0:834a312c0114 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author ulfschaefer
date Thu, 10 Dec 2015 09:22:39 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:834a312c0114
1 '''Filter VCFs on AD ratio.
2
3 Created on 24 Sep 2015
4
5 @author: alex
6 '''
7
8 import argparse
9 import logging
10
11 from phe.variant_filters import PHEFilterBase
12
13
14 class DP4Filter(PHEFilterBase):
15 '''Filter sites by AD ratio.'''
16
17
18 name = "DP4"
19 _default_threshold = 0.9
20 parameter = "dp4_ratio"
21
22 @classmethod
23 def customize_parser(self, parser):
24 arg_name = self.parameter.replace("_", "-")
25 parser.add_argument("--%s" % arg_name, type=float, default=self._default_threshold,
26 help="Filter sites below minimum dp4 ratio (default: %s)" % self._default_threshold)
27
28 def __init__(self, args):
29 """AD Ratio constructor."""
30 # This needs to happen first, because threshold is initialised here.
31 super(DP4Filter, self).__init__(args)
32
33 # Change the threshold to custom dp value.
34 self.threshold = self._default_threshold
35 if isinstance(args, argparse.Namespace):
36 self.threshold = args.ad_ratio
37 elif isinstance(args, dict):
38 try:
39 self.threshold = float(args.get(self.parameter))
40 except TypeError:
41 logging.error("Could not retrieve threshold from %s", args.get(self.parameter))
42 self.threshold = None
43
44
45 def __call__(self, record):
46 """Filter a :py:class:`vcf.model._Record`."""
47
48 if not record.is_snp:
49 return None
50
51 try:
52 record_dp = record.INFO.get("DP4")
53
54 # FIXME: when record length is > 2, what do you do?
55 assert len(record_dp) == 4, "DP4 data should have 4 datum POS: %i" % record.POS
56
57 depth = sum(record_dp)
58
59 ratio = float(sum(record_dp[2:])) / depth
60 except Exception:
61 logging.error("Could not calculate DP4 ratio from %s POS: %s", record_dp, record.POS)
62 ratio = None
63
64 if ratio is None or ratio < self.threshold:
65 # FIXME: When ratio is None, i.e. error, what do you do?
66 return ratio or False
67 else:
68 return None
69
70 def short_desc(self):
71 short_desc = self.__doc__ or ''
72
73 if short_desc:
74 short_desc = "%s (DP4 ratio > %s )" % (short_desc, self.threshold)
75
76 return short_desc