gene = "alr"
drug = "cycloserine"

#==========
# LIGPLUS
#===========
aa_ligplus_dcs = c(66, 64, 70, 112, 196
                   , 236, 237, 252, 253
                   , 254, 255, 388)

aa_ligplus_dcs_hbond = c(255, 254, 237, 66, 196)
aa_ligplus_dcs_other = aa_ligplus_dcs[!aa_ligplus_dcs%in%aa_ligplus_dcs_hbond]

c1 = length(aa_ligplus_dcs_other) ==  length(aa_ligplus_dcs) - length(aa_ligplus_dcs_hbond)

#==========
# PLIP
#===========
aa_plip_dcs = c(66, 70, 112, 196, 237
                , 252, 254, 255, 295
                , 314, 343)
aa_plip_dcs_hbond = c(66, 70, 196, 237
                      , 252, 254, 255, 295
                      , 314, 343)

aa_plip_dcs_other = aa_plip_dcs[!aa_plip_dcs%in%aa_plip_dcs_hbond]

c2 = length(aa_plip_dcs_other) == length(aa_plip_dcs) - length(aa_plip_dcs_hbond)


#==========
# Arpeggio
#===========
aa_arpeg_dcs = c(64, 66, 70, 112, 157, 164
                 , 194, 196, 200, 236, 237, 252, 253
                 , 254, 255, 256, 295, 314, 342, 343
                 , 344, 386, 388)

aa_arpeg_dcs_other = aa_arpeg_dcs[!aa_arpeg_dcs%in%c(aa_ligplus_dcs_other
                                                     , aa_plip_dcs_other)]

c3 = length(aa_arpeg_dcs_other) == length(aa_arpeg_dcs) - ( length(aa_ligplus_dcs_other) + length(aa_plip_dcs_other) )

#######################################################################
#NEW AFTER ADDING PLP to structure! huh
# ADDED: 18 Aug 2022
# PLIP server for co factor PLP (CONFUSING!)
#and 2019 lit:lys42, M319, and Y364 : OFFSET is 24
#K42: K66, Y271:Y295, M319:M343, W89: W113, W203: W227, H209:H233, Q321:Q345
aa_pos_paper = sort(unique(c(66,70,112,113,164,196,227,233,237,252,254,255,295,342,343,344,345,388)))
plp_pos_paper = sort(unique(c(66, 70, 112, 196, 227, 237, 252, 254, 255, 388)))

#active_aa_pos = sort(unique(c(aa_pos_paper, active_aa_pos)))
aa_pos_plp = sort(unique(c(plp_pos_paper, 66, 70, 112, 237, 252, 254, 255, 196)))

#######################################################################

#===============
# Active site aa
#===============
active_aa_pos = sort(unique(c(aa_ligplus_dcs
                              , aa_plip_dcs
                              , aa_arpeg_dcs
                              , aa_pos_paper
                              , aa_pos_plp)))
#=================
# Drug binding aa
#=================
aa_pos_dcs = sort(unique(c(aa_ligplus_dcs
                           , aa_plip_dcs
                           , aa_arpeg_dcs)))
aa_pos_drug = aa_pos_dcs


#===============
# Co-factor: PLP aa
#===============
aa_pos_plp = aa_pos_plp

#===============
# Hbond aa
#===============
aa_pos_dcs_hbond = sort(unique(c(aa_ligplus_dcs_hbond
                                 , aa_plip_dcs_hbond)))

#=======================
# Other interactions aa
#=======================
aa_pos_dcs_other = active_aa_pos[!active_aa_pos%in%aa_pos_dcs_hbond]

c3 = length(aa_pos_dcs_other) == length(active_aa_pos) - length(aa_pos_dcs_hbond) 

#######################################################################
if ( all(c1, c2, c3) ) {
    
    cat("\nPASS:All active site residues and interctions checked and identified for"
        , "\ngene:", gene
        , "\ndrug:", drug
        , "\n==================================================="
        , "\nActive site residues for:", length(active_aa_pos) 
        , "\n==================================================="
        , "\n"
        , active_aa_pos
        
        , "\n=================================================="
        , "\nDrug binding residues:", length(aa_pos_drug)
        , "\n==================================================="
        , "\n"
        #, aa_pos_dcs
        , aa_pos_drug

        , "\n==================================================="
        , "\nHbond residues:", length(aa_pos_dcs_hbond)
        , "\n==================================================="
        , "\n"
        , aa_pos_dcs_hbond
        
        , "\n=================================================="
        , "\nOther interaction residues:",  length(aa_pos_dcs_other)
        , "\n==================================================="
        , "\n"
        , aa_pos_dcs_other
        , "\n\nNO other co-factors or ligands present\n")
        
}
######################################################################
#NEW
# PLIP server for co factor PLP (CONFUSING!)
#and 2019 lit:lys42, M319, and Y364 : OFFSET is 24
#K42: K66, Y271:Y295, M319:M343, W89: W113, W203: W227, H209:H233, Q321:Q345
aa_pos_paper = sort(unique(c(66,70,112,113,164,196,227,233,237,252,254,255,295,342,343,344,345,388)))
plp_pos_paper = sort(unique(c(66, 70, 112, 196, 227, 237, 252, 254, 255, 388)))

active_aa_pos = sort(unique(c(aa_pos_paper, active_aa_pos)))
aa_pos_plp = sort(unique(c(plp_pos_paper, 66, 70, 112, 237, 252, 254, 255, 196)))

# add two key residues
#aa_pos_drug = sort(unique(c(319, 364, aa_pos_drug)))
#active_aa_pos = sort(unique(c(319, 364, active_aa_pos, aa_pos_plp)))

# FIXME: these should be populated!
aa_pos_lig1 = aa_pos_plp
aa_pos_lig2 = NULL
aa_pos_lig3 = NULL

tile_map=data.frame(tile=c("ALR","PLP"),
                    tile_colour=c("green","navyblue")) #darkslategrey