rdkit-data 201603.5-2 / usr / share / RDKit / Contrib / pzc / p_con.py

# coding=utf-8
# Copyright (c) 2014 Merck KGaA
from __future__ import print_function
import os,re,gzip,json,requests,sys, optparse,csv
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import SDWriter
from rdkit.Chem import Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors
from scipy import interp
from scipy import stats
from sklearn import cross_validation
from sklearn.ensemble import RandomForestClassifier
from sklearn import metrics
from sklearn.cross_validation import train_test_split
from sklearn.metrics import roc_curve, auc
from sklearn.metrics import precision_score,recall_score
from sklearn import preprocessing
import cPickle
from pickle import Unpickler
import numpy as np
import math
from pylab import *
from sklearn.metrics import make_scorer


kappa_template = '''\
%(kind)s Kappa Coefficient
--------------------------------
Kappa %(kappa)6.4f
ASE %(std_kappa)6.4f
%(alpha_ci)s%% Lower Conf Limit %(kappa_low)6.4f
%(alpha_ci)s%% Upper Conf Limit %(kappa_upp)6.4f

Test of H0: %(kind)s Kappa = 0

ASE under H0 %(std_kappa0)6.4f
Z %(z_value)6.4f
One-sided Pr > Z %(pvalue_one_sided)6.4f
Two-sided Pr > |Z| %(pvalue_two_sided)6.4f
'''

'''
Weighted Kappa Coefficient
--------------------------------
Weighted Kappa 0.4701
ASE 0.1457
95% Lower Conf Limit 0.1845
95% Upper Conf Limit 0.7558

Test of H0: Weighted Kappa = 0

ASE under H0 0.1426
Z 3.2971
One-sided Pr > Z 0.0005
Two-sided Pr > |Z| 0.0010
'''


def int_ifclose(x, dec=1, width=4):
    '''helper function for creating result string for int or float

only dec=1 and width=4 is implemented

Parameters
----------
x : int or float
value to format
dec : 1
number of decimals to print if x is not an integer
width : 4
width of string

Returns
-------
xint : int or float
x is converted to int if it is within 1e-14 of an integer
x_string : str
x formatted as string, either '%4d' or '%4.1f'
'''
    xint = int(round(x))
    if np.max(np.abs(xint - x)) < 1e-14:
        return xint, '%4d' % xint
    else:
        return x, '%4.1f' % x


class KappaResults(dict):

    def __init__(self, **kwds):
        self.update(kwds)
        if not 'alpha' in self:
            self['alpha'] = 0.025
            self['alpha_ci'] = int_ifclose(100 - 0.025 * 200)[1]

        self['std_kappa'] = np.sqrt(self['var_kappa'])
        self['std_kappa0'] = np.sqrt(self['var_kappa0'])

        self['z_value'] = self['kappa'] / self['std_kappa0']

        self['pvalue_one_sided'] = stats.norm.sf(self['z_value'])
        self['pvalue_two_sided'] = self['pvalue_one_sided'] * 2

        delta = stats.norm.isf(self['alpha']) * self['std_kappa']
        self['kappa_low'] = self['kappa'] - delta
        self['kappa_upp'] = self['kappa'] + delta

    def __str__(self):
        return kappa_template % self


def cohens_kappa(table, weights=None, return_results=True, wt=None):
    '''Compute Cohen's kappa with variance and equal-zero test

Parameters
----------
table : array_like, 2-Dim
square array with results of two raters, one rater in rows, second
rater in columns
weights : array_like
The interpretation of weights depends on the wt argument.
If both are None, then the simple kappa is computed.
see wt for the case when wt is not None
If weights is two dimensional, then it is directly used as a weight
matrix. For computing the variance of kappa, the maximum of the
weights is assumed to be smaller or equal to one.
TODO: fix conflicting definitions in the 2-Dim case for
wt : None or string
If wt and weights are None, then the simple kappa is computed.
If wt is given, but weights is None, then the weights are set to
be [0, 1, 2, ..., k].
If weights is a one-dimensional array, then it is used to construct
the weight matrix given the following options.

wt in ['linear', 'ca' or None] : use linear weights, Cicchetti-Allison
actual weights are linear in the score "weights" difference
wt in ['quadratic', 'fc'] : use linear weights, Fleiss-Cohen
actual weights are squared in the score "weights" difference
wt = 'toeplitz' : weight matrix is constructed as a toeplitz matrix
from the one dimensional weights.

return_results : bool
If True (default), then an instance of KappaResults is returned.
If False, then only kappa is computed and returned.

Returns
-------
results or kappa
If return_results is True (default), then a results instance with all
statistics is returned
If return_results is False, then only kappa is calculated and returned.

Notes
-----
There are two conflicting definitions of the weight matrix, Wikipedia
versus SAS manual. However, the computation are invariant to rescaling
of the weights matrix, so there is no difference in the results.

Weights for 'linear' and 'quadratic' are interpreted as scores for the
categories, the weights in the computation are based on the pairwise
difference between the scores.
Weights for 'toeplitz' are a interpreted as weighted distance. The distance
only depends on how many levels apart two entries in the table are but
not on the levels themselves.

example:

weights = '0, 1, 2, 3' and wt is either linear or toeplitz means that the
weighting only depends on the simple distance of levels.

weights = '0, 0, 1, 1' and wt = 'linear' means that the first two levels
are zero distance apart and the same for the last two levels. This is
the sampe as forming two aggregated levels by merging the first two and
the last two levels, respectively.

weights = [0, 1, 2, 3] and wt = 'quadratic' is the same as squaring these
weights and using wt = 'toeplitz'.

References
----------
Wikipedia
SAS Manual

'''
    table = np.asarray(table, float) #avoid integer division
    agree = np.diag(table).sum()
    nobs = table.sum()
    probs = table / nobs
    freqs = probs #TODO: rename to use freqs instead of probs for observed
    probs_diag = np.diag(probs)
    freq_row = table.sum(1) / nobs
    freq_col = table.sum(0) / nobs
    prob_exp = freq_col * freq_row[:, None]
    assert np.allclose(prob_exp.sum(), 1)
    #print prob_exp.sum()
    agree_exp = np.diag(prob_exp).sum() #need for kappa_max
    if weights is None and wt is None:
        kind = 'Simple'
        kappa = (agree / nobs - agree_exp) / (1 - agree_exp)

        if return_results:
            #variance
            term_a = probs_diag * (1 - (freq_row + freq_col) * (1 - kappa))**2
            term_a = term_a.sum()
            term_b = probs * (freq_col[:, None] + freq_row)**2
            d_idx = np.arange(table.shape[0])
            term_b[d_idx, d_idx] = 0 #set diagonal to zero
            term_b = (1 - kappa)**2 * term_b.sum()
            term_c = (kappa - agree_exp * (1-kappa))**2
            var_kappa = (term_a + term_b - term_c) / (1 - agree_exp)**2 / nobs
            #term_c = freq_col * freq_row[:, None] * (freq_col + freq_row[:,None])
            term_c = freq_col * freq_row * (freq_col + freq_row)
            var_kappa0 = (agree_exp + agree_exp**2 - term_c.sum())
            var_kappa0 /= (1 - agree_exp)**2 * nobs

    else:
        if weights is None:
            weights = np.arange(table.shape[0])
        #weights follows the Wikipedia definition, not the SAS, which is 1 -
        kind = 'Weighted'
        weights = np.asarray(weights, float)
        if weights.ndim == 1:
            if wt in ['ca', 'linear', None]:
                weights = np.abs(weights[:, None] - weights) / \
                           (weights[-1] - weights[0])
            elif wt in ['fc', 'quadratic']:
                weights = (weights[:, None] - weights)**2 / \
                           (weights[-1] - weights[0])**2
            elif wt == 'toeplitz':
                #assume toeplitz structure
                from scipy.linalg import toeplitz
                #weights = toeplitz(np.arange(table.shape[0]))
                weights = toeplitz(weights)
            else:
                raise ValueError('wt option is not known')
        else:
            rows, cols = table.shape
            if (table.shape != weights.shape):
                raise ValueError('weights are not square')
        #this is formula from Wikipedia
        kappa = 1 - (weights * table).sum() / nobs / (weights * prob_exp).sum()
        #TODO: add var_kappa for weighted version
        if return_results:
            var_kappa = np.nan
            var_kappa0 = np.nan
            #switch to SAS manual weights, problem if user specifies weights
            #w is negative in some examples,
            #but weights is scale invariant in examples and rough check of source
            w = 1. - weights
            w_row = (freq_col * w).sum(1)
            w_col = (freq_row[:, None] * w).sum(0)
            agree_wexp = (w * freq_col * freq_row[:, None]).sum()
            term_a = freqs * (w - (w_col + w_row[:, None]) * (1 - kappa))**2
            fac = 1. / ((1 - agree_wexp)**2 * nobs)
            var_kappa = term_a.sum() - (kappa - agree_wexp * (1 - kappa))**2
            var_kappa *= fac

            freqse = freq_col * freq_row[:, None]
            var_kappa0 = (freqse * (w - (w_col + w_row[:, None]))**2).sum()
            var_kappa0 -= agree_wexp**2
            var_kappa0 *= fac

    kappa_max = (np.minimum(freq_row, freq_col).sum() - agree_exp) / \
                (1 - agree_exp)

    if return_results:
        res = KappaResults( kind=kind,
                    kappa=kappa,
                    kappa_max=kappa_max,
                    weights=weights,
                    var_kappa=var_kappa,
                    var_kappa0=var_kappa0
                    )
        return res
    else:
        return kappa

def to_table(data, bins=None):
    '''convert raw data with shape (subject, rater) to (rater1, rater2)

    brings data into correct format for cohens_kappa

    Parameters
    ----------
    data : array_like, 2-Dim
        data containing category assignment with subjects in rows and raters
        in columns.
    bins : None, int or tuple of array_like
        If None, then the data is converted to integer categories,
        0,1,2,...,n_cat-1. Because of the relabeling only category levels
        with non-zero counts are included.
        If this is an integer, then the category levels in the data are already
        assumed to be in integers, 0,1,2,...,n_cat-1. In this case, the
        returned array may contain columns with zero count, if no subject
        has been categorized with this level.
        If bins are a tuple of two array_like, then the bins are directly used
        by ``numpy.histogramdd``. This is useful if we want to merge categories.

    Returns
    -------
    arr : nd_array, (n_cat, n_cat)
        Contingency table that contains counts of category level with rater1
        in rows and rater2 in columns.

    Notes
    -----
    no NaN handling, delete rows with missing values

    This works also for more than two raters. In that case the dimension of
    the resulting contingency table is the same as the number of raters
    instead of 2-dimensional.

    '''

    data = np.asarray(data)
    n_rows, n_cols = data.shape
    if bins is None:
        #I could add int conversion (reverse_index) to np.unique
        cat_uni, cat_int = np.unique(data.ravel(), return_inverse=True)
        n_cat = len(cat_uni)
        data_ = cat_int.reshape(data.shape)
        bins_ = np.arange(n_cat+1) - 0.5
        #alternative implementation with double loop
        #tt = np.asarray([[(x == [i,j]).all(1).sum() for j in cat_uni]
        #                 for i in cat_uni] )
        #other altervative: unique rows and bincount
    elif np.isscalar(bins):
        bins_ = np.arange(bins+1) - 0.5
        data_ = data
    else:
        bins_ = bins
        data_ = data


    tt = np.histogramdd(data_, (bins_,)*n_cols)

    return tt[0], bins_



class p_con:
    """Class to create Models to classify Molecules active or inactive
    using threshold for value in training-data"""

    def __init__(self,acc_id=None,proxy={}):
        """Constructor to initialize Object, use proxy if neccessary"""
        self.request_data={"acc_id":acc_id,"proxy":proxy}
        self.acc_id = acc_id
        self.proxy = proxy
        self.model = []
        self.verbous = False


    def __str__(self):
        """String-Representation for Object"""
        self.request_data["cmpd_count"] = len(self.sd_entries)
        retString = ""
        for key in self.request_data.keys():
            retString += "%s: %s\n" % (key,self.request_data[key])
        return retString.rstrip()


    def step_0_get_chembl_data(self):
        """Download Compound-Data for self.acc_id, these are available in self.sd_entries afterwards"""
        def looks_like_number(x):
            """Check for proper Float-Value"""
            try:
                float(x)
                return True
            except ValueError:
                return False

        if self.acc_id.find("CHEMBL") == -1:
            self.target_data = requests.get("https://www.ebi.ac.uk/chemblws/targets/uniprot/{}.json".format(self.acc_id),proxies=self.proxy).json()

        else:
            self.target_data = {}
            self.target_data['target'] = {}
            self.target_data['target']['chemblId'] = self.acc_id
        
        self.chembl_id = self.target_data['target']['chemblId']
        self.request_data["chembl_id"] = self.target_data['target']['chemblId']
#        print self.target_data
        self.bioactivity_data = requests.get("https://www.ebi.ac.uk/chemblws/targets/{}/bioactivities.json".format(self.target_data['target']['chemblId']),proxies=self.proxy).json()

        ic50_skip=0
        ki_skip=0
        inhb_skip=0

        count=0
        non_homo=0
        self.dr={}
        i = 0
        x = len(self.bioactivity_data['bioactivities'] )

        for bioactivity in [record for record in self.bioactivity_data['bioactivities'] if looks_like_number(record['value']) ] :
         
            if i%100 == 0:
                sys.stdout.write('\r' + str(i) + '/' +str(x) + ' >          <\b\b\b\b\b\b\b\b\b\b\b')
            elif (i%100)%10==0:
                sys.stdout.write('|')
            sys.stdout.flush()
            i += 1
#            if i > 5000: break
            if bioactivity['organism'] != 'Homo sapiens':
                non_homo+=1
                continue
            if re.search('IC50', bioactivity['bioactivity_type']):
                if bioactivity['units'] != 'nM':
                    ic50_skip+=1
                    continue
            elif re.search('Ki', bioactivity['bioactivity_type']):
                ki_skip+=1
                continue
            elif re.search('Inhibition', bioactivity['bioactivity_type']):
                inhb_skip+=1
            else:
                continue

            self.cmpd_data =  requests.get("https://www.ebi.ac.uk/chemblws/compounds/{}.json".format(bioactivity['ingredient_cmpd_chemblid']),proxies=self.proxy).json()

            my_smiles = self.cmpd_data['compound']['smiles']
            bioactivity['Smiles']=my_smiles
            self.dr[count] = bioactivity
            count+=1

        SDtags = self.dr[0].keys()
        cpd_counter=0
        self.sd_entries = []
        for x in range(len(self.dr)):
            entry = self.dr[x]
            cpd = Chem.MolFromSmiles(str(entry['Smiles']))
            AllChem.Compute2DCoords(cpd)
            cpd.SetProp("_Name",str(cpd_counter))
            cpd_counter += 1
            for tag in SDtags: cpd.SetProp(str(tag),str(entry[tag]))
            self.sd_entries.append(cpd)
        return True


    def step_1_keeplargestfrag(self):
        """remove all smaller Fragments per compound, just keep the largest"""
        result=[]

        for cpd in self.sd_entries:
            fragments = Chem.GetMolFrags(cpd,asMols=True)
            list_cpds_fragsize = []
            for frag in fragments:
                list_cpds_fragsize.append(frag.GetNumAtoms())
            largest_frag_index = list_cpds_fragsize.index(max(list_cpds_fragsize))
            largest_frag = fragments[largest_frag_index]
            result.append(largest_frag)

        self.sd_entries = result
        return True


    def step_2_remove_dupl(self):
        """remove duplicates from self.sd_entries"""
        result = []
        all_struct_dict = {}
        for cpd in self.sd_entries:
            Chem.RemoveHs(cpd)
            cansmi = Chem.MolToSmiles(cpd,canonical=True)
            if not cansmi in all_struct_dict.keys():
                all_struct_dict[cansmi] = []
            all_struct_dict[cansmi].append(cpd)
        
        for entry in all_struct_dict.keys():
            if len(all_struct_dict[entry])==1:
                all_struct_dict[entry][0].SetProp('cansmirdkit',entry)
                result.append(all_struct_dict[entry][0])
        
        self.sd_entries=result
        return True


    def step_3_merge_IC50(self):
        """merge IC50 of duplicates into one compound using mean of all values if:
        min(IC50) => IC50_avg-3*IC50_stddev && max(IC50) <= IC50_avg+3*IC50_stddev && IC50_stddev <= IC50_avg"""
        np_old_settings = np.seterr(invalid='ignore') #dirty way to ignore warnings from np.std
        def get_mean_IC50(mol_list):
            IC50 = 0
            IC50_avg = 0
            for bla in mol_list:
                try:
                    IC50 +=  float(bla.GetProp("value"))
                except Exception:
                    print("no IC50 reported",bla.GetProp("_Name"))
            IC50_avg = IC50 / len(mol_list)
            return IC50_avg

        def get_stddev_IC50(mol_list):
            IC50_list = []
            for mol in mol_list:
                try:
                    IC50_list.append(round(float(mol.GetProp("value")),2))
                except Exception:
                    print("no IC50 reported",mol.GetProp("_Name"))
            IC50_stddev = np.std(IC50_list,ddof=1)
            return IC50_stddev,IC50_list

        result = []
        IC50_dict = {}
        for cpd in self.sd_entries:
            if not "cansmirdkit" in cpd.GetPropNames():
                Chem.RemoveHs(cpd)
                cansmi = Chem.MolToSmiles(cpd,canonical=True)
                cpd.SetProp('cansmirdkit',cansmi)
            cansmi = str(cpd.GetProp("cansmirdkit"))
            IC50_dict[cansmi]={}

        for cpd in self.sd_entries:
            cansmi = str(cpd.GetProp("cansmirdkit"))
            try:
                IC50_dict[cansmi].append(cpd)
            except Exception:
                IC50_dict[cansmi] = [cpd]
        for entry in IC50_dict:
            IC50_avg = str(get_mean_IC50(IC50_dict[entry]))
            IC50_stddev,IC50_list = get_stddev_IC50(IC50_dict[entry])
            IC50_dict[entry][0].SetProp("value_stddev",str(IC50_stddev))
            IC50_dict[entry][0].SetProp("value",IC50_avg)
            minimumvalue = float(IC50_avg)-3*float(IC50_stddev)
            maximumvalue = float(IC50_avg)+3*float(IC50_stddev)
            
            if round(IC50_stddev,1) == 0.0:
                result.append(IC50_dict[entry][0])
            elif IC50_stddev > float(IC50_avg):
                runawaylist = []
                for e in IC50_dict[entry]:
                    runawaylist.append(e.GetProp("_Name"))
                    print("stddev larger than mean", runawaylist, IC50_list, IC50_avg,IC50_stddev)
            elif np.min(IC50_list) < minimumvalue or np.max(IC50_list) > maximumvalue:
                pass
            else:
                result.append(IC50_dict[entry][0])
        
        self.sd_entries=result
        np.seterr(over=np_old_settings['over'],divide=np_old_settings['divide'],invalid=np_old_settings['invalid'],under=np_old_settings['under'])
        return True

        
    def step_4_set_TL(self,threshold,ic50_tag="value"):
        """set Property "TL"(TrafficLight) for each compound:
        if ic50_tag (default:"value") > threshold: TL = 0, else 1"""
        result = []
        i,j = 0,0
        for cpd in self.sd_entries:
            if float(cpd.GetProp(ic50_tag))> float(threshold):
                cpd.SetProp('TL','0')
                i += 1
            else:
                cpd.SetProp('TL','1')
                j += 1
            result.append(cpd)

        self.sd_entries = result
        if self.verbous: print("## act: %d, inact: %d" % (j,i))
        return True


    def step_5_remove_descriptors(self):
        """remove list of Properties from each compound (hardcoded)
        which would corrupt process of creating Prediction-Models"""
        sd_tags = ['activity__comment','alogp','assay__chemblid','assay__description','assay__type','bioactivity__type','activity_comment','assay_chemblid','assay_description','assay_type','bioactivity_type','cansmirdkit','ingredient__cmpd__chemblid','ingredient_cmpd_chemblid','knownDrug','medChemFriendly','molecularFormula','name__in__reference','name_in_reference','numRo5Violations','operator','organism','parent__cmpd__chemblid','parent_cmpd_chemblid','passesRuleOfThree','preferredCompoundName','reference','rotatableBonds','smiles','Smiles','stdInChiKey','synonyms','target__chemblid','target_chemblid','target__confidence','target__name','target_confidence','target_name','units','value_avg','value_stddev'] + ['value']
        result = []

        for mol in self.sd_entries:
            properties = mol.GetPropNames()
            for tag in properties:
                if tag in sd_tags: mol.ClearProp(tag)
            result.append(mol)

        self.sd_entries = result
        return True

    def step_6_calc_descriptors(self):
        """calculate descriptors for each compound, according to Descriptors._descList"""
        nms=[x[0] for x in Descriptors._descList]
        calc = MoleculeDescriptors.MolecularDescriptorCalculator(nms)
        for i in range(len(self.sd_entries)):
            descrs = calc.CalcDescriptors(self.sd_entries[i])
            for j in range(len(descrs)):
                self.sd_entries[i].SetProp(str(nms[j]),str(descrs[j]))
        return True

    
    def step_7_train_models(self):
        """train models according to trafficlight using sklearn.ensamble.RandomForestClassifier
        self.model contains up to 10 models afterwards, use save_model_info(type) to create csv or html
        containing data for each model"""
        title_line = ["#","accuracy","MCC","precision","recall","f1","auc","kappa","prevalence","bias","pickel-File"]
        self.csv_text= [title_line]

        TL_list = []
        property_list_list = []
        directory = os.getcwd().split("/")[-2:]
        dir_string  = ';'.join(directory)
        for cpd in self.sd_entries:
            property_list = []
            property_name_list = []
            prop_name = cpd.GetPropNames()
            for property in prop_name:
                if property not in ['TL','value']:
                    try:
                        f = float(cpd.GetProp(property))
                        if math.isnan(f) or math.isinf(f):
                            print("invalid: %s" % property)
                        
                    except ValueError:
                        print("valerror: %s" % property)
                        continue
                    property_list.append(f)
                    property_name_list.append(property)
                elif property == 'TL':
                    TL_list.append(int(cpd.GetProp(property)))
                else:
                    print(property)
                    pass
            property_list_list.append(property_list)
        dataDescrs_array = np.asarray(property_list_list)
        dataActs_array   = np.array(TL_list)

        for randomseedcounter in range(1,11):
                if self.verbous: 
                    print("################################")
                    print("try to calculate seed %d" % randomseedcounter)
                X_train,X_test,y_train,y_test = cross_validation.train_test_split(dataDescrs_array,dataActs_array,test_size=.4,random_state=randomseedcounter)
#            try:
                clf_RF     = RandomForestClassifier(n_estimators=100,random_state=randomseedcounter)
                clf_RF     = clf_RF.fit(X_train,y_train)

                cv_counter = 5

                scores = cross_validation.cross_val_score( clf_RF, X_test,y_test, cv=cv_counter,scoring='accuracy')

                accuracy_CV = round(scores.mean(),3)
                accuracy_std_CV = round(scores.std(),3)
   
                calcMCC = make_scorer(metrics.matthews_corrcoef,greater_is_better=True,needs_threshold=False)
                scores = cross_validation.cross_val_score( clf_RF, X_test,y_test, cv=cv_counter,scoring=calcMCC) 

                MCC_CV = round(scores.mean(),3)
                MCC_std_CV = round(scores.std(),3)

                scores = cross_validation.cross_val_score( clf_RF, X_test,y_test, cv=cv_counter,scoring='f1')
                scores_rounded = [round(x,3) for x in scores]
                f1_CV = round(scores.mean(),3)
                f1_std_CV = round(scores.std(),3)

                scores = cross_validation.cross_val_score( clf_RF, X_test,y_test, cv=cv_counter,scoring='precision')
                scores_rounded = [round(x,3) for x in scores]
                precision_CV = round(scores.mean(),3)
                precision_std_CV = round(scores.std(),3)

                scores = cross_validation.cross_val_score( clf_RF, X_test,y_test, cv=cv_counter,scoring='recall')
                scores_rounded = [round(x,3) for x in scores]
                recall_CV = round(scores.mean(),3)
                recall_std_CV = round(scores.std(),3)

                scores = cross_validation.cross_val_score( clf_RF, X_test,y_test, cv=cv_counter,scoring='roc_auc')
                scores_rounded = [round(x,3) for x in scores]
                auc_CV = round(scores.mean(),3)
                auc_std_CV = round(scores.std(),3)

                y_predict = clf_RF.predict(X_test)
                conf_matrix = metrics.confusion_matrix(y_test,y_predict)
#                coh_kappa = cohenskappa.cohens_kappa(conf_matrix)
                coh_kappa = cohens_kappa(conf_matrix)
                kappa = round(coh_kappa['kappa'],3)
                kappa_stdev = round(coh_kappa['std_kappa'],3)
            
                tp = conf_matrix[0][0]
                tn = conf_matrix[1][1]
                fp = conf_matrix[1][0]
                fn = conf_matrix[0][1]
                n = tn+fp
                p = tp+fn
                kappa_prevalence = round(float(abs(tp-tn))/float(n),3)
                kappa_bias = round(float(abs(fp-fn))/float(n),3)

                if self.verbous:
                    print("test:")
                    print("\tpos\tneg")
                    print("true\t%d\t%d" % (tp,tn))
                    print("false\t%d\t%d" % (fp,fn))
                    print(conf_matrix)
                    print("\ntrain:")
                    y_predict2 = clf_RF.predict(X_train)
                    conf_matrix2 = metrics.confusion_matrix(y_train,y_predict2)
                    tp2 = conf_matrix2[0][0]
                    tn2 = conf_matrix2[1][1]
                    fp2 = conf_matrix2[1][0]
                    fn2 = conf_matrix2[0][1]
                    print("\tpos\tneg")
                    print("true\t%d\t%d" % (tp2,tn2))
                    print("false\t%d\t%d" % (fp2,fn2))
                    print(conf_matrix2)                    

                result_string_cut = [randomseedcounter,
                                     str(accuracy_CV)+"_"+str(accuracy_std_CV),
                                     str(MCC_CV)+"_"+str(MCC_std_CV),
                                     str(precision_CV)+"_"+str(precision_std_CV),
                                     str(recall_CV)+"_"+str(recall_std_CV),
                                     str(f1_CV)+"_"+str(f1_std_CV),
                                     str(auc_CV)+"_"+str(auc_std_CV),
                                     str(kappa)+"_"+str(kappa_stdev),
                                     kappa_prevalence,kappa_bias,"model_file.pkl"]


                self.model.append(clf_RF)
                self.csv_text.append(result_string_cut)

#            except Exception as e:
#                print "got %d models" % len(self.model)
#                print e
#                sys.exit(-1)
#                break
        return True if len(self.model)>0 else False

    def save_model_info(self,outfile,mode="html"):
        """create html- or csv-File for models according to mode (default: "html")"""
        if mode=="csv":
            if not outfile.endswith(".csv"): outfile += ".csv"
            csv_file = open(outfile,"wb")
            csv_file_writer = csv.writer(csv_file,delimiter=";",quotechar=' ')
            for line in self.csv_text: csv_file_writer.writerow(line)
            csv_file.flush()
            csv_file.close()
        elif mode=="html":
            if not outfile.endswith(".html"): outfile += ".html"
            def lines2list(lines):
                return lines

            def list2html(data,act,inact):
                html_head = """<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<title></title>
<style type="text/css">
table {
  max-width: 100%;
  background-color: transparent;
}

th {
  text-align: left;
}

.table {
  width: 100%;
  margin-bottom: 20px;
}

.table > thead > tr > th,
.table > tbody > tr > th,
.table > tfoot > tr > th,
.table > thead > tr > td,
.table > tbody > tr > td,
.table > tfoot > tr > td {
  padding: 8px;
  line-height: 1.428571429;
  vertical-align: top;
  border-top: 1px solid #dddddd;
}

.table > thead > tr > th {MSC1013123
  vertical-align: bottom;
  border-bottom: 2px solid #dddddd;
}

.table > caption + thead > tr:first-child > th,
.table > colgroup + thead > tr:first-child > th,
.table > thead:first-child > tr:first-child > th,
.table > caption + thead > tr:first-child > td,
.table > colgroup + thead > tr:first-child > td,
.table > thead:first-child > tr:first-child > td {
  border-top: 0;
}

.table > tbody + tbody {
  border-top: 2px solid #dddddd;
}

.table .table {
  background-color: #ffffff;
}

.table-condensed > thead > tr > th,
.table-condensed > tbody > tr > th,
.table-condensed > tfoot > tr > th,
.table-condensed > thead > tr > td,
.table-condensed > tbody > tr > td,
.table-condensed > tfoot > tr > td {
  padding: 5px;
}

.table-bordered {
  border: 1px solid #dddddd;
}

.table-bordered > thead > tr > th,
.table-bordered > tbody > tr > th,
.table-bordered > tfoot > tr > th,
.table-bordered > thead > tr > td,
.table-bordered > tbody > tr > td,
.table-bordered > tfoot > tr > td {
  border: 1px solid #dddddd;
}

.table-bordered > thead > tr > th,
.table-bordered > thead > tr > td {
  border-bottom-width: 2px;
}

.table-striped > tbody > tr:nth-child(odd) > td,
.table-striped > tbody > tr:nth-child(odd) > th {
  background-color: #f9f9f9;
}

.table-hover > tbody > tr:hover > td,
.table-hover > tbody > tr:hover > th {
  background-color: #f5f5f5;
}

table col[class*="col-"] {
  position: static;
  display: table-column;
  float: none;
}

table td[class*="col-"],
table th[class*="col-"] {
  display: table-cell;
  float: none;
}

.table > thead > tr > .active,
.table > tbody > tr > .active,
.table > tfoot > tr > .active,
.table > thead > .active > td,
.table > tbody > .active > td,
.table > tfoot > .active > td,
.table > thead > .active > th,
.table > tbody > .active > th,
.table > tfoot > .active > th {
  background-color: #f5f5f5;
}

.table-hover > tbody > tr > .active:hover,
.table-hover > tbody > .active:hover > td,
.table-hover > tbody > .active:hover > th {
  background-color: #e8e8e8;
}

.table > thead > tr > .success,
.table > tbody > tr > .success,
.table > tfoot > tr > .success,
.table > thead > .success > td,
.table > tbody > .success > td,
.table > tfoot > .success > td,
.table > thead > .success > th,
.table > tbody > .success > th,
.table > tfoot > .success > th {
  background-color: #dff0d8;
}

.table-hover > tbody > tr > .success:hover,
.table-hover > tbody > .success:hover > td,
.table-hover > tbody > .success:hover > th {
  background-color: #d0e9c6;
}

.table > thead > tr > .danger,
.table > tbody > tr > .danger,
.table > tfoot > tr > .danger,
.table > thead > .danger > td,
.table > tbody > .danger > td,
.table > tfoot > .danger > td,
.table > thead > .danger > th,
.table > tbody > .danger > th,
.table > tfoot > .danger > th {
  background-color: #f2dede;
}

.table-hover > tbody > tr > .danger:hover,
.table-hover > tbody > .danger:hover > td,
.table-hover > tbody > .danger:hover > th {
  background-color: #ebcccc;
}

.table > thead > tr > .warning,
.table > tbody > tr > .warning,
.table > tfoot > tr > .warning,
.table > thead > .warning > td,
.table > tbody > .warning > td,
.table > tfoot > .warning > td,
.table > thead > .warning > th,
.table > tbody > .warning > th,
.table > tfoot > .warning > th {
  background-color: #fcf8e3;
}

.table-hover > tbody > tr > .warning:hover,
.table-hover > tbody > .warning:hover > td,
.table-hover > tbody > .warning:hover > th {
  background-color: #faf2cc;
}

@media (max-width: 767px) {
  .table-responsive {
    width: 100%;
    margin-bottom: 15px;
    overflow-x: scroll;
    overflow-y: hidden;
    border: 1px solid #dddddd;
    -ms-overflow-style: -ms-autohiding-scrollbar;
    -webkit-overflow-scrolling: touch;
  }
  .table-responsive > .table {
    margin-bottom: 0;
  }
  .table-responsive > .table > thead > tr > th,
  .table-responsive > .table > tbody > tr > th,
  .table-responsive > .table > tfoot > tr > th,
  .table-responsive > .table > thead > tr > td,
  .table-responsive > .table > tbody > tr > td,
  .table-responsive > .table > tfoot > tr > td {
    white-space: nowrap;
  }
  .table-responsive > .table-bordered {
    border: 0;
  }
  .table-responsive > .table-bordered > thead > tr > th:first-child,
  .table-responsive > .table-bordered > tbody > tr > th:first-child,
  .table-responsive > .table-bordered > tfoot > tr > th:first-child,
  .table-responsive > .table-bordered > thead > tr > td:first-child,
  .table-responsive > .table-bordered > tbody > tr > td:first-child,
  .table-responsive > .table-bordered > tfoot > tr > td:first-child {
    border-left: 0;
  }
  .table-responsive > .table-bordered > thead > tr > th:last-child,
  .table-responsive > .table-bordered > tbody > tr > th:last-child,
  .table-responsive > .table-bordered > tfoot > tr > th:last-child,
  .table-responsive > .table-bordered > thead > tr > td:last-child,
  .table-responsive > .table-bordered > tbody > tr > td:last-child,
  .table-responsive > .table-bordered > tfoot > tr > td:last-child {
    border-right: 0;
  }
  .table-responsive > .table-bordered > tbody > tr:last-child > th,
  .table-responsive > .table-bordered > tfoot > tr:last-child > th,
  .table-responsive > .table-bordered > tbody > tr:last-child > td,
  .table-responsive > .table-bordered > tfoot > tr:last-child > td {
    border-bottom: 0;
  }
}
</style>
</head>
<body>
<p style="padding-left:10px;padding-top:10px;font-size:200&#37;">Data for Models</p>
<p style="padding-left:10px;padding-right:10px;">"""

                html_topPlot_start = """<table style="vertical-align:top; background-color=#CCCCCC">
<tr align="left" valign="top"><td><img src="pieplot.png"></td><td><H3>Distribution</H3><font color="#00C000">active %d</font><br><font color="#FF0000">inactive %d</td><td>"""




                html_topPlot_bottom="""</td></tr></table>"""

                html_tableStart="""<table class="table table-bordered table-condensed">
<thead>
<tr>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
</tr>
</thead>
<tbody>"""

                html_tElements ="""
<tr bgcolor = "%s">
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td><a href="%s">model.pkl</a></td>
</tr>"""

                html_bottomPlot = """</tbody>
</table>
<img src="barplot.png"><br>"""


                html_foot ="""
</p>
</body>
</html>"""

                html_kappa_table_head="""<table class="table table-bordered table-condensed">
<thead>
<tr>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
<th>%s</th>
</tr>
</thead>
<tbody>"""

                html_kappa_table_element="""<tr bgcolor = "%s">
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td>%s</td>
<td><a href="%s">model.pkl</a></td>
</tr>"""

                html_kappa_table_bottom="""</tbody>
</table>
<img src="barplot.png"><br>"""


                best,worst = findBestWorst(data)
                html = []
                html.append(html_head)
                html.append(html_topPlot_start % (act,inact))
                html.append(html_topPlot_bottom)
                html.append(html_tableStart % tuple(data[0]))
                i = 0
                for l in data[1:len(data)]:
                    l_replaced = []
                    for elem in l:
                        elem_string = str(elem)
                        if elem_string.find("pkl")==-1: l_replaced.append(elem_string.replace("_","±"))
                        else: l_replaced.append(elem_string)

                    c = ""
                    if i == best: c = "#9CC089"
                    if i == worst: c = "#FF3333"

                    html.append(html_tElements % tuple([c] + l_replaced))
                    i += 1
                html.append(html_bottomPlot)
                html.append(html_foot)
                createBarPlot(data)
                return html


            def writeHtml(html,outf):
                outf_h = open(outf,'w')
                for block in html:
                    outf_h.write(block)
                outf_h.flush()
                outf_h.close()
                return


            def findBestWorst(data):
                auc = [float(x[6].split("_")[0]) for x in data[1:]]
                max_index,min_index = auc.index(max(auc)),auc.index(min(auc))
                return (max_index,min_index)


            def createPiePlot(cpds):
                def getActInact(cpds):
                    act,inact=0,0
                    for cpd in cpds:
                        if int(cpd.GetProp('TL'))==0: inact+=1
                        else: act+=1
                    return act,inact

                act_count,inact_count = getActInact(cpds)
                print("act/inact from TL's %d/%d" % (act_count,inact_count))
                fig = plt.figure(figsize=(2,2))
                pie = plt.pie([inact_count,act_count],colors=('r','g'))
                fig.savefig("pieplot.png",transparent=True)
                return act_count,inact_count


            def createBarPlot(data):

                def getLists(data,col):
                    accList = []
                    errList = []
                    for x in data[1:]:
                        if x[col].find("_")==-1: continue
                        if x[col].find(".pkl")!=-1:continue
                        spl = x[col].split("_")
                        accList.append(float(spl[0]))
                        errList.append(float(spl[1]))
                    return accList,errList

                def plotLists(cnt):
                    result=[]
                    clr = ['#DD1E2F','#EBB035','#06A2CB','#218559','#D0C6B1','#192823','#DDAACC']
#                    print ticks, list,errList,width
#                    print ticks
                    for i in range(1,cnt):
                        list,errList = getLists(data,i)
#                        print i,cnt,list,errList
                        result.append(ax.bar(ticks+width*i,list,width,color=clr[i-1],yerr=errList))
                    return result

                fig,ax = plt.subplots()
                fig.set_size_inches(15,6)
                ticks = np.arange(0.0,12.0,1.2)
                if len(self.model)==1: ticks = np.arange(0.0,1.0,1.5)
                width = 0.15
                plots = plotLists(8)
                ax.set_xticks(ticks+0.75)
                ax.set_xticklabels([str(x) for x in range(1,11,1)])
                ax.set_ylabel("Accuracy")
                ax.set_xlabel("# model")
                ax.set_xlim(-0.3,14)
                ax.set_ylim(-0.1,1.2)
                ax.legend(tuple(plots),[x for x in data[0][1:8]],'upper right')
                best,worst = findBestWorst(data)
                if len(self.model)>1:
                    ax.annotate("best",xy=(ticks[best],0.85),xytext=(ticks[best]+0.25,1.1),color="green")
                    ax.annotate("worst",xy=(ticks[worst],0.85),xytext=(ticks[worst]+0.25,1.10),color="red")
                fig.savefig("barplot.png",transparent=True)
                return

            act,inact = createPiePlot(self.sd_entries)
            lines = self.csv_text
            data = lines2list(lines)
            html = list2html(data,act,inact)
            writeHtml(html,outfile)
        return True

    def load_mols(self,sd_file):
        """load SD-File from .sdf, .sdf.gz or .sd.gz"""
        if sd_file.endswith(".sdf.gz") or sd_file.endswith(".sd.gz"):
            SDFile = Chem.ForwardSDMolSupplier(gzip.open(sd_file))
        else:
            SDFile = Chem.SDMolSupplier(sd_file)
        self.sd_entries = [mol for mol in SDFile]
        return True

    def save_mols(self,outfile,gzip=True):
        """create SD-File of current molecules in self.sd_entries"""
        sdw = Chem.SDWriter(outfile+".tmp")
        for mol in self.sd_entries: sdw.write(mol)
        sdw.flush()
        sdw.close()
        if not gzip:
            os.rename(outfile+".tmp",outfile)
            return
        f_in = open(outfile+".tmp", 'rb')
        f_out = gzip.open(outfile, 'wb')
        f_out.writelines(f_in)
        f_out.flush()
        f_out.close()
        f_in.close()
        os.remove(outfile+".tmp")
        return

    def save_model(self,outfile,model_number=0):
        """save Model to file using cPickle.dump"""
        cPickle.dump(self.model[model_number],file(outfile,"wb+"))
        return
        
    def load_models(self,model_files):
        """load model or list of models into self.model"""
        if type(model_files)==str: model_files = [model_files]
        i = 0
        for mod_file in model_files:
            model = open(mod_file,'r')
            unPickled = Unpickler(model)
            clf_RF = unPickled.load()
            self.model.append(clf_RF)
            model.close()
            i += 1
        return i

    def predict(self,model_number):
        """try to predict activity of compounds using giving model-Number"""
        if len(self.model)<=model_number:
            sys.stderr.write("\nModel-Number %d doesn't exist, there are just %d Models\n" % (model_number,len(self.model)))
            sys.exit(-1)
        descriptors = []
        active,inactive = 0,0

        for D in Descriptors._descList:
            descriptors.append(D[0])
        calculator = MoleculeDescriptors.MolecularDescriptorCalculator(descriptors)

        clf_RF = self.model[model_number]


        for sample in self.sd_entries:
            use = False
            try:
                pattern = calculator.CalcDescriptors(sample)
                use = True
            except e:
                sys.stderr.write("Error computing descriptors for %s, skip" % sample)
            
            if use:
                dataDescrs_array = np.asarray(pattern)
                y_predict  = int(clf_RF.predict(dataDescrs_array)[0])
                if y_predict==0: inactive += 1
                if y_predict==1: active += 1
                sample.SetProp("TL_prediction",str(y_predict))
        return (active,inactive)


if __name__ == "__main__":
    def step_error(step):
        sys.stderr.write("Error in Step: %s" % step)

    usage = "usage: python master.py [--accession=<Acc_ID>] [--sdf=<sdf-File>] --dupl/--uniq [--rof] [--combine=<file1>,<file2>] [--IC50=<IC50_tag>] [--cutoff=<value>] [--remove_descr=<txt_file>] [--proxy=<https://user:pass@proxy.de:portnumber] [--verbous] [--check_models=<model.pkl>]"
    parser = optparse.OptionParser(usage=usage)
    parser.add_option('--accession',action='store',type='string',dest='accession',help="Accession ID of Protein (hint: P43088 is Vitamin_D_Receptor with ~200 compounds)",default='')
    parser.add_option('--rof',action='store_true',dest='onefile',help='remove obsolete Files',default=False)
    parser.add_option('--dupl',action='store_true',dest='dupl',help='use only duplicates',default=False)
    parser.add_option('--uniq',action='store_true',dest='uniq',help='use only uniques',default=False)
    parser.add_option('--combine',action='store',type='string',dest='combine',help='Combine 2 SDF/SDF.GZ Files',default='')
    parser.add_option('--IC50',action='store',type='string',dest='SD_tag',help='name of IC50 field, default is \'value\'',default='value')
    parser.add_option('--cutoff',action='store',type='int',dest='cutoff',help='cutoff-value for hERG-trafficlight, default is \'5000\'',default=5000)
    parser.add_option('--remove_descr',action='store',type='string',dest='remove_descr',help='file with SDtags2remove, line-wise default:<internal list>',default='')
    parser.add_option('--proxy',action='store',type='string',dest='proxy',help='Use this Proxy',default='')
    parser.add_option('--sdf',action='store',type='string',dest='sdf',help='load this SDF-File',default='')
    parser.add_option('--verbous',action='store_true',dest='verbous',help='verbous',default=False)
    parser.add_option('--check_models',action='store',type='string',dest='modelfile',help='check compounds with this model',default='')

    (options,args) = parser.parse_args()
    combineItems = options.combine.split(',')

    if   len(combineItems) == 1 and len(combineItems[0])>0:
        print('need 2 files to combine')
        print(usage)
        sys.exit(-1)
    elif len(combineItems) == 2 and len(combineItems[0])>0 and len(combineItems[1])>0:
        cur_file = _04.combine(combineItems[0],combineItems[1])
        print("File: %s" % cur_file)
        sys.exit(0)

    code = options.accession.split(':')
    if len(code)==1:
        accession = code[0]
    else:
        accession = code[1]

    if options.accession == '' and options.sdf == '':
        print("please offer Accession-Number or SDF-File")
        print("-h for help")
        sys.exit(-1)
        
        
    if options.dupl==False and options.uniq==False:
        print("Please select uniq or dupl -h for help")
        print("-h for help")
        sys.exit(-1)


    pco = p_con(accession,proxy=options.proxy)
    pco.verbous = options.verbous


    if options.sdf != '':
        print("load sdf from File: %s" % options.sdf)
        result = pco.load_mols(options.sdf)
        if not result:
            step_error("load SDF-File")
            sys.exit(-1)
    else:
        print("gather Data for Accession-ID \'%s\'" % accession)
        result = pco.step_0_get_chembl_data()
        if not result:
            step_error("download ChEMBL-Data")
            sys.exit(-1)



    result = pco.step_1_keeplargestfrag()
    if not result:
        step_error("keep largest Fragment")
        sys.exit(-1)

    if options.uniq:
        result = pco.step_2_remove_dupl()
        if not result:
            step_error("remove duplicates")
            sys.exit(-1)

    result = pco.step_3_merge_IC50()
    if not result:
        step_error("merge IC50-Values for same Smiles")
        sys.exit(-1)

    if options.modelfile != '':
        result = pco.load_models(options.modelfile.split(","))
        if not result:
            step_error("Load Model-Files")
            sys.exit(-1)
        print("\n#Model\tActive\tInactive")
        for i in range(len(pco.model)):
            act,inact = pco.predict(i)
            print("%d\t%d\t%d" % (i,act,inact))
        sys.exit(0)

    result = pco.step_4_set_TL(options.cutoff)
    if not result:
        step_error("set Trafficlight for cutoff")
        sys.exit(-1)

    result = pco.step_5_remove_descriptors()
    if not result:
        step_error("remove descriptors")
        sys.exit(-1)

    result = pco.step_6_calc_descriptors()
    if not result:
        step_error("calculate Descriptors")
        sys.exit(-1)

    result = pco.step_7_train_models()
    if not result:
        step_error("Training of Models")
        sys.exit(-1)

    pco.save_model_info("model_info.csv",mode="csv")
    pco.save_model_info("model_info.html",mode="html")

    for i in range(len(pco.model)):
        filename = "%s_%dnm_model_%d.pkl" % (accession,options.cutoff,i)
        pco.save_model(filename,i)
        print("Model %d saved into File: %s" % (i,filename))
        

    for i in range(len(pco.model)):
        act,inact = pco.predict(i)
        print("Model %d active: %d\tinactive: %d" % (i,act,inact))