gene = "pncA"
drug = "pyrazinamide"

#===================================
#Iron centre --> purple
#Catalytic triad --> yellow 
#Substrate binding --> teal and blue
#H-bond --> green
#====================================
metal_aa_pos = c(49, 51, 57, 71) 
catalytic_aa_pos = c(8, 96, 138)
substrate_aa_pos = c(13, 68, 103, 137)
hbond_aa_pos = c(133, 134, 8, 138)

active_aa_pos = sort(unique(c(metal_aa_pos
                  , catalytic_aa_pos
                  , substrate_aa_pos
                  , hbond_aa_pos)))

cat("\nNo. of active site residues for gene"
    , gene, ":"
    , length(active_aa_pos)
    , "\nThese are:\n"
    , active_aa_pos)

cat("\n==================================================="
    , "\nActive site residues for", gene, "comprise of..."
    , "\n==================================================="
    , "\nMetal coordination centre residues:"
    , metal_aa_pos
    , "\nCatalytic triad residues:"
    , catalytic_aa_pos
    , "\nSubstrate binding residues:"
    , substrate_aa_pos
    , "\nH-bonding residues:"
    , hbond_aa_pos)