# -*- coding: utf-8 -*-
################################################################################
# MacSyFinder - Detection of macromolecular systems in protein datasets #
# using systems modelling and similarity search. #
# Authors: Sophie Abby, Bertrand Néron #
# Copyright © 2014 Institut Pasteur (Paris) and CNRS. #
# See the COPYRIGHT file for details #
# #
# MacsyFinder is distributed under the terms of the GNU General Public License #
# (GPLv3). See the COPYING file for details. #
################################################################################
from itertools import groupby
from collections import namedtuple
from glob import glob
import os.path
import logging
_log = logging.getLogger('macsyfinder.' + __name__)
from subprocess import Popen
from macsypy_error import MacsypyError
[docs]def fasta_iter(fasta_file):
"""
:param fasta_file: the file containing all input sequences in fasta format.
:type fasta_file: file object
:author: http://biostar.stackexchange.com/users/36/brentp
:return: for a given fasta file, it returns an iterator which yields tuples
(string id, string comment, int sequence length)
:rtype: iterator
"""
# ditch the boolean (x[0]) and just keep the header or sequence since
# we know they alternate.
faiter = (x[1] for x in groupby(fasta_file , lambda line: line[0] == ">"))
for header in faiter:
# drop the ">"
header = header.next()[1:].strip()
header = header.split()
_id = header[0]
comment = ' '.join(header[1:])
seq = ''.join(s.strip() for s in faiter.next())
length = len(seq)
yield _id, comment, length
[docs]class Indexes(object):
"""
Handle the indexes for macsyfinder:
- find the indexes for hmmer, or build them using formatdb or makeblastdb external tools
- find the indexes required by macsyfinder to compute some scores, or build them.
"""
[docs] def __init__(self, cfg):
"""
The constructor retrieves the file of indexes in the case they are not present
or the user asked for build indexes (--idx)
Launch the indexes building.
:param cfg: the configuration
:type cfg: :class:`macsypy.config.Config` object
"""
self.cfg = cfg
self._fasta_path = cfg.sequence_db
self.name = os.path.basename(cfg.sequence_db)
self._hmmer_indexes = None #list of path
self._my_indexes = None #path
[docs] def build(self, force = False):
"""
Build the indexes from the sequence dataset in fasta format
:param force: If True, force the index building even if the index files are present in the sequence dataset folder
:type force: boolean
"""
hmmer_indexes = self.find_hmmer_indexes()
my_indexes = self.find_my_indexes()
###########################
# build indexes if needed #
###########################
index_dir = os.path.abspath(os.path.dirname(self.cfg.sequence_db))
if force or not hmmer_indexes or not my_indexes:
#formatdb create indexes in the same directory as the sequence_db
#so it must be writable
#if the directory is not writable, formatdb do a Segmentation fault
if not os.access(index_dir, os.R_OK|os.W_OK):
msg = "cannot build indexes, (%s) is not writable" % index_dir
_log.critical(msg)
raise IOError(msg)
if force or not hmmer_indexes:
#self._build_hmmer_indexes() is asynchron
hmmer_indexes_proc = self._build_hmmer_indexes()
if force or not my_indexes:
#self._build_my_indexes() is synchron
self._build_my_indexes()
#################################
# synchronization point between #
# hmmer_indexes and my_indexes #
#################################
if force or not hmmer_indexes:
hmmer_indexes_proc.wait()
if hmmer_indexes_proc.returncode == 127:
msg = "neither makeblastdb nor formatdb can be found, check your config or install makeblastb"
_log.critical( msg, exc_info = True )
raise RuntimeError(msg)
if hmmer_indexes_proc.returncode != 0:
msg = "an error occurred during databases indexation see formatdb.log"
_log.critical( msg, exc_info = True )
raise RuntimeError(msg)
self._hmmer_indexes = self.find_hmmer_indexes()
self._my_indexes = self.find_my_indexes()
assert self._hmmer_indexes , "failed to create hmmer indexes"
assert self._my_indexes, "failed create macsyfinder indexes"
[docs] def find_hmmer_indexes(self):
"""
:return: The hmmer index files.
If indexes are inconsistent (some file(s) missing), a Runtime Error is raised
:rtype: list of string
"""
suffixes = ('.phr', '.pin', '.psd', '.psi', '.psq', '.pal')
idx = []
file_nb = 0
for suffix in suffixes:
index_files = glob( "%s*%s" % (self._fasta_path, suffix))
nb_of_index = len(index_files)
if suffix != '.pal':
if file_nb and file_nb != nb_of_index:
msg = "some index files are missing. Delete all index files (*.phr, *.pin, *.psd, *.psi, *.psq, *.pal) and try to rebuild them."
_log.critical(msg)
raise RuntimeError(msg)
else:
if nb_of_index > 1:
msg = "too many .pal file . Delete all index files (*.phr, *.pin, *.psd, *.psi, *.psq, *.pal) and try to rebuild them."
_log.critical(msg)
raise RuntimeError(msg)
elif file_nb > 1 and nb_of_index == 0:
msg = "some index files are missing. Delete all index files (*.phr, *.pin, *.psd, *.psi, *.psq, *.pal) and try to rebuild them."
_log.critical(msg)
raise RuntimeError(msg)
elif file_nb == 1 and nb_of_index == 1:
msg = "a virtual index is detected (.pal) but there is only one file per index type. Delete all index files (*.phr, *.pin, *.psd, *.psi, *.psq, *.pal) and try to rebuild them."
_log.critical(msg)
raise RuntimeError(msg)
idx.extend(index_files)
file_nb = nb_of_index
return idx
[docs] def find_my_indexes(self):
"""
:return: the file of macsyfinder indexes if it exists in the dataset folder, None otherwise.
:rtype: string
"""
path = os.path.join( os.path.dirname(self.cfg.sequence_db), self.name + ".idx")
if os.path.exists(path):
return path
[docs] def _build_hmmer_indexes(self):
"""
build the index files for hmmer using the formatdb or makeblastdb tool
"""
index_dir = os.path.dirname(self.cfg.sequence_db)
if self.cfg.index_db_exe.find('makeblast') != -1:
command = "%s -title %s -in %s -dbtype prot -parse_seqids" % (self.cfg.index_db_exe,
self.name,
self.cfg.sequence_db)
elif self.cfg.index_db_exe.find('formatdb') != -1:
# -t Title for database file [String]
# -i Input file(s) for formatting [File In]
# -p T Type of file = protein
# -o T Parse SeqId and create indexes.
# -s T Create indexes limited only to accessions
command = "%s -t %s -i %s -p T -o T -s T" % ( self.cfg.index_db_exe,
self.name,
self.cfg.sequence_db
)
else:
raise MacsypyError("%s is not supported to index the sequence dataset. Please use makeblastdb or formatdb." % self.cfg.sequence_db)
_log.debug("hmmer index command: {0}".format(command))
err_path = os.path.join(index_dir, "formatdb.err")
with open(err_path, 'w') as err_file:
try:
formatdb = Popen( command ,
shell = True ,
stdout = err_file ,
stdin = None ,
stderr = err_file ,
close_fds = False ,
)
except Exception, err:
msg = "unable to index the sequence dataset : %s : %s" % (command, err)
_log.critical( msg, exc_info = True )
raise err
return formatdb
[docs] def _build_my_indexes(self):
"""
Build macsyfinder indexes. These indexes are stored in a file.
The file format is the following:
- one entry per line, with each line having this format:
- sequence id;sequence length;sequence rank
"""
try:
with open(self._fasta_path, 'r') as fasta_file:
with open(os.path.join(os.path.dirname(self.cfg.sequence_db), self.name + ".idx" ), 'w') as my_base:
f_iter = fasta_iter(fasta_file)
seq_nb = 0
for seqid, comment, length in f_iter:
seq_nb += 1
my_base.write("%s;%d;%d\n" % (seqid, length, seq_nb))
except Exception, err:
msg = "unable to index the sequence dataset: %s : %s" % (self.cfg.sequence_db, err)
_log.critical(msg, exc_info = True)
raise err
"""handle name, topology type, and min/max positions in the sequence dataset for a replicon"""
RepliconInfo = namedtuple('RepliconInfo', 'topology, min, max, genes')
[docs]class RepliconDB(object):
"""
Stores information (topology, min, max, [genes]) for all replicons in the sequence_db
the Replicon object must be instantiated only for sequence_db of type 'gembase' or 'ordered_replicon'
"""
ordered_replicon_name = 'UserReplicon'
[docs] def __init__(self, cfg):
"""
:param cfg: The configuration object
:type cfg: :class:`macsypy.config.Config` object
.. note ::
This class can be instanciated only if the db_type is 'gembase' or 'ordered_replicon'
"""
self.cfg = cfg
assert self.cfg.db_type in ('gembase', 'ordered_replicon')
idx = Indexes(self.cfg)
self.sequence_idx = idx.find_my_indexes()
self.topology_file = self.cfg.topology_file
self._DB = {}
if self.topology_file:
topo_dict = self._fill_topology()
else:
topo_dict = {}
if self.cfg.db_type == 'gembase':
self._fill_gembase_min_max(topo_dict, default_topology = self.cfg.replicon_topology)
else:
self._fill_ordered_min_max(self.cfg.replicon_topology)
[docs] def _fill_topology(self):
"""
Fill the internal dictionary with min and max positions for each replicon_name of the sequence_db
"""
topo_dict = {}
with open(self.topology_file) as topo_f:
for l in topo_f:
if l.startswith('#'):
continue
replicon_name, topo = l.split(':')
replicon_name = replicon_name.strip()
topo = topo.strip().lower()
topo_dict[replicon_name] = topo
return topo_dict
[docs] def _fill_ordered_min_max(self, default_topology = None):
"""
For the replicon_name of the ordered_replicon sequence base, fill the internal dict with RepliconInfo
:param default_topology: the topology provided by config.replicon_topology
:type default_topology: string
"""
_min = 1
#self.sequence_idx is a file with the following structure seq_id;seq_length;seq_rank\n
with open(self.sequence_idx) as idx_f:
_max = 0
genes = []
for l in idx_f:
seq_id, length, _rank = l.split(";")
genes.append((seq_id, length))
_max += 1
self._DB[self.ordered_replicon_name] = RepliconInfo(default_topology, _min, _max, genes)
[docs] def _fill_gembase_min_max(self, topology, default_topology):
"""
For each replicon_name of a gembase dataset, it fills the internal dictionary with a namedtuple RepliconInfo
:param topology: the topologies for each replicon
(parsed from the file specified with the option --topology-file)
:type topology: dict
:param default_topology: the topology provided by the config.replicon_topology
:type default_topology: string
"""
def grp_replicon(line):
"""
in gembase the idtentifier of fasta sequence follow the following schema: replicon_name_seq_name
so grp_replicon allow to group sequences belonging to the same replicon.
"""
return line.split('_')[0]
def parse_entry(entry):
"""
parse an entry in the index file (.idx)
an entry have the folowing format sequence_id;sequence lenght;sequence rank in replicon
"""
entry = entry.rstrip()
seq_id, length, rank = entry.split(';')
replicon_name , seq_name = seq_id.split('_')
return replicon_name, seq_name, length, int(rank)
with open(self.sequence_idx) as idx_f:
replicons = (x[1] for x in groupby(idx_f, grp_replicon))
for replicon in replicons:
genes = []
entry = replicon.next()
replicon_name, seq_name, seq_lenght, _min = parse_entry(entry)
genes.append((seq_name, seq_lenght))
for entry in replicon:
#pass all sequence of the replicon until the last one
_, seq_name, seq_lenght, _ = parse_entry(entry)
genes.append((seq_name, seq_lenght))
_, seq_name, seq_lenght, _max = parse_entry(entry)
genes.append((seq_name, seq_lenght))
if replicon_name in topology:
self._DB[replicon_name] = RepliconInfo(topology[replicon_name], _min, _max, genes)
else:
self._DB[replicon_name] = RepliconInfo(default_topology, _min, _max, genes)
[docs] def __contains__(self, replicon_name):
"""
:param replicon_name: the name of the replicon
:type replicon_name: string
:returns: True if replicon_name is in the repliconDB, false otherwise.
:rtype: boolean
"""
return replicon_name in self._DB
[docs] def __getitem__(self, replicon_name):
"""
:param replicon_name: the name of the replicon to get information on
:type replicon_name: string
:returns: the RepliconInfo for the provided replicon_name
:rtype: :class:`RepliconInfo` object
:raise: KeyError if replicon_name is not in repliconDB
"""
return self._DB[replicon_name]
[docs] def get(self, replicon_name, default = None):
"""
:param replicon_name: the name of the replicon to get informations
:type replicon_name: string
:param default: the value to return if the replicon_name is not in the RepliconDB
:type default: any
:returns: the RepliconInfo for replicon_name if replicon_name is in the repliconDB, else default. If default is not given, it is set to None, so that this method never raises a KeyError.
:rtype: :class:`RepliconInfo` object
"""
return self._DB.get(replicon_name, default)
[docs] def items(self):
"""
:return: a copy of the RepliconDB as a list of (replicon_name, RepliconInfo) pairs
"""
return self._DB.items()
[docs] def iteritems(self):
"""
:return: an iterator over the RepliconDB as a list (replicon_name, RepliconInfo) pairs
"""
return self._DB.iteritems()
[docs] def replicon_names(self):
"""
:return: a copy of the RepliconDB as a list of replicon_names
"""
return self._DB.keys()
[docs] def replicon_infos(self):
"""
:return: a copy of the RepliconDB as list of replicons info
:rtype: RepliconInfo instance
"""
return self._DB.values()