|
10
|
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
|
