maplot.py 9.74 KB
Newer Older
1
# coding=utf-8
fabrice's avatar
fabrice committed
2
3
4
5
6
"""
                            Input/Output aria_ec scripts
"""
from __future__ import absolute_import, division, print_function

7
import sys
fabrice's avatar
fabrice committed
8
9
10
import json
import logging
from .base import get_filename
11
from .reader import MapFileListReader
fabrice's avatar
fabrice committed
12
13
14
from .protmap import MapFilter
from .protein import Protein

15
LOG = logging.getLogger(__name__)
fabrice's avatar
fabrice committed
16
17
18


class AriaEcContactMap(object):
19
20
21
22
    """
    Contact maplot class
    """

fabrice's avatar
fabrice committed
23
24
25
26
    def __init__(self, settings):
        # TODO: check_type settings (AriaEcSettings)
        self.settings = settings
        self.protein = Protein(settings)
27
        self.file_reader = MapFileListReader(
fabrice's avatar
fabrice committed
28
29
30
31
32
            cont_def=settings.contactdef.config)
        self.filter = MapFilter(settings.setup.config)
        self.outprefix = ''
        self.allresmap = {}
        self.refmap = None
33
        self.refname = None
fabrice's avatar
fabrice committed
34
35
36
37
        self.reftype = ''

    def run(self):
        # Check input
38
39
40
        """
        Main method
        """
41
42
43
44
        LOG.debug("Settings:\n" + json.dumps(self.settings.contactmap.config,
                                             indent=4))
        LOG.debug("Args:\n" + json.dumps(self.settings.contactmap.args,
                                         indent=4))
fabrice's avatar
fabrice committed
45
46
        if not self.settings.contactmap.args.get("onlyreport", False):
            self.settings.make_infra()
fabrice's avatar
fabrice committed
47
        # ----------------------------- Input -------------------------------- #
48
49
50
51
        if self.settings.contactmap.args.get("prefix") and len(
                self.settings.contactmap.args.get("infiles")) == len(
                    self.settings.contactmap.args.get("prefix")):
                self.outprefix = self.settings.contactmap.args.get("prefix", "")
52
        else:
53
54
            self.outprefix = get_filename(self.settings.contactmap.args.get("seq",
                                                                            None))
fabrice's avatar
fabrice committed
55
56
57
58
59
        # Load Sequence file
        self.protein.set_aa_sequence(self.settings.contactmap.args.get("seq", None))
        # Load secondary structure prediction file
        self.protein.set_sec_struct(self.settings.contactmap.args.get("sspred",
                                                                      None),
60
                                    ssdist_filename=self.settings.ssdist,
fabrice's avatar
fabrice committed
61
62
63
64
65
66
                                    ssidx=self.settings.contactmap.args.get(
                                        "ssidx", False))
        # ---------------------------- Processing ---------------------------- #
        self.file_reader.read(self.settings.contactmap.args.get("infiles"),
                              filetypelist=self.settings.contactmap.args.get("contact_types"),
                              protein=self.protein,
67
                              groupby_method=self.settings.setup.config[
fabrice's avatar
fabrice committed
68
69
70
71
72
73
74
75
                                  'groupby_method'],
                              scsc=self.settings.scsc_min)

        for idx, fo in enumerate(self.file_reader.filelist):
            # fo need a contactmap in order to wite XML dist restraints
            # TODO: filter pour toutes les map de mapdict !! (fonction remove
            #  s'applique sur l'humanidx contenant les residus)
            if idx == 0:
76
                LOG.info("%s map set as reference", fo.filetype.capitalize())
fabrice's avatar
fabrice committed
77
78
                self.refmap = fo.mapdict
                self.reftype = fo.filetype
79
80
                self.refname = fo.filename if type(self.outprefix) != list \
                    else self.outprefix[idx]
81
82
83
84
85
            if (self.settings.contactmap.args.get("filter") and idx != 0) or \
                    (self.settings.contactmap.args.get("filter")
                     and len(self.file_reader.filelist) == 1):
                # Filtering all contactmap except the reference map
                LOG.info("Filtering %s map", fo.filetype)
fabrice's avatar
fabrice committed
86
87
                self.filter(fo.mapdict, fo.filetype, fo.contactlist,
                            self.protein, clashlist=fo.clashlist,
88
89
                            outprefix=self.outprefix[idx] if type(
                                self.outprefix) == list else self.outprefix,
fabrice's avatar
fabrice committed
90
                            outdir=self.settings.outdir)
91
92
93
94
            # else:
            #     Use only position filter
            #     self.filter(fo.mapdict, fo.filetype, fo.contactlist,
            #                 self.protein, clashlist=fo.clashlist,
95
            #                 outprefix=self.outprefix,
96
            #                 outdir=self.settings.outdir, mapfilters="pos")
97
98
99
            self.allresmap[(fo.filename if type(self.outprefix) != list else
                            self.outprefix[idx], fo.filetype,
                            fo.filepath)] = fo.mapdict
fabrice's avatar
fabrice committed
100

101
102
103
        try:
            refmap = self.refmap["contactmap"]
        except TypeError:
104
            LOG.error("First contact map should be a valid file")
105
106
            sys.exit(1)

107
108
        # nb_c = int(len(self.protein.aa_sequence.sequence) *
        #            self.settings.setup.config.get("n_factor"))
fabrice's avatar
fabrice committed
109
110
        plotparams = {k: self.settings.contactmap.config.get(k, None)
                      for k in ('size_fig', 'plot_ext', 'plot_dpi')}
111
112
        outdir = self.settings.outdir

fabrice's avatar
fabrice committed
113
114
        if self.settings.contactmap.config.get("save_fig") and not \
                self.settings.contactmap.args.get("onlyreport", False):
115
            refmap.saveplot(outdir=outdir,
116
117
                            outprefix=self.outprefix[0] if type(
                                self.outprefix) == list else self.outprefix,
118
                            **plotparams)
fabrice's avatar
fabrice committed
119

fabrice's avatar
fabrice committed
120
121
122
123
        if self.settings.contactmap.args.get("merge", None):
            # Combine merge maps with other maps
            mergelist = self.settings.contactmap.args.get("merge")
            for mergetype in mergelist:
124
                if mergetype in zip(self.allresmap.keys())[0]:
fabrice's avatar
fabrice committed
125
126
                    mergemaps = self.allresmap.pop(mergetype)
                    mergecontactmap = mergemaps.get("contactmap")
127
                    for mapname, mapt in self.allresmap.keys():
fabrice's avatar
fabrice committed
128
                        if mapt != self.reftype:
129
                            # TODO: DON'T WORK !!!!
130
131
                            LOG.info("Merging %s with %s map",
                                     mergetype, mapt)
fabrice's avatar
fabrice committed
132
                            up_map = self.allresmap[mapt]["contactmap"]
133
                            up_map[:] = up_map[:] + mergecontactmap[:]
fabrice's avatar
fabrice committed
134
135
136
137
                            mergekey = "%s_%s" % (mapt, mergetype)
                            self.allresmap[mergekey] = {}
                            self.allresmap[mergekey]["contactmap"] = up_map

138
        for mapname, mapt, mapath in self.allresmap.keys():
139

140
            prefix = "_".join((mapname, self.refname)).replace(".", "_")
141
142

            if mapname == self.refname:
143
                if self.settings.contactmap.args.get("onlyreport", False) is not False:
144
145
146
147
                    refmap.write_contacts(mapname,
                                          outdir=outdir,
                                          scoremap=self.refmap.get("scoremap",
                                                                   None))
148
149
                continue

150
151
            scoremap = self.allresmap[(mapname, mapt, mapath)].get(
                'scoremap', None)
152
153
154
155
156
157
158
159
160
            # if self.allresmap[mapt].get("contactmap") is not None and \
            #         self.allresmap[mapt].get("scoremap") is not None:
            #     Get top contact map/list
            #     cmpmap = self.allresmap[mapt]["contactmap"].topmap(
            #         self.allresmap[mapt]["scoremap"], nb_c)
            #     cmplist = self.allresmap[mapt]['scoremap'].sortedset(
            #         human_idx=True)[:nb_c]
            # elif self.allresmap[mapt].get("contactmap") is not None:
            #     If no score given, use all contact list
161
162
163
            cmpmap = self.allresmap[(mapname, mapt, mapath)]["contactmap"]
            cmplist = self.allresmap[(mapname, mapt, mapath)][
                'contactmap'].contact_list(human_idx=True)
164
            # else:
165
166
            #     LOG.warning("%s map doesn't have any score related. Can't "
            #                    "define top list related to this map", mapt)
167
            #     continue
fabrice's avatar
fabrice committed
168

169
            # TODO: only one function for output files
170
            # Write contact list in txt file
171
            refmap.report(cmpmap, scoremap=scoremap,
fabrice's avatar
fabrice committed
172
173
174
                          outdir=outdir,
                          plotag=not self.settings.contactmap.args.get(
                              "onlyreport"),
fabrice's avatar
fabrice committed
175
176
177
                          plotdir=self.settings.infra.get("graphics", outdir)
                          if not self.settings.contactmap.args.get(
                              "onlyreport") else self.settings.outdir)
178
179
            # Contact map comparison plot
            # TODO: elementwise error with compare method
fabrice's avatar
fabrice committed
180
181
            # Write cmp stats
            if not self.settings.contactmap.args.get("onlyreport", False):
182
                cmpmap.write_contacts(mapname,
fabrice's avatar
fabrice committed
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
                                      scoremap=scoremap,
                                      outdir=outdir)
                cmpmap.compare_contactmap(refmap, cmplist, prefix,
                                          distmap=self.refmap["distmap"],
                                          human_idx=True,
                                          outdir=outdir)
                refmap.compareplot(cmpmap, outprefix=prefix,
                                   outdir=outdir,
                                   save_fig=self.settings.contactmap.config.get(
                                       "save_fig"),
                                   alpha=self.settings.contactmap.config.get(
                                       "alpha"),
                                   **plotparams)
                # Contingency table
                # print(cmpmap.to_series())
198
                # LOG.info(pd.crosstab(cmpmap.values, refmap.values,
fabrice's avatar
fabrice committed
199
                #                         rownames=[mapt], colnames=[self.reftype]))