Commit a36c0e26 authored by Fabrice Allain's avatar Fabrice Allain
Browse files

Bug revision: filter.out merged when there is more than one contact file...

parent 6352972f
......@@ -21,6 +21,16 @@ from io import StringIO
LOG = logging.getLogger()
def addtup(tup, inc=1):
"""
Increment all values by 1 in a tuple
:param tup:
:param inc:
:return:
"""
return tuple((val + inc for val in tup))
def titleprint(outfile, progname='', desc=''):
"""
Init log file
......
No preview for this file type
......@@ -96,6 +96,8 @@ class AriaEcCommand(object):
default=None, help="configuration file")
parser.add_argument("--nolog", action="store_true",
default=False, help="Don't generate log files")
parser.add_argument("--debug", default=False, action='store_true',
help="Toggle on debug mode")
# Create subcommands
self._create_subparsers(parser.add_subparsers(dest="command"))
return parser
......@@ -228,6 +230,9 @@ class AriaEcCommand(object):
if self.args.output_directory:
LOG.info("Updating output directory")
settings.infra = self.args.output_directory
if self.args.debug:
LOG.info("Toggle on debug mode")
logging.getLogger().setLevel(logging.DEBUG)
return settings
def run(self):
......
No preview for this file type
......@@ -14,6 +14,7 @@ import os
import sys
import json
import re
import textwrap
import aria.legacy.AminoAcid as AminoAcid
from .base import get_filename
from .protein import Protein
......@@ -602,7 +603,11 @@ class AriaEcXMLConverter(AriaXMLConverter):
nb_c)
# Initial contact list start at 0
# pair_list = [(int(x[0]) + 1, int(x[1]) + 1) for x in pair_list]
LOG.info("Selecting %d contacts:\n%s", nb_c, pair_list)
hum_list = [(con1 + 1, con2 +1) for con1, con2 in pair_list]
LOG.info("Selecting %d contacts\n%s", nb_c,
'\n'.join(textwrap.wrap(' '.join(["(%2d, %2d)" % pair for
pair in hum_list]),
width=80)))
if self.settings.setup.config['evfold_weight'] and scoremap is not None:
weight_list = list(float(10.0 / (x + 1)) for x, v in enumerate(
......@@ -648,7 +653,8 @@ class AriaEcXMLConverter(AriaXMLConverter):
contrib_id = 0
for contactidx, contact in enumerate(pair_list):
LOG.debug("Contact %s" % str(contact))
LOG.debug("Writing restraint related to contact %s" %
str((contact[0] + 1, contact[1] + 1)))
# Add neighbors if neigh_flag
resx_idx = range(min_ind(contact[0] - 1),
......
No preview for this file type
No preview for this file type
No preview for this file type
......@@ -25,7 +25,7 @@ from matplotlib import pyplot as plt
import aria.ConversionTable as ConversionTable
import aria.legacy.AminoAcid as AminoAcid
from matplotlib.lines import Line2D
from .base import (tickmin, tickrot, titleprint)
from .base import (tickmin, tickrot, titleprint, addtup)
import sklearn.metrics as skm
LOG = logging.getLogger(__name__)
......@@ -1192,44 +1192,68 @@ class MapFilter(object):
if mapdict.get("scoremap") is not None:
scoremap = mapdict.get("scoremap")
# List of all cysteine
cys_list = [idx for idx, val in enumerate(scoremap.index.values)
if val.endswith('CYS')]
dis_bridge_list = {ss: scoremap.iat[ss] for ss
in list(itertools.combinations(cys_list, 2))}
if dis_bridge_list:
# Get list of all disulfure bridges sorted
sorted_ss = zip(*sorted(dis_bridge_list.items(),
key=operator.itemgetter(1),
reverse=True))[0]
else:
return {'clash': None, 'desc': None}
for dis_bridge in sorted_ss:
# foreach disulfure bridge
LOG.debug('Checking %s %s', addtup(dis_bridge),
scoremap.iat[dis_bridge])
if dis_bridge in clash_list:
# given contact already removed with previous filters
continue
else:
# Check for each cys in dis_bridge if they aready exists
# in unidisbridge_list
exdis = next((unidis for cys in dis_bridge for unidis in
unidisbridge_list if cys in unidis))
if exdis:
if scoremap.iat[dis_bridge] > scoremap.iat[exdis]:
# Better cys--cys contact
# List cys-cys contacts that will be removed in
# unidisbridge_list
remcys = (unidis for cys in dis_bridge for unidis
in unidisbridge_list if cys in unidis)
for dis in remcys:
# PB si un des dis supprime est
clashdisbridge_list.append(dis)
clashdisbridge_list.append(dis[::-1])
unidisbridge_list.remove(dis)
unidisbridge_list.append(dis_bridge)
else:
clashdisbridge_list.append(dis_bridge)
clashdisbridge_list.append(dis_bridge[::-1])
else:
# New cys-cys contact
unidisbridge_list.append(dis_bridge)
exdis = None
if unidisbridge_list:
for cys in dis_bridge:
for unidis in unidisbridge_list:
if cys in unidis:
exdis = unidis
break
if exdis:
LOG.debug("Ss bridge %s have at least one cys in "
"common with previous ss filtered (%s). ",
addtup(dis_bridge),
map(addtup, unidisbridge_list))
if scoremap.iat[dis_bridge] > scoremap.iat[exdis]:
# Better cys--cys contact
# List cys-cys contacts that will be removed in
# unidisbridge_list
LOG.debug("%s has better score than %s",
addtup(dis_bridge),
addtup(exdis))
remcys = (unidis for cys in dis_bridge for unidis
in unidisbridge_list if cys in unidis)
for dis in remcys:
# PB si un des dis supprime est
clashdisbridge_list.append(dis)
clashdisbridge_list.append(dis[::-1])
unidisbridge_list.remove(dis)
LOG.debug("Added %s to clashlist %s",
addtup(dis),
map(addtup, clashdisbridge_list))
unidisbridge_list.append(dis_bridge)
continue
else:
clashdisbridge_list.append(dis_bridge)
clashdisbridge_list.append(dis_bridge[::-1])
LOG.debug("Added %s to ss clashlist %s",
addtup(dis_bridge),
map(addtup, clashdisbridge_list))
continue
LOG.debug('Adding bridge %s to ss filtered list', dis_bridge)
unidisbridge_list.append(dis_bridge)
return {'clash': clashdisbridge_list, 'desc': desc}
else:
# If no score given, return empty list
......@@ -1481,14 +1505,17 @@ class MapFilter(object):
hum_list = [(x + 1, y + 1) for x, y in flt_res.get("clash")]
LOG.info(
"Removed %d contacts:\n%s", len(flt_res.get("clash")) / 2,
hum_list)
'\n'.join(
textwrap.wrap(' '.join(["(%2d, %2d)" % pair for
pair in hum_list]),
width=80)))
clash_dict[flt] = flt_res.get("clash")
if flt_res.get("desc"):
desc_dict.update(flt_res.get("desc"))
# write filter.out
self.write_filtout(clash_dict, desc_dict, contactlist, outdir=outdir,
outprefix=outprefix, clashlist=clashlist)
outprefix=outprefix + '.' + mtype, clashlist=clashlist)
# Contactmap always filtered
# TODO: could set a treshold instead of n_factor
......
No preview for this file type
......@@ -65,43 +65,7 @@ class MapFile(RegexFile):
"regex": re.compile(r"^(?P<res1_nb>\d+)\s+(?P<res1_name>\w)\s+"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>--?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>--?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>-?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>-?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>--?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>--?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>-?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>-?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>--?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>--?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>-?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>-?\d+\.\d+)\s*$"
r"(?P<res2_nb>\d+)\s+(?P<res2_name>\w)\s+"
r"(?P<mi_score>\d)\s+"
r"(?P<plm_score>-?\d+\.\d+)\s*$"),
r"(?P<plm_score>--?\d+\.\d+)\s*$"),
"score_field": "plm_score"
},
"evfold": {
......
No preview for this file type
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment