Source code for phenix.viewers.viewer_search_fit

# **************************************************************************
# *
# * Authors:     Roberto Marabini
# *              Marta Martinez (mmmtnez@cnb.csic.es)
# *
# * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
# *
# * This program is free software; you can redistribute it and/or modify
# * it under the terms of the GNU General Public License as published by
# * the Free Software Foundation; either version 2 of the License, or
# * (at your option) any later version.
# *
# * This program is distributed in the hope that it will be useful,
# * but WITHOUT ANY WARRANTY; without even the implied warranty of
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# * GNU General Public License for more details.
# *
# * You should have received a copy of the GNU General Public License
# * along with this program; if not, write to the Free Software
# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
# * 02111-1307  USA
# *
# *  All comments concerning this program package may be sent to the
# *  e-mail address 'scipion@cnb.csic.es'
# *
# **************************************************************************
import math
import os
import sqlite3

import matplotlib.pyplot as plt
from tkinter import messagebox
from phenix.protocols.protocol_search_fit import PhenixProtSearchFit
from pyworkflow.protocol.params import LabelParam, EnumParam, IntParam, FloatParam
from pyworkflow.viewer import DESKTOP_TKINTER, WEB_DJANGO, ProtocolViewer
from pwem.viewers import TableView, Chimera
from pwem import Domain
from phenix.protocols.protocol_search_fit import (DATAFILE,
                                                  TABLE)

[docs]def errorWindow(tkParent, msg): try: # if tkRoot is null the error message may be behind # other windows messagebox.showerror("Error", # bar title msg, # message parent=tkParent) except: print(("Error:", msg))
[docs]class PhenixProtRuSearchFitViewer(ProtocolViewer): """ Viewer for Phenix program DEARCHFIT """ _label = 'Search Fit Viewer' _environments = [DESKTOP_TKINTER, WEB_DJANGO] _targets = [PhenixProtSearchFit] SEARCHFITSUBPLOTSFILENAME = 'model_to_map_fit.py' def __init__(self, **kwargs): ProtocolViewer.__init__(self, **kwargs) def _defineParams(self, form): form.addSection(label="Volume and models") form.addParam('showMapModel', LabelParam, label="Show volume and atomic structures in ChimeraX", help="Display input volume and atom struct" " plus 5 better atom struct candidates.") form.addParam("numAtomStruct", IntParam, label="Max. Number Atom Structs.", default=1000, help="Number of atom structs to show ordered by model_to_map_fit\n") form.addParam("zone", FloatParam, label="Show Area around input atomic struct (A)", default=3, help="Limit the display to a zone around the input atomic structure.\n" "Units = A.") form.addParam('showPlot', LabelParam, label="Summary Plot", help="Plot showing 'model_to_map_fit' values") def _getVisualizeDict(self): return { 'showMapModel': self._showMapModel, 'showPlot': self._showPlot } def _showMapModel(self, e=None): bildFileName = self.protocol._getExtraPath("axis_output.bild") _inputVol = self.protocol.inputVolume.get() dim = _inputVol.getDim()[0] sampling = _inputVol.getSamplingRate() counter = 1 Chimera.createCoordinateAxisFile(dim, bildFileName=bildFileName, sampling=sampling) fnCmd = os.path.abspath(self.protocol._getExtraPath("chimera_output.cxc")) f = open(fnCmd, 'w') # change to workingDir # If we do not use cd and the project name has an space # the protocol fails even if we pass absolute paths f.write('cd %s\n' % os.getcwd()) # reference axis model = 0 f.write("open %s\n" % os.path.abspath(bildFileName)) # 1 f.write("cofr 0,0,0\n") # set center of coordinates # input 3D map counter += 1 # 2 vol = self._getInputVolume() fnVol = os.path.abspath(vol.getFileName()) f.write("open %s\n" % fnVol) x, y, z = vol.getOrigin(force=True).getShifts() sampling = vol.getSamplingRate() f.write("volume #%d style surface voxelSize %f\nvolume #%d origin " "%0.2f,%0.2f,%0.2f\n" % (counter, sampling, counter, x, y, z)) # input PDB (usually from coot) counter += 1 # 3 pdbFileName = os.path.abspath(self.protocol.inputStructure.get().getFileName()) f.write("open %s\n" % pdbFileName) # zone command f.write("volume zone #%d near #%d range %f newMap false\n" % (counter -1, counter, self.zone.get())) # open database and retrieve all files conn = sqlite3.connect(os.path.abspath(self.protocol._getExtraPath(DATAFILE))) c = conn.cursor() sqlCommand = """SELECT filename FROM %s WHERE model_to_map_fit != -1 ORDER BY model_to_map_fit DESC LIMIT %d""" % (TABLE, self.numAtomStruct) c.execute(sqlCommand) rows = c.fetchall() for row in rows: atomStructFn = row[0][:-4] + "_real_space_refined.cif" f.write("open %s\n" % atomStructFn) c.close() conn.close() f.close() # run in the background chimeraPlugin = Domain.importFromPlugin('chimera', 'Plugin', doRaise=True) chimeraPlugin.runChimeraProgram(chimeraPlugin.getProgram(), fnCmd + "&", cwd=os.getcwd()) return [] def _getInputVolume(self): if self.protocol.inputVolume.get() is None: fnVol = self.protocol.inputStructure.get().getVolume() else: fnVol = self.protocol.inputVolume.get() return fnVol def _showPlot(self, e=None): xList = [] yList = [] conn = sqlite3.connect(os.path.abspath(self.protocol._getExtraPath(DATAFILE))) c = conn.cursor() sqlCommand = """ SELECT id, model_to_map_fit FROM ( SELECT id, model_to_map_fit FROM %s WHERE model_to_map_fit != -1 ORDER BY model_to_map_fit desc LIMIT %d) AS a ORDER BY id """ % (TABLE, self.numAtomStruct) c.execute(sqlCommand) rows = c.fetchall() for row in rows: xList.append(float(row[0])) yList.append(float(row[1])) # compute avg sqlCommand = """SELECT AVG(model_to_map_fit) FROM (SELECT model_to_map_fit FROM %s WHERE model_to_map_fit != -1 ORDER BY model_to_map_fit desc )""" % (TABLE) c.execute(sqlCommand) rows = c.fetchone(); avg = rows[0] print("avg", avg) fromRelation = """(SELECT model_to_map_fit FROM %s WHERE model_to_map_fit != -1 ORDER BY model_to_map_fit desc ) AS mainTable""" % (TABLE) sqlCommand = """SELECT AVG((mainTable.model_to_map_fit - sub.a) * (mainTable.model_to_map_fit - sub.a)) as var FROM %s, (SELECT AVG(model_to_map_fit) AS a FROM %s WHERE model_to_map_fit != -1 ) AS sub WHERE model_to_map_fit != -1 """ % (fromRelation, fromRelation) c.execute(sqlCommand) rows = c.fetchone(); std_2 = rows[0] std = math.sqrt(std_2) c.close() conn.close() if not (xList or yList): errorWindow(self.getTkRoot(), "No data available") return title = 'avg (all data) = %f, std (all data) = %f'%(avg, std) plt.plot(xList, yList, 'x') plt.axis([0, max(xList) + 1.0, 0.0, max(yList)+0.1]) plt.title(title) plt.xlabel('#Atom Structs') plt.ylabel('Map Model Fit Score') plt.show() print(xList, yList)
# SELECT FLOOR(model_to_map_fit/5.00)*5 As Grade, # COUNT(*) AS [Grade Count] # FROM TableName # GROUP BY FLOOR(model_to_map_fit/5.00)*5 # ORDER BY 1