gene = "rpoB"
drug = "rifampicin"

#==========
# LIGPLUS
#===========
# Error! No atom records found!

#==========
# PLIP
#===========
aa_plip_rfp = c(429, 432, 491, 487)
aa_plip_rfp_hbond = c(429, 432, 487)

# chainC: equivalent with offset (-6 from 5uhc) accounted
aa_plip_5uhc_rfp = c(430, 452, 483
                 , 491, 432, 433
                 , 448, 450, 459, 487)
aa_plip_5uhc_rfp_hbond = c(432, 433, 448, 450, 459, 487)

#==========
# Arpeggio
#===========
# rfp: 1894
aa_arpeg_rfp = c(170, 428, 429, 430, 431, 432
                 , 433, 435, 445, 448, 450, 452
                 , 453, 458, 483, 487, 491, 604
                 , 607, 674)

##############################################################
active_aa_pos = sort(unique(c(aa_plip_rfp
                              , aa_plip_5uhc_rfp
                              , aa_arpeg_rfp)))
##############################################################
cat("\nNo. of active site residues for gene"
    , gene, ":"
    , length(active_aa_pos)
    , "\nThese are:\n"
    , active_aa_pos)
##############################################################
aa_pos_rfp = sort(unique(c(aa_plip_rfp
                           , aa_plip_5uhc_rfp
                           , aa_arpeg_rfp)))

aa_pos_drug = aa_pos_rfp

aa_pos_rfp_hbond = sort(unique(c(aa_plip_rfp_hbond
                           , aa_plip_5uhc_rfp_hbond)))

cat("\n==================================================="
    , "\nActive site residues for", gene, "comprise of..."
    , "\n==================================================="
    , "\nNo. of", drug, "binding residues:" , length(aa_pos_rfp), "\n"
    , aa_pos_rfp
    , "\n\nNO other co-factors or ligands present\n")

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