Mercurial > repos > ulfschaefer > vcfs2fasta

comparison phe/variant/MPileupVariantCaller.py @ 16:1d0bc21232ec draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author ulfschaefer
date Wed, 16 Dec 2015 07:32:22 -0500
parents f72039c5faa4
children
comparison
equal deleted inserted replaced
15:fd42dbe82c12 16:1d0bc21232ec
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