gene = "katG"
drug = "isoniazid"

#==========
# LIGPLUS
#===========
# hem (1500)
aa_ligplus_inh = c(107, 108, 137, 229, 230)
#aa_ligplus_inh_hbond # none

aa_ligplus_hem = c(94, 276, 315, 274, 270, 381, 273, 104, 314, 275, 
                   100, 101, 321, 103, 269, 107, 266, 230, 380, 275, 314) 
                   
aa_ligplus_hem_hbond = c(94, 276, 315, 274, 270, 381)
aa_ligplus_hem_other = aa_ligplus_hem[!aa_ligplus_hem%in%aa_ligplus_hem_hbond]

c1 = length(aa_ligplus_hem_other) ==  length(aa_ligplus_hem) - length(aa_ligplus_hem_hbond)

#==========
# PLIP
#===========
aa_plip_inh = c(104, 229, 230)
aa_plip_inh_hbond = c(104, 229, 230)

aa_plip_hem = c(104, 107, 248, 252, 265, 275, 321, 412, 274, 276, 315)
aa_plip_hem_hbond = c(274, 276, 315)
#aa_plip_hem_sb = c(104, 276)
#aa_plip_hem_pi = c(107)
aa_plip_hem_other = aa_plip_hem[!aa_plip_hem%in%aa_plip_hem_hbond]

c2 = length(aa_plip_hem_other) ==  length(aa_plip_hem) - length(aa_plip_hem_hbond)

#==========
# Arpeggio
#===========
aa_arpeg_inh = c(104, 107, 108, 136, 137, 228, 229, 230, 232, 315) 
aa_arpeg_inh_hbond = c(104, 137)

aa_arpeg_hem = c(94, 100, 101, 103, 104, 107, 230, 231, 232, 248
                 , 252, 265, 266, 269, 270, 272, 273, 274, 275, 276, 314, 315
                 , 317, 321, 378, 380, 408, 412)

#from here

##############################################################
#===============
# Active site aa
#===============
active_aa_pos = sort(unique(c(aa_ligplus_inh
                              , aa_plip_inh
                              , aa_arpeg_inh
                              
                              , aa_ligplus_hem
                              , aa_plip_hem
                              , aa_arpeg_hem
                              )))
cat("\nNo. of active site residues for gene"
    , gene, ":"
    , length(active_aa_pos)
    , "\nThese are:\n"
    , active_aa_pos)

#=================
# Drug binding aa
#=================
aa_pos_inh = sort(unique(c(  aa_ligplus_inh
                           , aa_plip_inh
                           , aa_arpeg_inh)))
aa_pos_drug = aa_pos_inh


#===============
# Hbond aa
#===============
aa_pos_inh_hbond = sort(unique(c( aa_plip_inh_hbond
                           , aa_arpeg_inh_hbond)))

#=======================
# Other interactions aa
#=======================



#---------------------------------------------

aa_pos_hem = sort(unique(c(  aa_ligplus_hem
                           , aa_plip_hem
                           , aa_arpeg_hem)))

aa_pos_hem_hbond = sort(unique(c(  aa_ligplus_hem_hbond
                           , aa_plip_hem_hbond
                           #, aa_arpeg_hem_hbond
                           )))
                           

cat("\n==================================================="
    , "\nActive site residues for", gene, "comprise of..."
    , "\n==================================================="
    , "\nNo. of", drug, "binding residues:" , length(aa_pos_inh) , "\n"
    , aa_pos_inh
    , "\nNo. of 'HEM' binding residues:"    , length(aa_pos_hem) , "\n"
    , aa_pos_hem, "\n")

##############################################################
# var for position customisation for plots
aa_pos_lig1 = aa_pos_hem
aa_pos_lig2 = NULL
aa_pos_lig3 = NULL
tile_map=data.frame(tile=c("INH","HEME"),
                    tile_colour=c("green","darkslategrey"))



#toString(aa_pos_hem)
#toString(aa_pos_drug)
#toString(active_aa_pos)