LSHTM_analysis/config/alr.R

176 lines
6.1 KiB
R

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)))
#######################################################################
# this is post inspection on chimera
#remove_pos = c(295, 314, 342, 343, 344)
remove_pos = c(0)
#select :295.A, 314.A, 342.A, 343.A, 344.A
#===============
# Active site aa
#===============
active_aa_pos = sort(unique(c(aa_ligplus_dcs
, aa_plip_dcs
, aa_arpeg_dcs
, aa_pos_plp)))
active_aa_pos = active_aa_pos[!active_aa_pos%in%remove_pos]
#=================
# Drug binding aa
#=================
aa_pos_dcs = sort(unique(c(aa_ligplus_dcs
, aa_plip_dcs
, aa_arpeg_dcs)))
aa_pos_dcs = aa_pos_dcs[!aa_pos_dcs%in%remove_pos]
aa_pos_drug = aa_pos_dcs
#===============
# Co-factor: PLP aa
#===============
aa_pos_plp = aa_pos_plp
#aa_pos_plp = aa_pos_plp[!aa_pos_plp%in%remove_pos]
#===============
# Hbond aa
#===============
aa_pos_dcs_hbond = sort(unique(c(aa_ligplus_dcs_hbond
, aa_plip_dcs_hbond)))
aa_pos_dcs_hbond = aa_pos_dcs_hbond[!aa_pos_dcs_hbond%in%remove_pos]
#=======================
# Other interactions aa
#=======================
aa_pos_dcs_other = active_aa_pos[!active_aa_pos%in%aa_pos_dcs_hbond]
aa_pos_dcs_other = aa_pos_dcs_other[!aa_pos_dcs_other%in%remove_pos]
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)))
#add_to_dcs = c(113, 227, 233, 345)
#add_to_plp = c(113, 227, 233, 345) # 227 not in plp and 227, 233 and 345 not with snp
#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_to_plp)))
aa_pos_plp = sort(unique(c(plp_pos_paper, 66, 70, 112, 237, 252, 254, 255, 196)))
#aa_pos_dcs = sort(unique(c(aa_pos_dcs, add_to_dcs)))
#aa_pos_drug = aa_pos_dcs
# 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("DCS","PLP"),
tile_colour=c("green","navyblue")) #darkslategrey
######
chain_suffix = ".A"
toString(paste0(aa_pos_drug, chain_suffix))
toString(paste0(aa_pos_plp, chain_suffix))
toString(paste0(active_aa_pos, chain_suffix))
common_pos = aa_pos_drug[aa_pos_drug%in%aa_pos_plp]
cat("\nCommon interacting partners:", length(common_pos))
common_pos
toString(paste0(common_pos, chain_suffix))