diff --git a/ariaec/econverter.pyc b/ariaec/econverter.pyc index 81fa84027c6739ce54a70a6ad347b328da8c30a4..e3481cf613ef97b7a6d67d480cabccf0baa8baa9 100644 Binary files a/ariaec/econverter.pyc and b/ariaec/econverter.pyc differ diff --git a/ariaec/ecsettings.pyc b/ariaec/ecsettings.pyc index 19f73a296b5c13e0ec0d16e4b6f932549dfd827b..b44fb90f2a6510c1985e58448fbbea72f8fd679c 100644 Binary files a/ariaec/ecsettings.pyc and b/ariaec/ecsettings.pyc differ diff --git a/ariaec/maplot.py b/ariaec/maplot.py index f5800cf9f92f6773b9475d2acddeca0ef7f34b59..2be61d14ded7358034d234cc43a882f70412052f 100644 --- a/ariaec/maplot.py +++ b/ariaec/maplot.py @@ -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: diff --git a/ariaec/maplot.pyc b/ariaec/maplot.pyc index 8fd4173edea62cd64849117317a17db95abb7157..9309021709176257fae3c5271192e2a9f44936cd 100644 Binary files a/ariaec/maplot.pyc and b/ariaec/maplot.pyc differ diff --git a/ariaec/protein.pyc b/ariaec/protein.pyc index 403f62a32dae04748fa0616abf20cdc067c654cc..5d811cfb9d7a9306750c04658f30e0269554a469 100644 Binary files a/ariaec/protein.pyc and b/ariaec/protein.pyc differ diff --git a/ariaec/protmap.py b/ariaec/protmap.py index 719dad83cd8a9fbfe9ee9f24891dfe219e18b12a..d5c475286c36d3885d5a5eab9523244e4711bc0e 100644 --- a/ariaec/protmap.py +++ b/ariaec/protmap.py @@ -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] diff --git a/ariaec/protmap.pyc b/ariaec/protmap.pyc index 8183bad2755bd1f3ffeb4920ab8b9412ecfc17e4..a7d3811e93dd82a249644f868314f5b3e3353fe5 100644 Binary files a/ariaec/protmap.pyc and b/ariaec/protmap.pyc differ diff --git a/ariaec/reader.pyc b/ariaec/reader.pyc index 3a3e015d302ea5a636aa5ec9cbe225a0659cbaea..5eeb4d61ab4ba9568ce24f11b49683d31a4434f0 100644 Binary files a/ariaec/reader.pyc and b/ariaec/reader.pyc differ