85 lines
No EOL
2.4 KiB
R
85 lines
No EOL
2.4 KiB
R
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)
|
|
########################################################### |