library(bio3d)
library(seqinr)
library(bios2mds)
library(protr)
#############################################################
#%% TASK
# use this to return df for AA index and mutation properties

source()

##############################################################
my_fasta_file = "~/git/Data/streptomycin/input/gid_complex.fasta"
my_mcsmf_snps = "~/git/Data/streptomycin/output/gid_mcsm_formatted_snps.csv"
###############################################################
#%% fasta as vector
gid_aa_seq_v= read.fasta(my_fasta_file
                    , seqtype = "AA"
                    , as.string = F)

gid_aa_v = as.character(gid_aa_seq_v[[1]]); gid_aa_v

#%% fasta as string
gid_aa_seq_s = read.fasta(my_fasta_file
                         , seqtype = "AA"
                         , as.string = T)

gid_aa_s = as.character(gid_aa_seq_s[[1]]); gid_aa_s
###############################################################
#===================
# AA indices
# https://www.genome.jp/aaindex/AAindex/list_of_indices
#===================
data(aa.index)

# default
aai_kd  = aa2index(gid_aa_v, index = "KYTJ820101") # Hydropathy, KD

aai_rv  = aa2index(gid_aa_v, index = "BIGC670101") # Residue volume, Bigelow, 1967
aai_rv2 = aa2index(gid_aa_v, index = "GOLD730102") # Residue volume (Goldsack-Chalifoux, 1973)
aai_b   = aa2index(gid_aa_v, index = "VENT840101") # Bitterness (Venanzi, 1984)

par(mfrow = c(1,1))
barplot(aai_kd)
barplot(aai_rv)
barplot(aai_rv2)
#barplot(aai_b, col = c("black", "yellow"))

##########################################################
#===================
# mutation matrices
#===================
data(sub.mat)
snps = read.csv(my_mcsmf_snps
                , header = 0)
snps
colnames(snps) <- "mutationinformation"

# run using all matrices
sub_mat_names = as.character(unlist(attributes(sub.mat)))
#sub_mat_names = "BLOSUM80"

for (j in sub_mat_names){
  print(j)
  snps[[j]] <- NA
for (i in 1:nrow(snps)) {
    curr_snp = snps$mutationinformation[i]
    m1 = str_match(curr_snp, "^([A-Z]{1})[0-9]*([A-Z]{1})")
    aa1 = m1[,2]
    aa2 = m1[,3]
    #snps$blosum_80[i]
    snps[[j]][i] = sub.mat[[j]][aa1,aa2]
  }

}
snps
##########################################################
gid_aac = extractAAC(gid_aa_s)
gid_dc = extractDC(gid_aa_s)
gid_tc = extractTC(gid_aa_s)

par(mfrow = c(1, 3))
barplot(gid_aac)
barplot(gid_dc)
barplot(gid_tc)
###########################################################