Mercurial > repos > ulfschaefer > vcfs2fasta
comparison phe/variant/MPileupVariantCaller.py @ 22:96f393ad7fc6 draft default tip
Uploaded
author | ulfschaefer |
---|---|
date | Wed, 23 Dec 2015 04:50:58 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
21:b09ffe50c378 | 22:96f393ad7fc6 |
---|---|
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 make_aux = kwargs.get("make_aux", False) | |
58 | |
59 if kwargs.get("vcf_file") is None: | |
60 kwargs["vcf_file"] = "variants.vcf" | |
61 | |
62 opts = {"ref": os.path.abspath(ref), | |
63 "bam": os.path.abspath(bam), | |
64 "all_variants_file": os.path.abspath(kwargs.get("vcf_file")), | |
65 "extra_cmd_options": self.cmd_options} | |
66 | |
67 if make_aux: | |
68 if not self.create_aux_files(ref): | |
69 logging.warn("Auxiliary files were not created.") | |
70 return False | |
71 | |
72 with tempfile.NamedTemporaryFile(suffix=".pileup") as tmp: | |
73 opts["pileup_file"] = tmp.name | |
74 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 | |
75 print cmd | |
76 self.last_command = cmd | |
77 if os.system(cmd) != 0: | |
78 logging.warn("Pileup creation was not successful.") | |
79 return False | |
80 | |
81 return True | |
82 | |
83 def create_aux_files(self, ref): | |
84 """Index reference with faidx from samtools. | |
85 | |
86 Parameters: | |
87 ----------- | |
88 ref: str | |
89 Path to the reference file. | |
90 | |
91 Returns: | |
92 -------- | |
93 bool: | |
94 True if auxiliary files were created, False otherwise. | |
95 """ | |
96 | |
97 success = os.system("samtools faidx %s" % ref) | |
98 | |
99 if success != 0: | |
100 logging.warn("Fasta index could not be created.") | |
101 return False | |
102 | |
103 return True |