library(Biostrings)
library(readr)
library(DECIPHER)
library(seqinr)


load("image_applsci.Rdata")

nsm <- nucleotideSubstitutionMatrix(match = 5, mismatch = -4, baseOnly = FALSE, type = "DNA")

#amplicons936 <- sequences of amplicons generated with the 936 type primers
#ampliconsC2 <- sequences of amplicons generated with the c2 type primers
#ampliconsP335 <- sequences of amplicons generated with the P335 type primers



#BETWEEN SAMPLES PHAGE DIVERSITY ANALYSIS

# perform the multiple alignment with DECIPHER library
aligned_936 <- AlignSeqs(amplicons936,
                        guideTree = NULL,
                        iterations = 1,
                        refinements = 1,
                        gapOpening = -1,
                        gapExtension = -2,
                        useStructures = FALSE,
                        FUN = AdjustAlignment,
                        levels = c(0.9, 0.7, 0.7, 0.4, 10, 5, 5, 2),
                        processors = 1,
                        verbose = TRUE,
                        substitutionMatrix = nsm)

aligned_c2 <- AlignSeqs(ampliconsC2,
                         guideTree = NULL,
                         iterations = 1,
                         refinements = 1,
                         gapOpening = -1,
                         gapExtension = -2,
                         useStructures = FALSE,
                         FUN = AdjustAlignment,
                         levels = c(0.9, 0.7, 0.7, 0.4, 10, 5, 5, 2),
                         processors = 1,
                         verbose = TRUE,
                         substitutionMatrix = nsm)

aligned_P335 <- AlignSeqs(ampliconsP335,
                        guideTree = NULL,
                        iterations = 1,
                        refinements = 1,
                        gapOpening = -1,
                        gapExtension = -2,
                        useStructures = FALSE,
                        FUN = AdjustAlignment,
                        levels = c(0.9, 0.7, 0.7, 0.4, 10, 5, 5, 2),
                        processors = 1,
                        verbose = TRUE,
                        substitutionMatrix = nsm)

# calculate a maximum likelihood tree
dm_936 <- DistanceMatrix(aligned_936,
                    correction = "Jukes-Cantor")

dm_c2 <- DistanceMatrix(aligned_c2,
                         correction = "Jukes-Cantor")

dm_P335 <- DistanceMatrix(aligned_P335,
                        correction = "Jukes-Cantor")

dend_936 <- IdClusters(dm_936,
                   method = "ML",
                   type = "dendrogram",
                   myXStringSet = aligned_936)
dend_936 <- dendrapply(dend_936,
                   FUN = function(n) {
                     if(is.leaf(n)) 
                       attr(n, "label") <- 
                         as.expression(substitute(italic(leaf),
                                                  list(leaf = attr(n, "label"))))
                     n
                   })

dend_c2 <- IdClusters(dm_c2,
                       method = "ML",
                       type = "dendrogram",
                       myXStringSet = aligned_c2)
dend_c2 <- dendrapply(dend_c2,
                       FUN = function(n) {
                         if(is.leaf(n)) 
                           attr(n, "label") <- 
                             as.expression(substitute(italic(leaf),
                                                      list(leaf = attr(n, "label"))))
                         n
                       })

dend_P335 <- IdClusters(dm_P335,
                      method = "ML",
                      type = "dendrogram",
                      myXStringSet = aligned_P335)
dend_P335 <- dendrapply(dend_P335,
                      FUN = function(n) {
                        if(is.leaf(n)) 
                          attr(n, "label") <- 
                            as.expression(substitute(italic(leaf),
                                                     list(leaf = attr(n, "label"))))
                        n
                      })

# display the phylogenetic trees
plot(dend_936,
     yaxt = "n",
     horiz = TRUE)

plot(dend_c2,
     yaxt = "n",
     horiz = TRUE)

plot(dend_P335,
     yaxt = "n",
     horiz = TRUE)


#GLOBAL PHAGE DIVERSITY ANALYSIS

#preparation of local database  
#subsetBc2 <- list.files(pattern="*_seqdump.txt")#make a list of files with BLAST hit results
#sBLASTc2 <- readDNAStringSet(file = subsetBc2, format="fasta",
#                             nrec=-1L, skip=0L, seek.first.rec=FALSE,
#                             use.names=TRUE, with.qualities=FALSE) #
sBLASTc2 <- OrientNucleotides(sBLASTc2)#reorient all the sequences in the same direction
sBLASTc2 <- unique(sBLASTc2)#remove duplicates
total_c2 <- append(ampliconsC2, sBLASTc2, after=length(ampliconsC2))#combine amplicons seqs and BLAST hits
total_c2 <- OrientNucleotides(total_c2)

sBLAST936 <- OrientNucleotides(sBLAST936)#reorient all the sequences in the same direction
sBLAST936 <- unique(sBLAST936)#remove duplicates
total_936 <- append(amplicons936, sBLAST936, after=length(amplicons936))#combine amplicons seqs and BLAST hits
total_936 <- OrientNucleotides(total_936)

sBLASTP335 <- OrientNucleotides(sBLASTP335)#reorient all the sequences in the same direction
sBLASTP335 <- unique(sBLASTP335)#remove duplicates
total_P335 <- append(ampliconsP335, sBLASTP335, after=length(ampliconsP335))#combine amplicons seqs and BLAST hits
total_P335 <- OrientNucleotides(total_P335)

#align all the seqs
total_aligned_c2 <- AlignSeqs(total_c2,
                        guideTree = NULL,
                        iterations = 1,
                        refinements = 1,
                        gapOpening = -1,
                        gapExtension = -2,
                        useStructures = FALSE,
                        FUN = AdjustAlignment,
                        levels = c(0.9, 0.7, 0.7, 0.4, 10, 5, 5, 2),
                        processors = 1,
                        verbose = TRUE,
                        substitutionMatrix = nsm)

total_aligned_936 <- AlignSeqs(total_936,
                              guideTree = NULL,
                              iterations = 1,
                              refinements = 1,
                              gapOpening = -1,
                              gapExtension = -2,
                              useStructures = FALSE,
                              FUN = AdjustAlignment,
                              levels = c(0.9, 0.7, 0.7, 0.4, 10, 5, 5, 2),
                              processors = 1,
                              verbose = TRUE,
                              substitutionMatrix = nsm)

total_aligned_P335 <- AlignSeqs(total_P335,
                              guideTree = NULL,
                              iterations = 1,
                              refinements = 1,
                              gapOpening = -1,
                              gapExtension = -2,
                              useStructures = FALSE,
                              FUN = AdjustAlignment,
                              levels = c(0.9, 0.7, 0.7, 0.4, 10, 5, 5, 2),
                              processors = 1,
                              verbose = TRUE,
                              substitutionMatrix = nsm)

# calculate a maximum likelihood tree
globalDM_c2 <- DistanceMatrix(total_aligned_c2,
                      correction="Jukes-Cantor")

globalDM_936 <- DistanceMatrix(total_aligned_936,
                              correction="Jukes-Cantor")

globalDM_P335 <- DistanceMatrix(total_aligned_P335,
                              correction="Jukes-Cantor")

#PCA analysis for global diversity comparison
library(PCAtools)


all(is.finite(globalDM_c2[,c(1:length(ampliconsC2))]))#check for infinite values
PCA_c2 <- pca(globalDM_c2[,c(1:length(ampliconsC2))],
             metadata = NULL,
             center = TRUE,
             scale = TRUE,
             rank = NULL, 
             removeVar = NULL)

all(is.finite(globalDM_936[,c(1:length(amplicons936))]))#check for infinite values
PCA_936 <- pca(globalDM_936[,c(1:length(amplicons936))],
              metadata = NULL,
              center = TRUE,
              scale = TRUE,
              rank = NULL, 
              removeVar = NULL)

all(is.finite(globalDM_P335[,c(1:length(ampliconsP335))]))#check for infinite values
PCA_P335 <- pca(globalDM_P335[,c(1:length(ampliconsP335))],
              metadata = NULL,
              center = TRUE,
              scale = TRUE,
              rank = NULL, 
              removeVar = NULL)

options(ggrepel.max.overlaps = Inf)

screeplot(PCA_c2)
biplot(PCA_c2, lab = PCA_c2$yvars) #for figure
pairsplot(PCA_c2)
plotloadings(PCA_c2)

screeplot(PCA_936)
biplot(PCA_936, lab = PCA_936$yvars) #for figure
pairsplot(PCA_936)
plotloadings(PCA_936)

screeplot(PCA_P335)
biplot(PCA_P335, lab = PCA_P335$yvars) #for figure
pairsplot(PCA_P335)
plotloadings(PCA_P335)

#for total data on plot
PCA_c2 <- pca(globalDM_c2,
              metadata = NULL,
              center = TRUE,
              scale = TRUE,
              rank = NULL, 
              removeVar = NULL)

PCA_936 <- pca(globalDM_936,
               metadata = NULL,
               center = TRUE,
               scale = TRUE,
               rank = NULL, 
               removeVar = NULL)

PCA_P335 <- pca(globalDM_P335,
                metadata = NULL,
                center = TRUE,
                scale = TRUE,
                rank = NULL, 
                removeVar = NULL)

library(ggpubr)
ggarrange(biplot(PCA_c2, selectLab = PCA_c2$yvars[1:length(ampliconsC2)], lab = PCA_c2$yvars), 
          biplot(PCA_936, selectLab = PCA_936$yvars[1:length(amplicons936)], lab = PCA_936$yvars), 
          biplot(PCA_P335, selectLab = PCA_P335$yvars[1:length(ampliconsP335)], lab = PCA_P335$yvars), 
          labels = c("A", "B", "C"),
          ncol = 3, nrow = 1)
 
