###################################################
### chunk number 1: 
###################################################
library(bgafun)
LDH <- read.alignment(file = system.file("sequences/LDH-MDH-PF00056.fasta", package = "bgafun"), format = "fasta")
class(LDH)


###################################################
### chunk number 2: 
###################################################
library(bgafun)
data(LDH)
LDH.amino=convert_aln_amino(LDH)
LDH.groups=rownames(LDH.amino)
LDH.groups[grep("LDH",LDH.groups)]="LDH"
LDH.groups[grep("MDH",LDH.groups)]="MDH"
LDH.groups=as.factor(LDH.groups)
LDH.groups


###################################################
### chunk number 3: 
###################################################
library(bgafun)
data(LDH)
data(LDH.groups)
LDH.amino=convert_aln_amino(LDH)
dim(LDH.amino)
LDH.amino.gapless=remove_gaps_groups(LDH.amino,LDH.groups)
dim(LDH.amino.gapless)


###################################################
### chunk number 4: 
###################################################
library(bgafun)
data(LDH.amino.gapless)
LDH.pseudo=LDH.amino.gapless+1
dim(LDH.pseudo)


###################################################
### chunk number 5: 
###################################################
library(bgafun)
data(LDH.amino.gapless)
LDH.pseudo=add_pseudo_counts(LDH.amino.gapless,LDH.groups)
dim(LDH.pseudo)


###################################################
### chunk number 6: 
###################################################
library(bgafun)
data(LDH)
data(LDH.groups)
LDH.aap=convert_aln_AAP(LDH)
dim(LDH.aap)
LDH.aap.ave=average_cols_aap(LDH.aap,LDH.groups)
dim(LDH.aap.ave)


###################################################
### chunk number 7: 
###################################################
library(bgafun)
data(LDH)
data(LDH.groups)
data(LDH.amino.gapless)
data(LDH.aap.ave)
LDH.aap.ave.bga=run_between_pca(LDH.amino.gapless,LDH.aap.ave,LDH.groups)
class(LDH.aap.ave.bga)



###################################################
### chunk number 8: 
###################################################
library(bgafun)
data(LDH.groups)
data(LDH.amino.gapless)
LDH.binary.bga=bga(t(LDH.amino.gapless+1),LDH.groups)


###################################################
### chunk number 9: PCAplot
###################################################
plot(LDH.aap.ave.bga)


###################################################
### chunk number 10: 
###################################################
top_res=top_residues_2_groups(LDH.binary.bga)
names(top_res)=sub("X","",names(top_res))


###################################################
### chunk number 11: 
###################################################
LDH.profiles=create_profile_strings(LDH.amino,LDH.groups)
LDH.profiles[, colnames(LDH.profiles) %in% names(top_res)]