Mercurial > repos > ulfschaefer > vcfs2fasta

comparison phe/variant/MPileupVariantCaller.py @ 22:96f393ad7fc6 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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