xlogpdescriptor.java

化学图形处理软件

/*  $Revision: 9057 $ $Author: egonw $ $Date: 2007-10-14 20:35:55 +0200 (Sun, 14 Oct 2007) $
 *
 *  Copyright (C) 2005-2007  Christian Hoppe <[e-mail]>
 *
 *  Contact: [e-mail]
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public License
 *  as published by the Free Software Foundation; either version 2.1
 *  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 Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
 */
package org.openscience.cdk.qsar.descriptors.molecular;

import org._3pq.jgrapht.graph.SimpleGraph;
import org.openscience.cdk.CDKConstants;
import org.openscience.cdk.aromaticity.HueckelAromaticityDetector;
import org.openscience.cdk.exception.CDKException;
import org.openscience.cdk.graph.BFSShortestPath;
import org.openscience.cdk.graph.MoleculeGraphs;
import org.openscience.cdk.interfaces.*;
import org.openscience.cdk.isomorphism.UniversalIsomorphismTester;
import org.openscience.cdk.isomorphism.matchers.*;
import org.openscience.cdk.isomorphism.matchers.smarts.AnyOrderQueryBond;
import org.openscience.cdk.isomorphism.matchers.smarts.AromaticAtom;
import org.openscience.cdk.isomorphism.matchers.smarts.AromaticQueryBond;
import org.openscience.cdk.isomorphism.mcss.RMap;
import org.openscience.cdk.qsar.DescriptorSpecification;
import org.openscience.cdk.qsar.DescriptorValue;
import org.openscience.cdk.qsar.IMolecularDescriptor;
import org.openscience.cdk.qsar.result.DoubleResult;
import org.openscience.cdk.qsar.result.IDescriptorResult;
import org.openscience.cdk.ringsearch.AllRingsFinder;
import org.openscience.cdk.ringsearch.SSSRFinder;
import org.openscience.cdk.tools.manipulator.RingSetManipulator;

import java.util.ArrayList;
import java.util.Iterator;
import java.util.List;

/**
 * <p>Prediction of logP based on the atom-type method called XLogP. <b>Requires
 * all hydrogens to be explicit</b>.
 * <p>For description of the methodology see Ref. @cdk.cite{WANG97} and @cdk.cite{WANG00}
 * or <a href="http://www.chem.ac.ru/Chemistry/Soft/XLOGP.en.html">http://www.chem.ac.ru/Chemistry/Soft/XLOGP.en.html</a>.
 * Actually one molecular factor is missing (presence of para Hs donor pair).
 *
 * <p>This descriptor uses these parameters:
 * <table border="1">
 *   <tr>
 *     <td>Name</td>
 *     <td>Default</td>
 *     <td>Description</td>
 *   </tr>
 *   <tr>
 *     <td>checkAromaticity</td>
 *     <td>false</td>
 *     <td>True is the aromaticity has to be checked</td>
 *   </tr>
 *   <tr>
 *     <td>salicylFlag</td>
 *     <td>false</td>
 *     <td>True is to use the salicyl acid correction factor</td>
 *   </tr>
 * </table>
 *
 * <p>changed 2005-11-03 by chhoppe<br>
 *  -Internal hydrogen bonds are implemented<br>
 * CDK IDescriptor was validated against xlogp2.1<br>
 * As mentioned in the xlogP tutorial don't use charges, always draw bonds. To some extend we can support charges
 * but not in every case.
 * <p>CDK follows the program in following points (which is not documented in the paper):<br>
 * 	-Atomtyp 7 is -0.137<br>
 *  -Atomtype 81 is -0.447<br>
 *  -pi system does not consider P or S<br>
 *  -ring system >3<br>
 *  -aromatic ring systems >=6<br>
 *  -N atomtypes: (ring) is always (ring)c<br>
 *  -F 83 is not 0.375, the program uses 0.512 [2005-11-21]<br>
 *  -hydrophobic carbon is 1-3 relationship not 1-4 [2005-11-22]<br>
 *  -Atomtyp C 34/35/36 perception corrected [2005-11-22]; before Atomtyp perception ring perception is done -> slows run time<br>
 *
 *
 *
 *  <p>In question:<br>
 *  	-Correction factor for salicylic acid (in paper, but not used by the program)<br>
 *  	-Amid classification is not consequent (in 6 rings (R2)N-C(R)=0 is eg 46 and in !6 membered rings it is amid)<br>
 *  		-sometimes O=C(R)-N(R)-C(R)=O is an amid ... sometimes not<br>
 *		-Value for internal H bonds is in paper 0.429 but for no454 it is 0.643<br>
 *		-pi system defintion, the neighbourhood is unclear<br>
 *
 * <p>changed 2005-11-21 by chhoppe<br>
 * 	-added new parameter for the salicyl acid correction factor<br>
 *  -Corrected P and S perception for charges<br>
 *
 *
 *@author         mfe4, chhoppe
 *@cdk.created    2004-11-03
 *@cdk.module     qsar
 *@cdk.set        qsar-descriptors
 * @cdk.dictref qsar-descriptors:xlogP
 *
 * @cdk.keyword XLogP
 * @cdk.keyword descriptor
 */
public class XLogPDescriptor implements IMolecularDescriptor {

    private boolean checkAromaticity = false;
    private boolean salicylFlag=false;
    private SSSRFinder ssrf=null;
    /**
     *  Constructor for the XLogPDescriptor object.
     */
    public XLogPDescriptor() { }


    /**
     *  Gets the specification attribute of the XLogPDescriptor object.
     *
     *@return    The specification value
     */
    public DescriptorSpecification getSpecification() {
        return new DescriptorSpecification(
                "http://www.blueobelisk.org/ontologies/chemoinformatics-algorithms/#xlogP",
                this.getClass().getName(),
                "$Id: XLogPDescriptor.java 9057 2007-10-14 18:35:55Z egonw $",
                "The Chemistry Development Kit");
    }


    /**
     *  Sets the parameters attribute of the XLogPDescriptor object.
     *
     *@param  params            The new parameters value
     *@exception  CDKException  Description of the Exception
     *@see #getParameters
     */
    public void setParameters(Object[] params) throws CDKException {
        if (params.length != 2) {
            throw new CDKException("XLogPDescriptor expects two parameter");
        }
        if (!(params[0] instanceof Boolean)) {
            throw new CDKException("The first parameter must be of type Boolean");
        }else if(!(params[1] instanceof Boolean)) {
            throw new CDKException("The second parameter must be of type Boolean");
        }
        // ok, all should be fine
        checkAromaticity = ((Boolean) params[0]).booleanValue();
        salicylFlag=((Boolean) params[1]).booleanValue();
    }


    /**
     *Gets the parameters attribute of the XLogPDescriptor object.
     *
     *@return    The parameters value [boolean checkAromaticity, boolean salicylFlag]
     *@see #setParameters
     */
    public Object[] getParameters() {
        // return the parameters as used for the descriptor calculation
        Object[] params = new Object[2];
        params[0] = new Boolean(checkAromaticity);
        params[1] = new Boolean(salicylFlag);
        return params;
    }


    /**
     *  Calculates the xlogP for an atom container.
     *
     *  If checkAromaticity is true, the method check the aromaticity, if false, means that the aromaticity has
     *  already been checked. It is necessary to use before the call of this mehtod the
     *  addExplicitHydrogensToSatisfyValency method (HydrogenAdder classe).
     *
     *@param  atomContainer               AtomContainer
     *@return XLogP is a double
     *@exception CDKException  Possible Exceptions
     */

    public DescriptorValue calculate(IAtomContainer atomContainer) throws CDKException {
        IAtomContainer ac;
        try {
            ac = (IAtomContainer) atomContainer.clone();
        } catch (CloneNotSupportedException e) {
            throw new CDKException("Error during clone");
        }

        IRingSet rs = (IRingSet) (new AllRingsFinder()).findAllRings(ac);
        IRingSet atomRingSet=null;
        if (checkAromaticity) {
            HueckelAromaticityDetector.detectAromaticity(ac, rs, true);
        }
        double xlogP = 0;
//		SmilesParser sp = new SmilesParser(DefaultChemObjectBuilder.getInstance());
        String symbol = "";
        int bondCount = 0;
        int atomCount = ac.getAtomCount();
        int hsCount = 0;
        double xlogPOld=0;
        double maxBondOrder = 0;
        List hBondAcceptors=new ArrayList();
        List hBondDonors=new ArrayList();
        int checkAminoAcid=1;//if 0 no check, if >1 check
        IAtom atomi = null;
        for (int i = 0; i < atomCount; i++) {
            atomi = (IAtom)ac.getAtom(i);
            //			Problem fused ring systems
            atomRingSet= rs.getRings(atomi);
            atomi.setProperty("IS_IN_AROMATIC_RING", new Boolean(false));
            atomi.setProperty(CDKConstants.PART_OF_RING_OF_SIZE, new Integer(0));
            //logger.debug("atomRingSet.size "+atomRingSet.size());
            if (atomRingSet.getAtomContainerCount()>0){
                if (atomRingSet.getAtomContainerCount()>1){
                	Iterator containers = RingSetManipulator.getAllAtomContainers(atomRingSet).iterator();
                	atomRingSet = rs.getBuilder().newRingSet();
                	while (containers.hasNext()) {
                		ssrf = new SSSRFinder((IAtomContainer)containers.next());
                		atomRingSet.add(ssrf.findEssentialRings());
                	}
                    //logger.debug(" SSSRatomRingSet.size "+atomRingSet.size());
                }
                for (int j=0;j<atomRingSet.getAtomContainerCount();j++){
                    if (j==0){
                        atomi.setProperty(CDKConstants.PART_OF_RING_OF_SIZE, new Integer(((IRing)atomRingSet.getAtomContainer(j)).getRingSize()));
                    }

                    if (((IRing)atomRingSet.getAtomContainer(j)).contains(atomi)){
                        if (((IRing)atomRingSet.getAtomContainer(j)).getRingSize()>=6 && atomi.getFlag(CDKConstants.ISAROMATIC)){
                            atomi.setProperty("IS_IN_AROMATIC_RING", new Boolean(true));
                        }
                        if (((IRing)atomRingSet.getAtomContainer(j)).getRingSize()<((Integer)atomi.getProperty(CDKConstants.PART_OF_RING_OF_SIZE)).intValue()){
                            atomi.setProperty(CDKConstants.PART_OF_RING_OF_SIZE, new Integer(((IRing)atomRingSet.getAtomContainer(j)).getRingSize()));
                        }
                    }
                }
            }//else{
            //logger.debug();
            //}
        }


        for (int i = 0; i < atomCount; i++) {
            atomi = (IAtom)ac.getAtom(i);
            if (xlogPOld==xlogP & i>0 & !symbol.equals("H")){
                //logger.debug("\nXlogPAssignmentError: Could not assign atom number:"+(i-1));
            }

            xlogPOld=xlogP;
            symbol = atomi.getSymbol();
            bondCount = ac.getConnectedBondsCount(atomi);
            hsCount = getHydrogenCount(ac, atomi);
            maxBondOrder = ac.getMaximumBondOrder(atomi);
            if (!symbol.equals("H")){
                //logger.debug("i:"+i+" Symbol:"+symbol+" "+" bondC:"+bondCount+" Charge:"+atoms[i].getFormalCharge()+" hsC:"+hsCount+" maxBO:"+maxBondOrder+" Arom:"+atoms[i].getFlag(CDKConstants.ISAROMATIC)+" AtomTypeX:"+getAtomTypeXCount(ac, atoms[i])+" PiSys:"+getPiSystemsCount(ac, atoms[i])+" C=:"+getDoubleBondedCarbonsCount(ac, atoms[i])+" AromCc:"+getAromaticCarbonsCount(ac,atoms[i])+" RS:"+((Integer)atoms[i].getProperty(CDKConstants.PART_OF_RING_OF_SIZE)).intValue()+"\t");
            }
            if (symbol.equals("C")) {
                if (bondCount == 2) {
                    // C sp
                    if (hsCount >= 1) {
                        xlogP += 0.209;
                        //logger.debug("XLOGP: 38		 0.209");
                    } else {
                        if (maxBondOrder == 2.0) {
                            xlogP += 2.073;
                            //logger.debug("XLOGP: 40		 2.037");
                        } else if (maxBondOrder == 3.0) {
                            xlogP += 0.33;
                            //logger.debug("XLOGP: 39		 0.33");
                        }
                    }
                }
                if (bondCount == 3) {
                    // C sp2
                    if (((Boolean)atomi.getProperty("IS_IN_AROMATIC_RING")).booleanValue()) {
                        if (getAromaticCarbonsCount(ac, atomi) >= 2 && getAromaticNitrogensCount(ac,atomi)==0) {
                            if (hsCount == 0) {