####Distancni matice
#degenerovane baze (ambiguous character/bases) predstavuji problem pri vypoctu dist.matic
#ve vypoctu je nekolik moznosti, jak s nimi nalozit:
#1) berou se jako uplne jiny znak (tzn.: Y neni ani C, ani T, ale samostatny znak). Toto dela Splitstree, ale mozna to je proto, ze nacte sekvence jako protein?!
#2) mohou byt jedna i druha baze (tzn.: Y muze byt jak C, tak T). Tento pristup snizuje diskriminaci a stahuje hybridy k rodicum (v podstate hybrid je jak rodic1, tak soucasne rodic2). Toto delaji v R napr.: dist.ml
#3) jsou dusledkem sekvencniho polymorfizmu=konsenzu mezi dvemi sekvencne odlisnymi alelami. Degenerovana baze znamena, ze nese prislusny podil informace (v Y je 1/2 C a 1/2 T). Hybridi jsou z poloviny rodic1 a z druhe poloviny rodic2. Toto by mela delat dist.p  
#4) degenerovany kod muze byt i ignorovan, ale pak vznkaji fatalni chyby.


#install.packages("ape")
#install.packages("phylotools")
#install.packages("phangorn")
#install.packages("maps")
#install.packages("phytools")
#install.packages("TreeTools")

library(ape)
library(phylotools)
library(phangorn)
library(maps)
library(phytools)
library(TreeTools)


###vypocet distancni matice
## phyDat ma defoltne nastavene, ze bere vsechny baze i degenerovany kod!

##nacteni dat
dat <- read.phyDat("C:/Users/kuba/Plocha/ITS-sekvence.txt", format="fasta")

#FASTA: dist.ml
#dist.ml nefunguje zcela spravne, kdyz se vyskytuji degenerovane baze (ambiguous character/bases)
#exclude= zda ma odstrani gapy, degenerovane baze apod.
DistMat1 <- dist.ml(dat, model = "JC69", exclude = "none")
write.nexus.dist(DistMat1, file="C:/Users/kuba/Plocha/OUTPUT1")

#FASTA: dist.p
DistMat2 <- dist.p(dat, cost = "polymorphism", ignore.indels = TRUE)
write.nexus.dist(DistMat2, file="C:/Users/kuba/Plocha/OUTPUT2")

DistMat1
DistMat2

#Heatmapa
##https://www.datanovia.com/en/lessons/heatmap-in-r-static-and-interactive-visualization/
##https://www.r-graph-gallery.com/215-the-heatmap-function.html
heatmap(x=as.matrix(DistMat1), Rowv=TRUE, Colv=TRUE, symm=TRUE)
heatmap(x=as.matrix(DistMat2), Rowv=TRUE, Colv=TRUE, symm=TRUE)
heatmap(x=as.matrix(DistMat2), Rowv=TRUE, Colv=NA, symm=TRUE)


###Jiný typ heatmapy
##Vytvoreni stromu z distancni matice
#NJ
myTreeNJ <- nj(DistMat2)
plotTree(myTreeNJ,ftype="i",fsize=0.6,lwd=1)

#preskladani NJ stromu
NewTree1 <- root(myTreeNJ, "Pal-3P")
plotTree(NewTree1,ftype="i",fsize=0.6,lwd=1)
NewTree2 <- SortTree(NewTree1, how = "TipLabels", order =c("Can-3A", "Can-2C", "X1", "X3", "X2","X5", "X6", "X4", "Pal-3P", "Pal-1Q"))
plotTree(NewTree2,ftype="i",fsize=0.6,lwd=1)


#UPGMA
myTreeUPGMA <- upgma(DistMat2)
plotTree(myTreeUPGMA,ftype="i",fsize=0.6,lwd=1)

#preskladani UPGMA stromu
NewTree3 <- root(myTreeUPGMA, "Pal-3P")
plotTree(NewTree3,ftype="i",fsize=0.6,lwd=1)
NewTree4 <- SortTree(myTreeUPGMA, how = "TipLabels", order =c("Can-3A", "Can-2C", "X1", "X3", "X2","X5", "X6", "X4", "Pal-3P", "Pal-1Q"))
plotTree(NewTree4,ftype="i",fsize=0.6,lwd=1)

##heatmapa a strom ML
#prevod distancni matice na matici
DM_ml <- as.matrix (DistMat1)

#heatmapa ML
phylo.heatmap(NewTree2, X=DM_ml, split=c(0.7,0.7),fsize=c(0.5,0.5,0.5), colors=NULL, standardize=FALSE, lwd=1, ftype="reg", pts=TRUE, labels=TRUE)

# uprava poradi vzorku dle stromu (vzorky jsou v matici dle abecedniho poradi - je potřeba je preskladat dle stromu)
DM_ml2 <- DM_ml[c(1,2,5,7,6,9,10,8,3,4), c(1,2,5,7,6,9,10,8,3,4)]

#heatmapa ML - oprava
phylo.heatmap(NewTree2, X=DM_ml2, split=c(0.7,0.7),fsize=c(0.5,0.5,0.5), colors=NULL, standardize=FALSE, lwd=1, ftype="reg", pts=TRUE, labels=TRUE)
phylo.heatmap(NewTree2, X=DM_ml2, split=c(0.7,0.7),fsize=c(0.5,0.5,0.5), colors=NULL, standardize=FALSE, lwd=1, ftype="reg", pts=TRUE, labels=FALSE)


##heatmapa a strom UPGMA
#prevod distancni matice na matici
DM_UPGMA <- as.matrix (DistMat2)
              
#heatmapa UPGMA
phylo.heatmap(NewTree4, X=DM_UPGMA, split=c(0.7,0.7),fsize=c(0.5,0.5,0.5), colors=NULL, standardize=FALSE, lwd=1, ftype="reg", pts=TRUE, labels=TRUE)
              
# uprava poradi vzorku dle stromu (vzorky jsou v matici dle abecedniho poradi - je potřeba je preskladat dle stromu)
DM_UPGMA2 <- DM_UPGMA[c("Can-3A", "Can-2C", "X1", "X3", "X2","X5", "X6", "X4", "Pal-3P", "Pal-1Q"), c("Can-3A", "Can-2C", "X1", "X3", "X2","X5", "X6", "X4", "Pal-3P", "Pal-1Q")]
              
#heatmapa UPGMA - oprava
phylo.heatmap(NewTree4, X=DM_UPGMA2, split=c(0.7,0.7),fsize=c(0.5,0.5,0.5), colors=NULL, standardize=FALSE, lwd=1, ftype="reg", pts=TRUE, labels=TRUE)
phylo.heatmap(NewTree4, X=DM_UPGMA2, split=c(0.7,0.7),fsize=c(0.5,0.5,0.5), colors=NULL, standardize=FALSE, lwd=1, ftype="reg", pts=TRUE, labels=FALSE)

