Commit 50b28391 authored by fabrice's avatar fabrice

Bug fix: sc_sc contactdef ignored

parent ac092778
No preview for this file type
No preview for this file type
......@@ -72,17 +72,25 @@ class AriaEcContactMap(object):
nb_c = int(len(self.protein.aa_sequence.sequence) *
self.settings.setup.config.get("n_factor"))
refmap = self.refmap["contactmap"]
plotparams = {k: self.settings.contactmap.config.get(k, None)
for k in ('size_fig', 'plot_ext', 'plot_dpi')}
for mapt in self.allresmap.keys():
outdir = self.settings.outdir
prefix = "_".join((self.outprefix, mapt, self.reftype))
if mapt == self.reftype:
refmap.write_contacts("_".join((self.outprefix, mapt)),
outdir=outdir)
refmap.saveplot(outdir=outdir,
outprefix="_".join((self.outprefix, mapt)),
**plotparams)
continue
refmap = self.refmap["contactmap"]
if self.allresmap[mapt].get("contactmap") is not None and \
self.allresmap[mapt].get("scoremap") is not None:
prefix = "_".join((self.outprefix, mapt, self.reftype))
outdir = self.settings.outdir
if not self.settings.contactmap.args.get("nofilter"):
# Get filtered contact map/list
cmpmap = self.allresmap[mapt]["contactmap"].topmap(
self.allresmap[mapt]["scoremap"], nb_c)
cmplist = self.allresmap[mapt]['scoremap'].sortedset(
......@@ -101,9 +109,11 @@ class AriaEcContactMap(object):
# TODO: elementwise error with compare method
refmap.compare(cmpmap, outprefix=prefix,
out_dir=outdir,
**{k: self.settings.contactmap.config.get(k, None)
for k in ('save_fig', 'size_fig',
'plot_ext', 'plot_dpi', 'alpha')})
save_fig=self.settings.contactmap.config.get(
"save_fig"),
alpha=self.settings.contactmap.config.get(
"alpha"),
**plotparams)
logger.info(pd.crosstab(cmpmap.to_series(), refmap.to_series(),
rownames=[mapt],
colnames=[self.reftype]))
......@@ -123,9 +133,11 @@ class AriaEcContactMap(object):
# TODO: elementwise error with compare method
refmap.compare(cmpmap, outprefix=prefix,
out_dir=outdir,
**{k: self.settings.contactmap.config.get(k, None)
for k in ('save_fig', 'size_fig',
'plot_ext', 'plot_dpi', 'alpha')})
save_fig=self.settings.contactmap.config.get(
"save_fig"),
alpha=self.settings.contactmap.config.get(
"alpha"),
**plotparams)
logger.info(pd.crosstab(cmpmap.to_series(), refmap.to_series(),
rownames=[mapt], colnames=[self.reftype]))
else:
......
No preview for this file type
No preview for this file type
......@@ -168,6 +168,7 @@ class ProteinMap(Map):
kwargs["columns"] = col
super(ProteinMap, self).__init__(**kwargs)
self.contact_flags = flaglist if flaglist else None
self._maplot = None
def _constructor_expanddim(self):
super(ProteinMap, self)._constructor_expanddim()
......@@ -176,6 +177,34 @@ class ProteinMap(Map):
def sequence(self):
raise NotImplementedError
@property
def maplot(self):
# Contact map Plot
if not self._maplot:
minticks = tickmin(self, shift=1) # Nb graduations
self._maplot = sns.heatmap(self, square=True, cbar=False,
linewidths=1, vmax=1, vmin=-1,
cmap=sns.diverging_palette(20, 220, n=7,
as_cmap=True),
xticklabels=minticks[0],
yticklabels=minticks[1])
return self._maplot
def saveplot(self, outdir='', outprefix="protein", size_fig=10,
plot_ext="pdf", plot_dpi=200):
f, ax = plt.subplots(figsize=(12, 9))
tickrot(self.maplot.axes, self.maplot.figure,
rotype='horizontal')
self.maplot.figure.set_size_inches(size_fig, size_fig)
map_title = "%s contacts map" % outprefix
self.maplot.set_title(map_title)
self.maplot.figure.tight_layout()
f.tight_layout()
self.maplot.figure.savefig(os.path.join(outdir, "%s.contactmap.%s" % (
outprefix, plot_ext)), dpi=plot_dpi)
def contactset(self, human_idx=False):
# Remove duplicate in contact_list
if self.contact_list():
......@@ -215,8 +244,7 @@ class ProteinMap(Map):
# Indexation matrice (tous les atomes ou tous les residus)
raise NotImplementedError
def compare(self, protmap, save_fig=True, out_dir='', outprefix="protein",
size_fig=10, plot_ext="pdf", plot_dpi=200, alpha=None):
def compare(self, protmap, save_fig=True, alpha=None, **kwargs):
# Contact map plot
if getattr(protmap, "shape") and self.shape != protmap.shape:
logging.error("Cant't compare %s map with %s" % (
......@@ -224,14 +252,14 @@ class ProteinMap(Map):
return None
else:
# Contact map Plot
minticks = tickmin(self, shift=1) # Nb graduations
f, ax = plt.subplots(figsize=(12, 9))
contactmaplot = sns.heatmap(self, square=True, cbar=False,
linewidths=1, vmax=1, vmin=-1,
cmap=sns.diverging_palette(20, 220, n=7,
as_cmap=True),
xticklabels=minticks[0],
yticklabels=minticks[1])
# minticks = tickmin(self, shift=1) # Nb graduations
# f, ax = plt.subplots(figsize=(12, 9))
# contactmaplot = sns.heatmap(self, square=True, cbar=False,
# linewidths=1, vmax=1, vmin=-1,
# cmap=sns.diverging_palette(20, 220, n=7,
# as_cmap=True),
# xticklabels=minticks[0],
# yticklabels=minticks[1])
cmplist = protmap.contact_list()
ymax = len(self.sequence) - 1
......@@ -250,9 +278,9 @@ class ProteinMap(Map):
color = pal[i]
mark = Line2D.filled_markers[i]
for x, y in zip(xind, yind):
contactmaplot.axes.scatter(x, y, s=20, c=color,
linewidths=0,
alpha=alpha, marker=mark)
self.maplot.axes.scatter(x, y, s=20, c=color,
linewidths=0, alpha=alpha,
marker=mark)
else:
xind = [x + .5 for x in
zip(*cmplist)[0] + zip(*cmplist)[1]]
......@@ -262,27 +290,24 @@ class ProteinMap(Map):
# width = [0.3 for _ in xind]
# for x, y, h in zip(xind, yind, width):
for x, y in zip(xind, yind):
contactmaplot.axes.scatter(x, y, s=20, c=color,
linewidths=0,
alpha=alpha)
# contactmaplot.axes.add_artist(ptc.Circle(xy=(x, y), alpha=alpha,
# color="green",
# radius=h))
self.maplot.axes.scatter(x, y, s=20, c=color, linewidths=0,
alpha=alpha)
if save_fig:
tickrot(contactmaplot.axes, contactmaplot.figure,
rotype='horizontal')
contactmaplot.figure.set_size_inches(size_fig,
size_fig)
map_title = "%s contacts map" % outprefix
contactmaplot.set_title(map_title)
contactmaplot.figure.tight_layout()
f.tight_layout()
contactmaplot.figure.savefig(os.path.join(out_dir,
"%s.contactmap.%s" % (
outprefix,
plot_ext)),
dpi=plot_dpi)
self.saveplot(**kwargs)
# tickrot(contactmaplot.axes, contactmaplot.figure,
# rotype='horizontal')
# contactmaplot.figure.set_size_inches(size_fig,
# size_fig)
# map_title = "%s contacts map" % outprefix
# contactmaplot.set_title(map_title)
# contactmaplot.figure.tight_layout()
#
# f.tight_layout()
# contactmaplot.figure.savefig(os.path.join(out_dir,
# "%s.contactmap.%s" % (
# outprefix,
# plot_ext)),
# dpi=plot_dpi)
def compare_contactmap(self, cmpmap, contactlist, outprefix, outdir="",
distmap=None,
......@@ -438,6 +463,7 @@ class ResAtmMap(ProteinMap):
if self.dtype == bool:
# If self is already a contact map
return None
# TODO: issue with sc_sc treshold !!!!!
logger.info("Generate contact map using contact definition %s" %
contactdef)
# Contact map initialize to a boolean matrix filled with False
......@@ -457,12 +483,14 @@ class ResAtmMap(ProteinMap):
atm_list = set(self.index.get_level_values(1))
atms_list = set([(a, b) for a in atm_list for b in atm_list])
for pair in contactdef.keys():
if pair == "default_cutoff":
continue
treshold = contactdef[pair]
pair = tuple(pair.upper().split("_"))
if pair not in atms_list:
# Already applied a treshold for this pair
continue
elif pair == ("SC", "SC"):
logger.info(
"Filtering values in matrix related to %s (%s)" %
(str(pair), str(treshold)))
if pair in (("SC", "SC"), ("sc", "sc")):
# Use scsc_min to apply treshold only for selected atom
# sidechain
idx_list = []
......@@ -480,13 +508,14 @@ class ResAtmMap(ProteinMap):
pair[1]) and x[1] == atm2))
mask = ([any(tup) for tup in zip(*idx_list)],
[any(tup) for tup in zip(*col_list)])
elif pair not in atms_list:
logger.error("Pair %s doesn't exist ..." % str(pair))
# Already applied a treshold for this pair
continue
else:
logger.debug("Apply treshold for %s" % str(pair))
atms_list.discard(pair)
# Selecting rows for each atom
logger.info(
"Filtering values in matrix related to %s (%s)" %
(str(pair), str(treshold)))
mask = (self.index.get_level_values(1) == pair[0],
self.index.get_level_values(1) == pair[1])
tmp = self.loc[mask].apply(lambda x: x < float(treshold))
......@@ -946,7 +975,7 @@ class MapFilter:
if mapfilters == "all":
mapfilters = self.filter_types
elif not mapfilters:
mapfilters = ["pos", ]
mapfilters = ["pos"]
else:
if type(mapfilters) == list:
mapfilters = [elm for elm in mapfilters if elm in self.filter_types]
......
No preview for this file type
No preview for this file type
Markdown is supported
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