Mercurial > repos > ulfschaefer > filter_vcf
comparison phe/variant/MPileupVariantCaller.py @ 0:834a312c0114 draft
Uploaded
author | ulfschaefer |
---|---|
date | Thu, 10 Dec 2015 09:22:39 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:834a312c0114 |
---|---|
1 ''' | |
2 Created on 22 Sep 2015 | |
3 | |
4 @author: alex | |
5 ''' | |
6 from collections import OrderedDict | |
7 import logging | |
8 import os | |
9 import subprocess | |
10 import tempfile | |
11 | |
12 from phe.variant import VariantCaller | |
13 | |
14 | |
15 class MPileupVariantCaller(VariantCaller): | |
16 """Implemetation of the Broad institute's variant caller.""" | |
17 | |
18 name = "mpileup" | |
19 """Plain text name of the variant caller.""" | |
20 | |
21 _default_options = "-m -f GQ" | |
22 """Default options for the variant caller.""" | |
23 | |
24 def __init__(self, cmd_options=None): | |
25 """Constructor""" | |
26 if cmd_options is None: | |
27 cmd_options = self._default_options | |
28 | |
29 super(MPileupVariantCaller, self).__init__(cmd_options=cmd_options) | |
30 | |
31 self.last_command = None | |
32 | |
33 def get_info(self, plain=False): | |
34 d = {"name": self.name, "version": self.get_version(), "command": self.last_command} | |
35 | |
36 if plain: | |
37 result = "mpileup(%(version)s): %(command)s" % d | |
38 else: | |
39 result = OrderedDict(d) | |
40 | |
41 return result | |
42 | |
43 def get_version(self): | |
44 | |
45 p = subprocess.Popen(["samtools", "--version"], stdout=subprocess.PIPE, stderr=subprocess.STDOUT) | |
46 (output, _) = p.communicate() | |
47 | |
48 # first line is the version of the samtools | |
49 version = output.split("\n")[0].split(" ")[1] | |
50 | |
51 return version | |
52 | |
53 def make_vcf(self, *args, **kwargs): | |
54 ref = kwargs.get("ref") | |
55 bam = kwargs.get("bam") | |
56 | |
57 if kwargs.get("vcf_file") is None: | |
58 kwargs["vcf_file"] = "variants.vcf" | |
59 | |
60 opts = {"ref": os.path.abspath(ref), | |
61 "bam": os.path.abspath(bam), | |
62 "all_variants_file": os.path.abspath(kwargs.get("vcf_file")), | |
63 "extra_cmd_options": self.cmd_options} | |
64 | |
65 with tempfile.NamedTemporaryFile(suffix=".pileup") as tmp: | |
66 opts["pileup_file"] = tmp.name | |
67 cmd = "samtools mpileup -t DP,DV,DP4,DPR,SP -Auf %(ref)s %(bam)s | bcftools call %(extra_cmd_options)s > %(all_variants_file)s" % opts | |
68 print cmd | |
69 self.last_command = cmd | |
70 if os.system(cmd) != 0: | |
71 logging.warn("Pileup creation was not successful.") | |
72 return False | |
73 | |
74 return True | |
75 | |
76 def create_aux_files(self, ref): | |
77 """Index reference with faidx from samtools. | |
78 | |
79 Parameters: | |
80 ----------- | |
81 ref: str | |
82 Path to the reference file. | |
83 | |
84 Returns: | |
85 -------- | |
86 bool: | |
87 True if auxiliary files were created, False otherwise. | |
88 """ | |
89 | |
90 success = os.system("samtools faidx %s" % ref) | |
91 | |
92 if success != 0: | |
93 logging.warn("Fasta index could not be created.") | |
94 return False | |
95 | |
96 return True |