jquinone.py

import pandas as pd
import numpy as np
import re #regular expression
import itertools

import matplotlib.pyplot as plt
from mpltools import color

from collections import OrderedDict

from rdkit import Chem

#The following libraries are generated by Sungjin Kim
import jchem
import jutil

def gen_14BQ_OH():
    """
    return 1,4BQ species with OH functionals.
    """
    q_smiles_base = {}
    q_smiles_mid = {}
    q_smiles_base['1,4-BQ,2-OH'] = '[H]OC1=C([H])C(=O)C([H])=C([H])C1=O'
    q_smiles_base['1,4-BQ,Full-OH'] = 'OC1=C(O)C(=O)C(O)=C(O)C1=O'
    q_smiles_base['1,4-BQ'] = 'O=C1C=CC(=O)C=C1'

    q_smiles_mid['1,4-BQ'] = 'O=C1C=CC(=O)C=C1'
    q_smiles_mid['1,4-BQ,2-OH'] = 'OC1=CC(=O)C=CC1=O'
    q_smiles_mid['1,4-BQ,2,3-OH'] = 'OC1=C(O)C(=O)C=CC1=O'
    q_smiles_mid['1,4-BQ,2,3,5-OH'] = 'OC1=CC(=O)C(O)=C(O)C1=O'
    q_smiles_mid['1,4-BQ,Full-OH'] = 'OC1=C(O)C(=O)C(O)=C(O)C1=O'    

    return q_smiles_base, q_smiles_mid

def gen_910AQ_SO3H():
    """
    return 9,10AQ species with SO3H functionals.
    """
    q_smiles_base = {}
    q_smiles_mid = {}

    q_smiles_base['9,10AQ'] = 'O=C1C2C=CC=CC2C(=O)C2=C1C=CC=C2'
    q_smiles_base['9,10AQ,1-OH'] = 'OS(=O)(=O)C1=CC=CC2C1C(=O)C1=C(C=CC=C1)C2=O'
    q_smiles_base['9,10AQ,2-OH'] = 'OS(=O)(=O)C1=CC2C(C=C1)C(=O)C1=C(C=CC=C1)C2=O'
    q_smiles_base['9,10AQ,Full-OH'] = 'OS(=O)(=O)C1=C(C(=C(C2C1C(=O)C1=C(C2=O)C(=C(C(=C1S(O)(=O)=O)S(O)(=O)=O)S(O)(=O)=O)S(O)(=O)=O)S(O)(=O)=O)S(O)(=O)=O)S(O)(=O)=O'

    q_smiles_mid['9,10AQ'] = 'O=C1C2C=CC=CC2C(=O)C2=C1C=CC=C2'
    q_smiles_mid['9,10AQ,1-OH'] = 'OS(=O)(=O)C1=CC=CC2C1C(=O)C1=C(C=CC=C1)C2=O'
    q_smiles_mid['9,10AQ,2-OH'] = 'OS(=O)(=O)C1=CC2C(C=C1)C(=O)C1=C(C=CC=C1)C2=O'
    q_smiles_mid['9,10AQ,1,2-OH'] = 'OS(=O)(=O)C1=C(C2C(C=C1)C(=O)C1=C(C=CC=C1)C2=O)S(O)(=O)=O'
    q_smiles_mid['9,10AQ,Full-OH'] = 'OS(=O)(=O)C1=C(C(=C(C2C1C(=O)C1=C(C2=O)C(=C(C(=C1S(O)(=O)=O)S(O)(=O)=O)S(O)(=O)=O)S(O)(=O)=O)S(O)(=O)=O)S(O)(=O)=O)S(O)(=O)=O'

    return q_smiles_base, q_smiles_mid

def gen_smiles_quinone( quinone = '9,10AQ', r_group = 'SO3H'):
    if quinone == '1,4BQ' and r_group == 'OH':
        return gen_14BQ_OH()
    elif quinone == '9,10AQ' and r_group == 'SO3H':
        return gen_910AQ_SO3H()


class AQDS_OH():
    def __init__(self, fname = 'oh_subs_csv.csv'):
        self.pdr = pd.read_csv( fname)

    def get_Frag6_D( self, dfr = None):

        if dfr == None:
            dfr = self.pdr

        ri_vec = [1, 3, 4, 5, 6, 8]
        R = []
        HOH201 = {'H': 0, 'OH': 1}

        for ri in ri_vec:
            s = 'R{}'.format( ri)
            rv = dfr[s].tolist()
            rv_01 = [HOH201[x] for x in rv]
            R.append( rv_01)
        RM = np.mat( R).T

        Frag6_L = []   
        OHSYMB = {'H': '', 'OH': '(O)'}
        for ri in ri_vec:
            s = 'R{}'.format( ri)
            rv = dfr[s].tolist()
            fr_01 = [OHSYMB[x] for x in rv]
            Frag6_L.append( fr_01)
        #print Frag6_L

        Frag6_D = []
        for ii in range( len(Frag6_L[0])):
            Frag6_D.append({})
            
        for ii, frag in enumerate(Frag6_L):
            for ix, fr in enumerate(frag):
                dict_key = '{B%d}' % ii
                Frag6_D[ix][dict_key] = fr

        return Frag6_D

    def _gen_27aqds_with_oh_r1( self, Frag6_D, show = True):
        """
        2,7-AQDS with OH attachment are performed 
        using smiles interpolation
        """
        mol_smiles_list = []
        for ix, mol_symb in enumerate(Frag6_D):
            mol = bq14_oh2 = Chem.MolFromSmiles( 'C1(=O)c2c{B3}c{B4}c(S(=O)(=O)O)c{B5}c2C(=O)c2c{B0}c(S(=O)(=O)O)c{B1}c{B2}c21', replacements=mol_symb)
            mol_smiles = Chem.MolToSmiles( mol)
            mol_smiles_list.append( mol_smiles)

            if show:
                print(ix+1, mol_smiles)
                jchem.show_mol( mol_smiles)

        return mol_smiles_list

    def gen_27aqds_with_R( self, Frag6_D, r_gr, show = True):
        """
        2,7-AQDS with OH attachment are performed 
        using smiles interpolation
        """
        mol_smiles_list = []
        for ix, mol_symb in enumerate(Frag6_D):
            # r_gr = 'S(=O)(=O)O' #[N+]([O-])=O
            base_smiles = 'C1(=O)c2c{B3}c{B4}c(%s)c{B5}c2C(=O)c2c{B0}c(%s)c{B1}c{B2}c21' % (r_gr, r_gr)
            mol = bq14_oh2 = Chem.MolFromSmiles( base_smiles, replacements=mol_symb)
            mol_smiles = Chem.MolToSmiles( mol)
            mol_smiles_list.append( mol_smiles)

            if show:
                print(ix+1, mol_smiles)
                jchem.show_mol( mol_smiles)

        return mol_smiles_list

    def gen_27aqds_with_oh( self, Frag6_D, show = True):
        r_gr = 'S(=O)(=O)O' 
        return self.gen_27aqds_with_R( Frag6_D, r_gr, show = show)


    def gen_27aqds_with_no2( self, Frag6_D, show = True):
        r_gr = '[N+]([O-])=O'
        return self.gen_27aqds_with_R( Frag6_D, r_gr, show = show)

class HAQDS_OH():
    def __init__(self, fname = 'oh_subs_csv.csv'):
        self.pdr = pd.read_csv( fname)

    def get_Frag6_D( self, dfr = None):

        if dfr == None:
            dfr = self.pdr

        ri_vec = [1, 3, 4, 5, 6, 8]
        R = []
        HOH201 = {'H': 0, 'OH': 1}

        for ri in ri_vec:
            s = 'R{}'.format( ri)
            rv = dfr[s].tolist()
            rv_01 = [HOH201[x] for x in rv]
            R.append( rv_01)
        RM = np.mat( R).T

        Frag6_L = []   
        OHSYMB = {'H': '', 'OH': '(O)'}
        for ri in ri_vec:
            s = 'R{}'.format( ri)
            rv = dfr[s].tolist()
            fr_01 = [OHSYMB[x] for x in rv]
            Frag6_L.append( fr_01)
        #print Frag6_L

        Frag6_D = []
        for ii in range( len(Frag6_L[0])):
            Frag6_D.append({})
            
        for ii, frag in enumerate(Frag6_L):
            for ix, fr in enumerate(frag):
                dict_key = '{B%d}' % ii
                Frag6_D[ix][dict_key] = fr

        return Frag6_D

    def _gen_27aqds_with_oh_r1( self, Frag6_D, show = True):
        """
        2,7-AQDS with OH attachment are performed 
        using smiles interpolation
        """
        mol_smiles_list = []
        for ix, mol_symb in enumerate(Frag6_D):
            mol = bq14_oh2 = Chem.MolFromSmiles( 'C1(O)c2c{B3}c{B4}c(S(=O)(=O)O)c{B5}c2C(O)c2c{B0}c(S(=O)(=O)O)c{B1}c{B2}c21', replacements=mol_symb)
            mol_smiles = Chem.MolToSmiles( mol)
            mol_smiles_list.append( mol_smiles)

            if show:
                print(ix+1, mol_smiles)
                jchem.show_mol( mol_smiles)

        return mol_smiles_list

    def gen_27aqds_with_R( self, Frag6_D, r_gr, show = True):
        """
        2,7-AQDS with OH attachment are performed 
        using smiles interpolation
        """
        mol_smiles_list = []
        for ix, mol_symb in enumerate(Frag6_D):
            # r_gr = 'S(=O)(=O)O' #[N+]([O-])=O
            base_smiles = 'C1(O)c2c{B3}c{B4}c(%s)c{B5}c2C(O)c2c{B0}c(%s)c{B1}c{B2}c21' % (r_gr, r_gr)
            mol = Chem.MolFromSmiles( base_smiles, replacements=mol_symb)
            mol_smiles = Chem.MolToSmiles( mol)
            mol_smiles_list.append( mol_smiles)

            if show:
                print(ix+1, mol_smiles)
                jchem.show_mol( mol_smiles)

        return mol_smiles_list

    def gen_27aqds_with_oh( self, Frag6_D, show = True):
        r_gr = 'S(=O)(=O)O' 
        return self.gen_27aqds_with_R( Frag6_D, r_gr, show = show)


    def gen_27aqds_with_no2( self, Frag6_D, show = True):
        r_gr = '[N+]([O-])=O'
        return self.gen_27aqds_with_R( Frag6_D, r_gr, show = show)        

def get_r_list( N_Rgroup = 4, so3h = '(S(O)(=O)=O)', disp = False, pdForm = True):
    pdr_id, pdr_index, pdr_rgroups, pdr_no_r = [], [], [], []
    
    N_max_bin = '0b' + '1' * N_Rgroup
    for pos in range( int(N_max_bin, 2) + 1):
        pos_bin = bin( pos)[2:].rjust( N_Rgroup, '0')
        so_int_l = [int(x) for x in pos_bin]
        so_l = [so3h if x == 1 else '' for x in so_int_l ]
        no_r = sum( so_int_l)
        
        pdr_id.append( pos + 1)
        pdr_no_r.append( no_r)
        pdr_index.append( so_int_l)
        pdr_rgroups.append( so_l)
        
        if disp: print(pos, no_r, so_int_l, '==>', so_l)
        
    if pdForm:
        pdr = pd.DataFrame()
        pdr['ID'] = pdr_id
        pdr['Rgroup'] = [so3h] * len( pdr_id)   
        pdr['NoOfR'] = pdr_no_r
        pdr['Index'] = pdr_index
        pdr['Rgroups'] = pdr_rgroups
        return pdr
    else:
        return so_l

def get_multi_r_list( N_positions = 4, r_l = ['', '(S(O)(=O)=O)', '(O)'], disp = False, pdForm = True):
    """
    Multiple R-groups will be attached. 
    The '' attachment should be involved in the list of R-groups since 
    it is also one of the possible case.
    """
    pdr_id, pdr_index, pdr_rgroups, pdr_no_r = [], [], [], []
    
    # The number of the possible elements in product operation is length of R-groups    
    # N_positions reprensents the number of attachment positions. 
    Nr = len( r_l)
    so_int_l_all = itertools.product( list(range( Nr)), repeat = N_positions)
    for pos, so_int_l in enumerate(so_int_l_all):
        so_l = [ r_l[x] for x in so_int_l]
        no_r = jutil.count( so_int_l, 0, inverse = True)

        pdr_id.append( pos + 1)
        pdr_no_r.append( no_r)
        pdr_index.append( so_int_l)
        pdr_rgroups.append( so_l)

        if disp: print(pos, no_r, so_int_l, '==>', so_l)

    if pdForm:
        pdr = pd.DataFrame()
        pdr['ID'] = pdr_id
        if len( r_l) == 2 and '' in r_l:
            """
            If r_l consists of one R-group and blank, 
            R-group becomes the one R-group. 
            The empty position can be 0 or 1, which is support for 
            generalization although it usually located in 0. 
            """
            if r_l.index( '') == 0:
                pdr['Rgroup'] = [ r_l[1]] * len( pdr_id)   
            else:
                pdr['Rgroup'] = [ r_l[0]] * len( pdr_id)   
        else:
            pdr['Rgroup'] = ['Mix'] * len( pdr_id)   
        pdr['NoOfR'] = pdr_no_r
        pdr['Index'] = pdr_index
        pdr['Rgroups'] = pdr_rgroups
        return pdr
    else:
        return so_l     

def gen_r_attach( mol = 'Oc1nc(O)c2nc3c{0}c{1}c{2}c{3}c3nc2n1', so3h = '(S(O)(=O)=O)', disp = False, graph = False):
    """
    generate molecules with R group fragment
    """
    N_group = len( re.findall( '{[0-9]*}', mol)) # find number of R group positions

    pdr = get_r_list( N_group, so3h, disp = disp, pdForm = True)
    so_l = pdr['Rgroups'].tolist()

    aso_l = []
    for so in so_l:        
        aso = mol.format(*so)
        aso_l.append( aso)
        if disp: print(so, aso)
        if graph: jchem.show_mol( aso)

    pdr['SMILES'] = aso_l
    pdr['BaseMol'] = [aso_l[0]] * len( aso_l)
    pdr['BaseStr'] = [mol] * len( aso_l)

    return pdr

# [nH]1c2[nH]c3ccccc3[nH]c2c(=O)[nH]c1=O

def gen_rl_attach( mol = 'Oc1nc(O)c2nc3c{0}c{1}c{2}c{3}c3nc2n1', r_l = ['', '(S(O)(=O)=O)'], disp = False, graph = False):
    """
    generate molecules with R group fragment
    get_r_list becomes get_multi_r_list so as to generate molecules with multiple R-group attached.
    """
    N_group = len( re.findall( '{[0-9]*}', mol)) # find number of R group positions

    pdr = get_multi_r_list( N_group, r_l, disp = disp, pdForm = True)
    so_l = pdr['Rgroups'].tolist()

    aso_l = []
    for so in so_l:        
        aso = mol.format(*so)
        aso_l.append( aso)
        if disp: print(so, aso)
        if graph: jchem.show_mol( aso)

    pdr['SMILES'] = aso_l
    pdr['BaseMol'] = [aso_l[0]] * len( aso_l)
    pdr['BaseStr'] = [mol] * len( aso_l)

    return pdr  

def gen_rl_2attach( mol, mol_nH, r_l = ['', '(S(O)(=O)=O)'], disp = False, graph = False):
    """
    generate molecules with R group fragment
    get_r_list becomes get_multi_r_list so as to generate molecules with multiple R-group attached.
    Reduced (or hydrated) SMILES strings will be generated as well. 
    """
    N_group = len( re.findall( '{[0-9]*}', mol)) # find number of R group positions

    pdr = get_multi_r_list( N_group, r_l, disp = disp, pdForm = True)
    so_l = pdr['Rgroups'].tolist()

    aso_l = []
    aso_nH_l = []
    for so in so_l:        
        aso = mol.format(*so)
        aso_l.append( aso)

        aso_nH = mol_nH.format(*so)
        aso_nH_l.append( aso_nH)

        if disp: print(so, aso, aso_nH)
        if graph: 
            print("Oxidated molecule:")
            jchem.show_mol( aso)

            print("Hydrated molecule:")
            jchem.show_mol( aso_nH)

    # Storing canonical smiles strings
    pdr['SMILES'] = jchem.csmiles_l( aso_l)
    pdr['R-SMILES'] = jchem.csmiles_l( aso_nH_l)
    pdr['BaseMol'] = [aso_l[0]] * len( aso_l)
    pdr['BaseStr'] = [mol] * len( aso_l)

    return pdr  


def gen_r_attach_Alloxazine_R123457( so3h = '(S(O)(=O)=O)', disp = False, graph = False):
    """
    generate molecules with R group fragment
    """

    # n1{R5}c2nc3c{R1}c{R2}c{R3}c{R4}c3nc2c(=O)n{R7}c1=O
    # 
    N_group = 6 #R1234 5 7 -> 0123 4 5

    pdr = get_r_list( N_group, so3h, disp = disp, pdForm = True)
    so_l = pdr['Rgroups'].tolist()

    aso_l = []
    mol_l = []
    for so in so_l:        
        if so[4] != '' and so[5] != '':
            aso = 'n1{4}c2nc3c{0}c{1}c{2}c{3}c3nc2c(=O)n{5}c1=O'.format(*so)
            mol_l.append('n1{4}c2nc3c{0}c{1}c{2}c{3}c3nc2c(=O)n{5}c1=O')
        elif so[4] == '' and so[5] == '':
            aso = '[nH]1c2nc3c{0}c{1}c{2}c{3}c3nc2c(=O)[nH]c1=O'.format(*so[:4])
            mol_l.append('[nH]1c2nc3c{0}c{1}c{2}c{3}c3nc2c(=O)[nH]c1=O')
        elif so[4] == '':
            aso = '[nH]1c2nc3c{0}c{1}c{2}c{3}c3nc2c(=O)n{4}c1=O'.format(so[0],so[1],so[2],so[3], so[5])
            mol_l.append('[nH]1c2nc3c{0}c{1}c{2}c{3}c3nc2c(=O)n{4}c1=O')
        else: #so[5] == '':
            aso = 'n1{4}c2nc3c{0}c{1}c{2}c{3}c3nc2c(=O)[nH]c1=O'.format(*so[:5])
            mol_l.append('n1{4}c2nc3c{0}c{1}c{2}c{3}c3nc2c(=O)[nH]c1=O')

        aso_l.append( aso)
        if disp: print(so, aso)
        if graph: jchem.show_mol( aso)

    pdr['SMILES'] = aso_l
    pdr['BaseMol'] = [aso_l[0]] * len( aso_l)
    pdr['BaseStr'] = mol_l

    return pdr

def gen_r_attach_Flavins( mol = 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O', so3h = '(S(O)(=O)=O)', disp = False, graph = False):
    # jchem.show_mol( 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O'.format( so, so, so, so, so, so))
    return gen_r_attach( mol = mol, so3h = so3h, disp = disp, graph = graph)

def r_attach_Flavins( idx = [0,0,0,0,0,0], mol = 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O', Rg = '(S(O)(=O)=O)'):
    """
    It attachs a functional group in a specific position.
    """
    idx_R = [ Rg if x == 1 else '' for x in idx]
    return mol.format( *idx_R)

def r_2attach_Flavins( idx = [0,0,0,0,0,0], mol = 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O', mol_nH = '[nH]1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3[nH]c2c(=O)[nH]{0}c1=O', 
                        Rg = '(S(O)(=O)=O)'):
    """
    It attachs a functional group in a specific position.
    """
    idx_R = [ Rg if x == 1 else '' for x in idx]
    return mol.format( *idx_R), mol_nH.format( *idx_R)

# 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O'
def gen_rl_attach_Flavins( mol = 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O', 
                            r_l = ['', '(S(O)(=O)=O)', '(O)'], disp = False, graph = False):
    # jchem.show_mol( 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O'.format( so, so, so, so, so, so))
    return gen_rl_attach( mol = mol, r_l = r_l, disp = disp, graph = graph)

def gen_rl_2attach_Flavins( mol = 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O', mol_nH = '[nH]1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3[nH]c2c(=O)[nH]{0}c1=O',
                            r_l = ['', '(S(O)(=O)=O)', '(O)'], disp = False, graph = False):
    # jchem.show_mol( 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O'.format( so, so, so, so, so, so))
    return gen_rl_2attach( mol = mol, mol_nH = mol_nH, r_l = r_l, disp = disp, graph = graph)

def gen_rl_2attach_Alloxazine_R10( mol = 'O=C1N{0}C(=O)C2=NC3=C{1}C{2}=C{3}C{4}=C3N=C2N1{5}', mol_nH = 'O=C1N{5}C2NC3=C{4}C{3}=C{2}C{1}=C3NC2C(=O)N1{0}',
                            r_l = ['', '(S(O)(=O)=O)', '(O)'], disp = False, graph = False):
    """
    The Alloxazine is used as a frame molecules. 
    """
    return gen_rl_2attach( mol = mol, mol_nH = mol_nH, r_l = r_l, disp = disp, graph = graph)
    
def gen_rl_2attach_Alloxazine( mol = 'O=C1N{1}C(=O)C2=NC3=C{2}C{3}=C{4}C{5}=C3N=C2N1{0}', mol_nH = 'O=C1N{0}C2NC3=C{5}C{4}=C{3}C{2}=C3NC2C(=O)N1{1}',
                            r_l = ['', '(S(O)(=O)=O)', '(O)'], disp = False, graph = False):
    """
    The Alloxazine is used as a frame molecules. 
    R1 is NH fragment, which was r10 previously. 
    """
    return gen_rl_2attach( mol = mol, mol_nH = mol_nH, r_l = r_l, disp = disp, graph = graph)

def get_R_group_dict():
    R_group_dict = OrderedDict( [
            ('no group', '[H]'), 
            ('amine', 'N'),
            ('hydroxyl', 'O'), 
            ('methyl', 'C'),
            ('fluoro', 'F'),
            ('phosphonic acid', 'P(O)(O)(=O)'),
            ('sulfonic acid', 'S(=O)(=O)(O)'), 
            ('carboxylic acid', 'C(O)(=O)'),
            ('nitro', '[N+]([O-])(=O)')])

    return R_group_dict

def gen_rl_2attach_ThiophenoquinonesA( 
        mol_ox = "{0}C1=C({1})C2=C(S1)C(=O)C({3})=C({2})C2=O", 
        mol_rd = "OC1=C({3})C({2})=C(O)C2=C1SC({0})=C2{1}",
        r_l = None, 
        disp = False, graph = False):
    # jchem.show_mol( 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O'.format( so, so, so, so, so, so))

    if r_l is None:
        """
        If Rgroup list is not specified, a default list is used. 
        """
        r_l = list(get_R_group_dict().values())

    pdw = gen_rl_2attach( mol = mol_ox, mol_nH = mol_rd, r_l = r_l, disp = disp, graph = graph)

    # generate_sfp : position * no of r_group
    Index = pdw.Index.values

    N_positions = len( re.findall( '{[0-9]*}', mol_ox)) # find number of R group positions
    N_rgroups = len( r_l)
    
    sfp_list = list()
    for idx in Index:
        sfp = list()
        for np in idx:
            sfp_each = [0] * N_rgroups
            sfp_each[ np] = 1
            sfp.extend( sfp_each)
        sfp_list.append( sfp)   

    pdw['SimpleFingerprint'] = sfp_list

    return pdw


def gen_2attach_Alloxazine( fname = 'sheet/Alloxazine2_64.csv', rg_l = ['(O)', '(S(O)(=O)=O)', '(OC)', '(C)']):
    pdr_l = []
    for rg in rg_l:
        pdr = jquinone.gen_rl_2attach_Alloxazine( r_l = ['', rg])
        fname_each = fname[:-4] + '{}'.format( rg) + '.csv'
        print('The results are saved to', fname)
        pdr.to_csv( fname_each, index = False)
        pdr_l.append( pdr)

    pdr_all = pd.concat( pdr_l, ignore_index = True)
    pdr_all['ID'] = list(range(1, pdr_all.shape[0] + 1))        
    pdr_all.to_csv( fname, index = False)
    return pdr_all


def show_Alloxazine():
    print('Original Alloxazine')
    jchem.show_mol( 'O=C1NC(=O)C2=NC3=CC=CC=C3N=C2N1')
    jchem.show_mol( 'O=C1N{0}C(=O)C2=NC3=C{1}C{2}=C{3}C{4}=C3N=C2N1{5}'.format( '(O)','(O)','(O)','(O)','(O)','(O)'))

    print('Hydro Alloxazine')
    jchem.show_mol( 'O=C1NC2NC3=CC=CC=C3NC2C(=O)N1')
    jchem.show_mol( 'O=C1N{5}C2NC3=C{4}C{3}=C{2}C{1}=C3NC2C(=O)N1{0}'.format( '(O)','(O)','(O)','(O)','(O)','(O)'))


# [nH]1c2[nH]c3ccccc3[nH]c2c(=O)[nH]c1=O
def gen_r_attach_Flavins_nH( mol = '[nH]1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3[nH]c2c(=O)[nH]{0}c1=O', so3h = '(S(O)(=O)=O)', disp = False, graph = False):
    # jchem.show_mol( 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O'.format( so, so, so, so, so, so))
    return gen_r_attach( mol = mol, so3h = so3h, disp = disp, graph = graph)

def gen_r_attach_Alloxazine( mol = '[nH]1{5}c2nc3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O', so3h = '(S(O)(=O)=O)', disp = False, graph = False):
    # '[nH]1{5}c2nc3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O'
    return gen_r_attach( mol = mol, so3h = so3h, disp = disp, graph = graph)

def gen_r_attach_lowpot_Flavins( disp = False, graph = False):
    oh = '(O)'
    h = ''
    oc = '(OC)'

    rl = []
    rl.append(([h,oh, oh, oh, h, h], -0.47))
    rl.append(([oh, oh, h,h,h,h], -0.47))
    rl.append(([oh, oh, oh, oh, oh, h], -0.47))
    rl.append(([oh, oh, oh, oh, h, h], -0.51))
    rl.append(([h, oh, h, oh, h, h], -0.50))
    rl.append(([h, oh, h, h, h, h], -0.45))
    rl.append(([oh, oh, h, oh, oh, h], -0.50))
    rl.append(([h, oh, h, oh, oh, h], -0.46))
    rl.append(([oh, oh, h, oh, h, h], -0.53))

    rl.append(([h, oc, oc, oc, h, h], -0.48))
    rl.append(([oc, oc, oc, oc, h, h], -0.48))
    rl.append(([oc, oc, h, oc, h, h], -0.47))
    rl.append(([h, oc, h, oc, oc, h], -0.46))
    rl.append(([oc, oc, h, oc, oc, h], -0.50))

    BaseStr = 'n1c2[nH]{5}c3c{4}c{3}c{2}c{1}c3nc2c(=O)[nH]{0}c1=O'
    N_group = len( re.findall( '{[0-9]*}', BaseStr))
    emptyR = [''] * N_group
    BaseMol = BaseStr.format( *emptyR)

    smiles_l = [ BaseStr.format(*r[0]) for r in rl]

    pdr = pd.DataFrame()
    pdr['ID'] = list(range( 1, len( smiles_l) + 1))

    R_group_l = []
    Index_l = []
    NoOfR_l = []
    for r in rl:
        # Whether it is oh or oc family is determined
        r_oh_test = [ x == oh for x in r[0]]
        print(r[0], '-->', r_oh_test, '-->', any(r_oh_test))

        if any(r_oh_test): 
            r_type = oh
        else:
            r_type = oc
        R_group_l.append( r_type)

        r_groups = [ 0 if x == '' else 1 for x in r[0]]
        Index_l.append( r_groups)
        NoOfR_l.append( np.sum( r_groups))

    pdr['Rgroup'] = R_group_l # This is newly included.
    pdr['NoOfR'] = NoOfR_l 
    pdr['Index'] = Index_l
    pdr['Rgroups'] = [ r[0] for r in rl]
    pdr['SMILES'] = smiles_l
    pdr['BaseMol'] = [BaseMol] * len(rl)
    pdr['BaseStr'] = [BaseStr] * len(rl)
    pdr['RedoxPotential'] = [ r[1] for r in rl]

    for ix, s in enumerate( smiles_l):
        if disp: print(ix+1, s)
        if graph:
            jchem.show_mol( s)

    return pdr

def aq1x( x = '(S(=O)(=O)O)', disp = False):
    """
    It generates new quinone molecules with mono functional group attachment.
    For anthraquinone(AQ), only two attachment positions are unique and
    All the other position attachment generate a duplicated result.
    """

    MaxAttach = 2
    cs_l = []
    en_mol = [''] * MaxAttach
    for ix in range( MaxAttach):
        en_mol[ix] = x
        s = 'O=C1c2c{0}c{1}ccc2C(=O)c2ccccc21'.format( *en_mol)
        cs_l.append( jchem.csmiles( s))
        en_mol[ix] = ''

        if disp:
            print(ix)
            jchem.show_mol( s)


    return cs_l

def bq1x( x = '(S(=O)(=O)O)', disp = False):
    """
    It generates new quinone molecules with mono functional group attachment.
    For anthraquinone(AQ), only two attachment positions are unique and
    All the other position attachment generate a duplicated result.
    """

    MaxAttach = 1
    cs_l = []
    en_mol = [''] * MaxAttach
    for ix in range( MaxAttach):
        en_mol[ix] = x
        s = 'C1=CC(=O)C{0}=CC1=O'.format( *en_mol)
        cs_l.append( jchem.csmiles( s))
        en_mol[ix] = ''

        if disp:
            print(ix)
            jchem.show_mol( s)


    return cs_l 

def mol1x( mol = 'C1=CC(=O)C{0}=CC1=O', MaxAttach = 1, x = '(S(=O)(=O)O)', disp = False):
    """
    It generates new quinone molecules with mono functional group attachment.
    For anthraquinone(AQ), only two attachment positions are unique and
    All the other position attachment generate a duplicated result.
    """

    #MaxAttach = 1
    cs_l = []
    en_mol = [''] * MaxAttach
    for ix in range( MaxAttach):
        en_mol[ix] = x
        s = mol.format( *en_mol)
        cs_l.append( jchem.csmiles( s))
        en_mol[ix] = ''

        if disp:
            print(ix)
            jchem.show_mol( s)

    return cs_l 

def plot_coef( coef, intercept, shape = (9, 4) , loc = 2, R_g = None):
    print(coef[0].shape)
    W = np.reshape( coef[0], shape)
    print(W.shape)

    print("I. X-axis follows functional groups")
    plt.figure()
    for ii in range( shape[1]):
        plt.plot( W[:,ii], 'o', label = "R{}".format( ii+1))
    
    plt.plot( [0, shape[0]-1], [intercept, intercept])
    plt.xlabel( "Functional groups")
    plt.ylabel( "Property")
    plt.legend( loc = loc)
    plt.show()

    print(R_g)

    print("II. X-axis follows attached positions")
    plt.figure()
    if R_g:
        assert (len(R_g) == shape[0]), 'len(R_g) is not equal to coef.shape[0]'
        for ii in range( shape[0]):
            color = float(ii) / shape[0]
            plt.plot( list(range( shape[1])), W[ii, :], 'o', c = (1,color,1,1), ms = 10, label = R_g[ii])
            #plt.plot( x, y , '.', c = (1,0.1,1,1))
    else:
        for ii in range( shape[0]):
            plt.plot( W[ii, :], 'o', label = "F#{}".format( ii+1))
    
    plt.plot( [0, shape[1]-1], [intercept, intercept])
    plt.xlabel( "Attached positions")
    plt.ylabel( "Property")
    plt.legend( loc = loc)
    plt.show()


def anal_group_position_ridge( xM, yV, alpha = 0.5, shape = (9,4), loc = 1, R_g = None):
    """
    Show the relationship between groups & positions vs property

    """
    print("1. Testing regression with the given alpha using Ridge regression.")
    jutil.mlr_val_vseq_ridge_rand( xM, yV, alpha = alpha, rate = 5)

    print("2. Obtained regression coefficient with the same alpha.")
    coef, intercept = jutil.mlr3_coef_ridge( xM, yV, alpha = alpha)

    print("3. Showing the analyzed results.")
    plot_coef( coef, intercept, shape = shape, loc = loc, R_g = R_g)

    return coef, intercept

def pd_anal_ridge( pdr, y_id, alpha = 0.5, shape = (9,4), loc = 1, x_id = "SimpleFingerprint", R_g = None):
    """
    Analysis the property with respect to the functional groups and their positions. 

    """
    x = pdr[x_id].values
    xM = np.mat( list(map( eval, x)))

    yV = np.mat( pdr[y_id].values).T

    if R_g is None and shape[0] == 9:
        R_g = list(get_R_group_dict().keys())
        print(R_g)

    return anal_group_position_ridge( xM, yV, alpha = alpha, shape = shape, loc = loc, R_g = R_g)

def gen_pd_from_coef( coef0, Rgc):
    df = pd.DataFrame()
    r_l, p_l, c_l = [], [], []
    ci = 0
    for pi in range( 4):
        for ri in range( 9):
            r_l.append( Rgc[ri])
            p_l.append( pi + 1)
            c_l.append( coef0[ci])
            ci += 1   
    df['Position'] = p_l
    df['R-group'] = r_l    
    df['Coefficient'] = c_l
    return df