source("~/git/LSHTM_analysis/scripts/Header_TT.R")
source("~/git/LSHTM_analysis/scripts/functions/ed_pfm_data.R")

# data msa: mut
my_data = read.csv("/home/tanu/git/Misc/practice_plots/pnca_msa_eg2.csv", header = F) #15 cols only
msaSeq_mut = my_data$V1
msa_seq = msaSeq_mut

# data msa: wt
gene = "pncA"
drug = "pyrazinamide"
indir = paste0("~/git/Data/", drug , "/input/")

in_filename_fasta = paste0(tolower(gene), "2_f2.fasta")
infile_fasta = paste0(indir, in_filename_fasta)
cat("\nInput fasta file for WT: ", infile_fasta, "\n")  

msa2 = read.csv(infile_fasta, header = F)
head(msa2)
cat("\nLength of WT fasta:", nrow(msa2))
wt_seq = msa2$V1
head(wt_seq)  
msaSeq_wt = msa2$V1
wt_seq = msaSeq_wt

################################
# DataED_PFM():
# script: ed_pfm_data.R
################################
data_ed = DataED_PFM(msa_seq, wt_seq)
names(data_ed)

#par(mfrow = c(2,1))
#logomaker(msa_seq, type = "EDLogo")
ggseqlogo(data_ed[['combED_mutM']]
          , method = "custom")