### CV1 #########################################################

setwd("E:/Prednaska_cvika_statistika_2021")
dat<-read.table("CCv1Festuca.txt", header=T, sep = '\t')
names(dat)
str(dat)

#install.packages("openxlsx")
#library(openxlsx)
#dat <- read.xlsx("Cv1Festuca.xlsx", sheet = 1, range="A1:BC183")

aggregate(klasek.del.med ~ taxon, data = dat, min)
aggregate(cbind(klasek.del.med,klasek.del.max) ~ taxon, data = dat, min)


aggregate(. ~ taxon, data = dat[,-1], mean, na.rm=TRUE, na.action=NULL)
aggregate(. ~ taxon, data = dat[,-1], function(x) sd(x,na.rm=T)/mean(x,na.rm=T), na.action=NULL) #spocte CV

agg.quantile<-aggregate(. ~ taxon, data = dat[,-1], quantile, na.rm=TRUE, na.action=NULL) #spocita min, 25%, median, 75% a max
write.table(agg.quantile,"smazat.xls")


boxplot(klasek.del.med ~ taxon, data = dat)

par(mfrow=c(1,3))
boxplot(klasek.del.med ~ taxon, data = dat,main="Medián délky klásku (mm)", col=c("red","gold","dodgerblue","red","dodgerblue","red","gold"))
boxplot(osina.del.med ~ taxon, data = dat,main="Medián délky osin (mm)",col=c("red","gold","dodgerblue","red","dodgerblue","red","gold"))
boxplot(chlupy.list.delka.max ~ taxon, data = dat,main="Maximální délka chlupů na listech (mm)",col=c("red","gold","dodgerblue","red","dodgerblue","red","gold"))


hist(dat$steblo.del.med)
hist(dat$steblo.del.med,xlab="Medián délky stébla (cm)",ylab="Počet vzorků",main="",col="gold")


#Dva histogramy v jednom - popsano pekne zde: https://www.dataanalytics.org.uk/plot-two-overlapping-histograms-on-one-chart-in-r/
ax<-pretty(range(dat$osina.del.max,na.rm=TRUE),n=20)
hist.pa<-hist(dat[which(dat$taxon=="pa"),"osina.del.max"], breaks=ax,plot=FALSE)
hist.sc<-hist(dat[which(dat$taxon=="sc"),"osina.del.max"], breaks=ax,plot=FALSE)

plot(hist.pa, col="dodgerblue",xlab="Délka osiny (mm)",ylab="Počet vzorků",main="")
col2rgb("red")
plot(hist.sc, col=rgb(255,0,0, alpha=80, maxColorValue = 255),add=TRUE)

### CV2 ########################################################
setwd("G:/prednaskaMetody/Cvika/Cvika2021")
dat<-read.table("Cv1Festuca.txt", header=T, sep = '\t')


#_____ Normalita _____
dat$taxon<-as.factor(dat$taxon)
levels(dat$taxon)
dat$taxon2<-dat$taxon
levels(dat$taxon2)<-c("pal","psam","pal","pal","pal","pal","psam")
levels(dat$taxon2)

#vizualne
osina.pal<-dat[which(dat$taxon2=="pal"),"osina.del.med"]
plucha.pal<-dat[which(dat$taxon2=="pal"),"plucha.del.med"]
osina.psam<-dat[which(dat$taxon2=="psam"),"osina.del.med"]
plucha.psam<-dat[which(dat$taxon2=="psam"),"plucha.del.med"]

par(mfrow=c(2,2))
hist(osina.pal)
hist(osina.psam)
hist(plucha.pal)
hist(plucha.psam)

qqnorm(osina.pal)
qqline(osina.pal)
qqnorm(osina.psam)
qqline(osina.psam)
qqnorm(plucha.pal)
qqline(plucha.pal)
qqnorm(plucha.psam)
qqline(plucha.psam)

library(car) # v tomhle balicku to u qqplotu vykresli i konfidencni intervaly
qqPlot(osina.pal)
qqPlot(osina.psam)
qqPlot(plucha.pal)
qqPlot(plucha.psam)

shapiro.test(osina.pal) #prisny Shapiro-Wilk test
shapiro.test(osina.psam)
shapiro.test(plucha.pal)
shapiro.test(plucha.psam)

library(nortest)
lillie.test(osina.pal) #Lilliefors (Kolmogorov-Smirnov) test normality
lillie.test(osina.psam)
lillie.test(plucha.pal)
lillie.test(plucha.psam)


#library(nortest)
pearson.test(osina.pal) #Chi-square test (Pearsonuv) normality
pearson.test(osina.psam)
pearson.test(plucha.pal)
pearson.test(plucha.psam)


#_____ Testy rozdil delka osin _____

res<-vector(length=6) #na zacatek si udelam tabulku, kam si budu psat vysledky 
names(res)<-c("t-test","ANOVA","linear model","Mann-Whitney","Kruskal-Wallis","Chi-square")


tt<-t.test(osina.del.med~taxon2, data=dat) #t test
tt
res["t-test"]<-tt$p.value

an<-aov(osina.del.med~taxon2, data=dat) #ANOVA
anova(an)
res["ANOVA"]<-anova(an)[1,5]

li<-lm(osina.del.med~taxon2, data=dat) #linearni model
anova(li)
res["linear model"]<-anova(li)[1,5]

mw<-wilcox.test(osina.del.med~taxon2, data=dat) #Mann-Whitney neparametricky test
mw
res["Mann-Whitney"]<-mw$p.value

kw<-kruskal.test(osina.del.med~taxon2, data=dat) #Kruskal-Wallis neparametricka ANOVA
kw
res["Kruskal-Wallis"]<-kw$p.value

osina.mean<-mean(dat$osina.del.med,na.rm=TRUE);osina.mean
dat$osina.del.01<-ifelse(dat$osina.del.med>osina.mean,"long","short")
table (dat$osina.del.01,dat$taxon2) #tohle neni k testu nutny, jen to ukazuje, jakou tabulku pak test pocita
chi<-chisq.test(dat$osina.del.01,dat$taxon2)	#Chi-square test
chi
res["Chi-square"]<-chi$p.value

print(res)

#_____ ANOVA delka pluch _____
plch<-lm(plucha.del.med~taxon, data=dat) #radsi pouzivam lm nez aov - lm mi spocita odhady prumeru pro skupiny pomoci summary
anova(plch)
summary(plch)
plch1<-lm(plucha.del.med~taxon-1, data=dat) #to same, jen bez interceptu abych pak videl lip odhady pro taxony
summary(plch1)

boxplot(plucha.del.med~taxon, data=dat)
points(1:7, summary(plch1)$coefficients[,1],pch=15,col="red") #prida do boxplotu odhadovane prumery (boxplot tak kresli defaultne jen mediany)

#_____ post hoc + boxplot _____
plch2<-aov(plucha.del.med~taxon-1, data=dat) #tady je pro zmenu potreba pouzit aov, jinak se odmita pocitat post-hoc test
posthoc <- TukeyHSD(plch2, "taxon", conf.level=0.95,data=dat) # Tukey test pro rozdíly mezi skupinami
posthoc

library(multcompView
) #knihovnu je potreba nainstalovat
letters<- multcompLetters(posthoc$taxon[,4]) #přiřazení písmen jednotlivým vzorkům podle shody v Tukey testu 
letters$Letters<-letters$Letters[levels(dat$taxon)]

#nakreslení boxplotu
plot(plucha.del.med~taxon, data=dat,col=c("tomato2","gold","dodgerblue","tomato2","dodgerblue","tomato2","gold"),ylab="",xlab="",main="Medián délky pluch (mm)",ylim=c(3.3,6.5))
text(1:length(levels(dat$taxon)), 3.3, labels=letters$Letters, srt=90)
legend("topright", c("psamofilní diploidi", "pallens s.l. diploidi", "pallens s.l. tetraploidi"), col = c("gold","dodgerblue","tomato2"), pch = 15,pt.cex=1.8)

## To same jako predtim, je se serazenymi skupinami podle velikosti 
dato<-dat
dato$taxon<-with(dato, reorder(taxon, plucha.del.med, function(x) mean(x,na.rm=TRUE))) #tu funkci je potreba napsat sloziteji, aby fungovala i s NA hodnotama

plch3<-aov(plucha.del.med~taxon-1, data=dato)
posthoc <- TukeyHSD(plch3, "taxon", conf.level=0.95,data=dato)
letters<- multcompLetters(posthoc$taxon[,4])
letters$Letters<-letters$Letters[levels(dato$taxon)]

plot(plucha.del.med~taxon, data=dato,col=c("gold","gold","dodgerblue","tomato2","tomato2","tomato2","dodgerblue"),ylab="",xlab="",main="Medián délky pluch (mm)",ylim=c(3.3,6.5))
text(1:length(levels(dato$taxon)), 3.3, labels=letters$Letters, srt=90)
legend("topleft", c("psamofilní diploidi", "pallens s.l. diploidi", "pallens s.l. tetraploidi"), col = c("gold","dodgerblue","tomato2"), pch = 15,pt.cex=1.8)


#_____ Permutacni test _____
levels(dat$taxon)
dat$taxon2<-dat$taxon
levels(dat$taxon2)<-c("pal","psam","pal","pal","pal","pal","psam")
levels(dat$taxon2)

dat<- dat[complete.cases(dat$cev.svazky.max), ] #tohle vyhodi radky, kde je v delce klasku NA - lip se s tim pak bude pocitat
pal<-dat[which(dat$taxon2=="pal"),"cev.svazky.max"]
psam<-dat[which(dat$taxon2=="psam"),"cev.svazky.max"]
mean.pal=mean(pal, na.rm=TRUE)
mean.psam=mean(psam, na.rm=TRUE)
mean.diff=abs(mean.pal-mean.psam)
mean.pal; mean.psam; mean.diff
plot(cev.svazky.max~taxon2, data=dat)

pal.psam<-c(pal,psam) 
diffs<-vector(length=1000)

for (i in 1:1000)
{
pal.psam<-sample(pal.psam) #tohle nahodne zamicha data
rand.pal<-pal.psam[1:length(pal)]
rand.psam<-pal.psam[(1+length(pal)):length(pal.psam)]
diffs[i]<-abs(mean(rand.pal)-mean(rand.psam))
}

hist(diffs,breaks=20)
abline(v=mean.diff,col="red")
p.random=length(which(diffs>mean.diff))/1000; p.random #pravdepodobnost stanovena nasi randomizaci

anova(aov(cev.svazky.max~taxon2, data=dat)) #ANOVA test



#_____ Potreba Bonferoniho korekce _____
res<-vector(length=100)
a<-runif(100) #vygeneruje 100 nahodnych cisel od 0 do 1
for (i in 1:100)
{
m<-lm(a~runif(100))
res[i]<-anova(m)[1,5] #vypise pravdepodobnost korelace a napise ji do vektoru res
}
sort(res) #seradi hodnoty podle velikosti
length(which(res<0.05)) #vypise, kolik je hodnot mensi nebo vetsich nez zadane cislo




### CV3 ########################################################

#_____ Pocet opakovani test _____
a<-rnorm(1000,mean=20,sd=5)
as<-sample(a)
as.diff<-vector(length=1000)

for(i in (1:1000))
{as.diff[i]=abs(mean(as[1:i])-mean(as))}
plot(as.diff)

#ylim=max(c(abs(max(a)-mean(a)),abs(min(a)-mean(a)))) #největší rozdíl jakéhokoliv čísla v souboru a průměru

#a teď to celé zopakujem pro 100 takových náhodných výběrů  
lim=100 # tady si dáme maximum počtu, co bychom chtěli nakonec měřit - můžem i těch 1000, ale pak se to může počítat trochu dýl
plot(a,xlim=c(0,lim),ylim=c(0,sd(a)),type="n")
x=0
while(x<100){
  x<-x+1
  as<-sample(a)
  as.diff<-vector(length=lim)
  for(i in (1:lim)) 
    {as.diff[i]=abs(mean(as[1:i])-mean(as))}
  points(as.diff,cex=0.5)
  }

#_____ Korelace _____
setwd("G:/prednaskaMetody/Cvika/Cvika2021")
stom<-read.table("pruduchy.txt", header=T, sep = '\t')
stom<-stom[complete.cases(stom),] #odstraní z dat všechny řádky, ve kterých je NA
str(stom)

plot(l.stom~X2C,pch=21,cex=1.5,bg=grey(CV.2C/max(CV.2C)),data=stom) #neni moc normalne rozložené - obě proměnné proto zlogaritmujeme
plot(log(l.stom)~log(X2C),pch=21,cex=1.5,bg=grey(CV.2C/max(CV.2C)),data=stom)

cor.test(stom$X2C,stom$l.stom, method="pearson")
cor.test(~X2C+l.stom, method="pearson",data=stom) #alternativní zápis téhož
cor.test(~X2C+l.stom, method="spearman",data=stom) #Spearman
cor.test(~X2C+l.stom, method="kendall",data=stom) #Kendall
cor.test(~log(X2C)+log(l.stom), method="pearson",data=stom) # Pearson s logaritmovanými daty

#_____ LS (Model I) Regrese _____
m<-lm(l.stom~X2C, data=stom)
  par(mfrow=c(2,2))
  plot(m)

  anova(m); summary(m)
  par(mfrow=c(1,1))
  plot(l.stom~X2C,pch=21,cex=1.5,bg=grey(CV.2C/max(CV.2C)),data=stom)
  abline(m,col="red")


mg<-lm(log(l.stom)~log(X2C), data=stom) # s logaritmovanými daty
  par(mfrow=c(2,2))
  plot(mg)

  anova(mg); summary(mg)
  par(mfrow=c(1,1))
  plot(log(l.stom)~log(X2C),pch=21,cex=1.5,bg=grey(CV.2C/max(CV.2C)),data=stom)
  abline(mg,col="blue")

mg2<-lm(log(X2C)~log(l.stom), data=stom) # s logaritmovanými daty ale v opačném pořadí proměnných - t.j., že pruduchy jsou ta proměnná, která nemá variabilitu
  summary(mg);summary(mg2)
  #z výsledku mg2: log(X2C)=2.4955 + 1.9114*log(l.stom) si vyjadrime log(l.stom)=log(X2C)/1.9114 -2.4955/1.9114 a vykreslime to do grafu se standardní regresí abline(mg)
  abline(a=-2.4955/1.9114, b=1/1.9144,col="red")  #je z toho vidět, ža v OLS regresi na tom, co je závislá (bez variability) a nezávislá proměnná hodně závisí


mgw<-lm(log(l.stom)~log(X2C), weights=1/CV.2C,data=stom) # s logaritmovanými daty
  anova(mgw); summary(mgw)
  par(mfrow=c(1,1))
  points(log(l.stom)~log(X2C),pch=21,cex=1.5,bg=grey(CV.2C/max(CV.2C)),data=stom) #body obarvené podle hodnoty CV velikosti genomu
  abline(mgw,col="gold")
  legend("topleft", c("l.stom  predicted unweighted", "X2C predicted unweighted", "l.stom predicted weighted"), col = c("blue","red","gold"),lty=1)

  

#_____ Model II Regrese _____
library(lmodel2)
vignette("mod2user") #otevře detailní pdf s návodem co, jak a proč

mr2<-lmodel2(log(l.stom)~log(X2C), range.x = "interval",range.y = "interval",data=stom, nperm=999)
# pro spočtení RMA regrese je potřeba specifikovat charkter dat - jestli je minumium dat opravdová 0 a nemůže to být ani náhodou záporné (např. teplota ve stupních Kelvina, abundance) - pak je range "relative"; nebo jestli je minimum jen nějaká zvolená hodnota z dat, která může být ale nemusí být nutně záporná nebo nula - potom je range "interval"; pokud ten range nezadáme, spočítají se ostatní typy regreseregrese bez RMA 

mr2

par(mfrow=c(2,2))
plot(mr2,method="OLS")
plot(mr2,method="MA")
plot(mr2,method="SMA")
plot(mr2,method="RMA")


#_____ Kvantilová regrese _____
library(quantreg)
vignette("rq") #otevře detailní pdf s návodem co, jak a proč
mq<-rq(log(l.stom)~log(X2C), data=stom, tau = 0.5)
mq
summary(mq)

plot(log(l.stom)~log(X2C),pch=21,cex=1.5,data=stom)
abline(mq,col="blue")
abline(lm(log(l.stom)~log(X2C), data=stom),lty=2,col="red") # normální LS regrese pro srovnání

ts <- c(0.05,0.1,0.25,0.75,0.90,0.95)
for( i in 1:length(ts))
  {abline(rq(log(l.stom)~log(X2C), data=stom, tau=ts[i]),col="grey")}

mq.fit<-summary(rq(log(l.stom)~log(X2C), data=stom,tau=2:98/100))
plot(mq.fit,mfrow=c(1,2))



### CV4 ########################################################

setwd("G:/prednaskaMetody/Cvika/Cvika2021")
mat<-read.table("Festuca_matice.txt", header=T, sep = '\t')
str(mat)

#v matici nejsou žádná NA data. Pokud by byla dají se vyloučit, nebo nahradit průměrem daného znaku:
mat1<-na.omit(mat) #odstraní všechny řádky s NA daty
mat1<-droplevels(matn1) #někdy po odstranění dat zbydou nepoužívané hodnoty faktorů - ty je pak tímhle taky potřeba odstranit
mat2<-na.meansubst(mat) #nahradí NA průměrnou hodnotou znaku - to může a nemusí být užitečné nebo škodlivé

#___Deskriptivní statistika______________
pdf("histogramy_a_boxploty.pdf") #udělá pdf, kam nakreslí obrázky pro všechny proměnné s testem normality
for (i in names(mat)[4:length(mat)])
  {par(mfrow=c(2,1))
   hist(mat[,i],main=i,xlab="")
   legend("topright",paste("SW p = ",round(shapiro.test(mat[,i])$p.value,6)))
   boxplot(mat[,i]~taxon,main=i,data=mat)
  }     
dev.off()

### Pokračování pro pokročilejší:  ona ta nenormalita by ale mohla být jednoduše tím, že je tam jeden hodně odlišný taxon - ten by bylo potřeba odlišit a ne ten jeho rozdíl smazat. Chtělo by to udělat test pro jednotlivé skupiny zvlášť - a nebo je všechny převést na stejný průměr a pak znnovu otestovat:

#Asi se to dá udělat i jednoduššeji, ale tohle funguje: pro každý znak w se zjistí jeho celkový průměr (w.mean), pak se pro každý řádek x znaku w zjistí, co je ten řádek za taxon (tax), pak se spočítá pro daný znak w průměr všech řádků, které jsou označeny jako stejný taxon (tax.mean) a nakonec se hodnota v každé buňce vydělí průmernou hodnotou znaku dw daného taxonu (což je centrování na průměr) a převede se to zpátky na hodnoty v původních jednotách vynásobením celkovou průměrnou hodnotou znaku (w.mean) (jinak by ty menší hodnoty než průměr byly od 0 do 1 a větší od 1 do nekonečna - t.j nesymetricky). Matematické operace je potřeba dělat s původní maticí mat a výsledky psát do té nové, jinak by se v průběhu výpočtu měnilo tax.mean.

mat.centr<-mat
for (w in names(mat)[4:length(mat)]) {
   w.mean=mean(mat[,w])
   for (x in 1:length(mat$taxon)) {
     tax=mat[x,"taxon"]
     tax.mean=mean(mat[which(mat$taxon==tax),w])
     mat.centr[x,w]=w.mean*mat[x,w]/tax.mean
  }
}

#a teď se to všechno znovu vykreslí, jako jsme to dělali s neupravenými daty
pdf("histogramy_a_boxploty_centred.pdf") #udělá pdf, kam nakreslí obrázky pro všechny proměnné s testem normality pro data centrovaná na průměr skupin
for (i in names(mat.centr)[5:length(mat.centr)])
  {par(mfrow=c(2,1))
   hist(mat.centr[,i],main=i,xlab="")
   legend("topright",paste("SW p = ",round(shapiro.test(mat.centr[,i])$p.value,6)))
   boxplot(mat.centr[,i]~taxon,main=i,data=mat.centr)
  }     
dev.off()

# Podle výsledků bych nakonec logaritmoval maxlaty, maxramen, max1N, medklas, maxpochv, minskldole, maxskldole, medcevsv, maxsklbok; 0/1 znaky bych nechal jak jsou, procenta bych převedl na empirické logity s tím, že je to většinou měření ze tří listů.

logaritmovat<-c("maxlaty", "maxramen", "max1N", "medklas", "maxpochv", "maxlaty","minskldole", "maxskldole", "medcevsv", "maxsklbok") 
for (q in logaritmovat)
   {mat[,paste("ln",q,sep=".")]<-log(mat[,q])} #Tohle udělá novou zlogaritmovanou prměnou se svým vlastním názvem

mat$lgt.prerus<-log((3*mat$X.prerus+1/2)/(3-3*mat$X.prerus+1/2)) #log (pocet nastani/nenastani), kde se k poctu pricte 1/2 jako nejistota mereni - nedeli se pak nikdy nulou: třeba str. 10 tady: https://data.princeton.edu/wws509/notes/c3.pdf nebo v kinzkach Stano Pekára

#teď bych zkontroloval, jak moc ta logaritmizace pomohla (zase pomocí histogramů a boxpoltů, předpokládejme, že všude to pomohlo) a nakonec by ty redundantní proměnné vyhodil:
mat2<-mat[,!(names(mat) %in% logaritmovat)]
 
#___sumarizace a kontrola dat s načtenými MorphoTools funkcemi:
md<-read.morphodata("Festuca_matice.txt")
summary.variables<-descr.all(md)
summary.variables
summary.variables.taxa<-descr.tax(md)
summary.variables.taxa

charhist(md) #je potreba klikat pro zobrazeni grafu dalsi promenne nebo to cele ulozit do pdf
charhist(md,"mm")


#___Korelační matice znaků__________________
cor.variables.P<-cor(mat[,-(1:3)],use="pairwise.complete.obs",method="pearson")
cor.variables.P
write.table(cor.variables.P,"Pearson_correlation_matrix_variables.xls",sep="\t") #Tohle potom otevřít v Excelu (první řádek je potřeba posunout o jednu buňku vpravo) a nechat si automaticky obarvit hodnoty korelací pomocí podmíněného formátování

cor.variables.K<-cor(mat[,-(1:3)],use="pairwise.complete.obs",method="kendall") # To samé jen pro Kendalovu korelaci
write.table(cor.variables.K,"Kendall_correlation_matrix_variables.xls",sep="\t")

#to samé pomocí MrphoTools funkce cormat
cor.variables.S<-cormat(md,method="spearman")
cor.variables.S
export.res(cor.variables.S) #tohle to exportuje do clippboardu:otevřete MS Excel a do buňky A1 dejte vložit (Ctrl+V) - zkuste si to

#___Standardizace a Euklidovská vzdálenost s normálními daty__________________
mat.sd.scale<-mat
names(mat.sd.scale)
mat.sd.scale[5:37]<-lapply(mat.sd.scale[5:37],scale) #na směrodatnou odchylku

mat.range.scale<-mat
mat.range.scale[5:37]<-lapply(mat.range.scale[5:37],function(x) (x-min(x))/(max(x)-min(x))) #na rozsah od 0 do 1

E.dist.sd<-dist(mat.sd.scale[,5:37])
E.dist.sd<-as.matrix(E.dist.sd, labels=TRUE)
colnames(E.dist.sd) <- rownames(E.dist.sd) <- mat.sd.scale[,"combicislo"]
head(E.dist.sd)
write.table(E.dist.sd,"Euclidian_distance_sd_scale_transformed_data.xls",sep="\t") #Tohle potom otevřít v Excelu (první řádek je potřeba posunout o jednu buňku vpravo) a nechat si automaticky obarvit hodnoty korelací pomocí podmíněného formátování

E.dist.range<-dist(mat.range.scale[,5:37]) #analogie předchozího pro range transformovaná data
E.dist.range<-as.matrix(E.dist.range, labels=TRUE)
colnames(E.dist.range) <- rownames(E.dist.range) <- mat.range.scale[,"combicislo"]
head(E.dist.range)
write.table(E.dist.range,"Euclidian_distance_range_scaled_data.xls",sep="\t")


#___Vzdálenosti s 0/1 (AFLP) daty__________________
aflp<-read.csv("Cirsium_AFLP.csv",header=F,row.names=1)

library(vegan)
J.dist.aflp<-vegdist(aflp, method="jaccard",diag=TRUE,upper=TRUE)
write.table(as.matrix(J.dist.aflp,labels=T),"Jaccard_distance_Cirsium_AFLP_data.xls",sep="\t")


#___Vzdálenosti s ATGC (sekvenčními) daty__________________
library(ape)
rbcl.long<-read.dna("Alignment_Angiosperms_rbcL.FAS",format = "fasta")
rbcl<-rbcl.long
row.names(rbcl)<-c("Nymphaea_odorata", "Ranunculus_acris", "Geranium_divaricatum", "Poa_trivialis", "Ginkgo_biloba", "Cabomba_caroliniana", "Amborella_trichopoda", "Schisandra_sphenanthera", "Tasmannia_insipida", "Ceratophyllum_demersum", "Oryza_sativa", "Sagittaria_latifolia", "Lilium_superbum", "Platanus_occidentalis", "Nelumbo_lutea", "Spinacia oleracea", "Pisum sativum", "Aster_bakeranus")

K80.dist.rbcl<-dist.dna(rbcl, model="K80", gamma=TRUE, as.matrix=T)
K80.dist.rbcl
write.table(K80.dist.rbcl,"Kimura_1980_distance_rbcL_data.xls",sep="\t")

JC69.dist.rbcl<-dist.dna(rbcl, model="JC69", gamma=TRUE, as.matrix=T)
JC69.dist.rbcl
write.table(JC69.dist.rbcl,"Jukes_Cantor_1969_distance_rbcL_data.xls",sep="\t")



### CV5 ########################################################

setwd("G:/prednaskaMetody/Cvika/Cvika2021")

#_____PCoA (MDS) s distanční maticí Festuca
Edist.scaled<-read.table("Euclidian_distance_sd_scale_transformed_data.xls",header=T,sep="\t")
Edist.scaled<-as.dist(Edist.scaled)

mat<-read.table("Festuca_matice.txt", header=T, sep = '\t')
rownames(mat)<-mat$combicislo
tax.col<-as.factor(mat$taxon)
levels(tax.col)<-c("tomato","gold","green3","dodgerblue","darkred","violet") #barvičky pro taxony

pcoa <- cmdscale(Edist.scaled)

plot(pcoa,pch=16,cex=1.8,col=as.character(tax.col),xlab="PCoA_1",ylab="PCoA_2")
text(pcoa,mat$combicislo, cex=0.6)
legend("bottomleft", levels(as.factor(mat$taxon)), col=levels(tax.col), pch = 16,pt.cex=1.8)


#_____PCoA s distanční maticí Cirsium
Jdist.aflp<-read.table("Jaccard_distance_Cirsium_AFLP_data.xls",header=T,sep="\t")
Jdist.aflp<-as.dist(Jdist.aflp)

aflp<-read.csv("Cirsium_AFLP.csv",header=F,row.names=1)
Cirsium.sp<-gsub('[0-9]+', '', rownames(aflp))
Cirsium.sp<-as.factor(Cirsium.sp)

pcoa <- cmdscale(Jdist.aflp,k=3) #Spočtem skóre pro první tři osy: to specifikuje parametr k

plot(pcoa,pch=16,col=as.factor(Cirsium.sp),cex=1.8,xlab="PCoA_1",ylab="PCoA_2")
text(pcoa, rownames(aflp), cex=0.6)
legend("topleft", ncol=3, levels(as.factor(Cirsium.sp)), col=1:length(levels(as.factor(Cirsium.sp))), pch = 16,pt.cex=1.8)


# 3D plot třeba pomocí těchto balíků - existuje jich spousta
library("scatterplot3d")
scatterplot3d(pcoa,type="h",pch=16,cex.symbols=1.8,color=as.numeric(as.factor(Cirsium.sp)),grid=T)

library(lattice)
cloud(pcoa[,1] ~ pcoa[,2] + pcoa[,3], groups=Cirsium.sp, pch=16, cex=1.8, col=1:length(levels(as.factor(Cirsium.sp))),xlab="PCoA 1",ylab="PCoA 2",zlab="PCoA 3")
legend("topright", bg="white",levels(as.factor(Cirsium.sp)), col=1:length(levels(as.factor(Cirsium.sp))), pch = 16,pt.cex=1.8)

#_____PCA s daty + biplot

#__Načtení a nachystání dat

mat<-read.table("Festuca_matice.txt", header=T, sep = '\t')
rownames(mat)<-mat$combicislo

# Data je potřeba napřed upravit podle výsledků z minulého scvičení:
logaritmovat<-c("maxlaty", "maxramen", "max1N", "medklas", "maxpochv", "maxlaty","minskldole", "maxskldole", "medcevsv", "maxsklbok") 
for (q in logaritmovat)
   {mat[,q]<-log(mat[,q])} #ln tranformace promenych

mat$X.prerus<-log((3*mat$X.prerus+1/2)/(3-3*mat$X.prerus+1/2)) #logit transformace

write.table(mat,"Festuca_matice_transformovana.txt", sep = '\t',row.names=F) #uložíme takto upravená data pro další použití


#__PCA
# Pokud budeme chtít udělat PCA a biplot, budeme muset počítat ne z matice podobností ale ze surových dat, kde jsou informace k jednotlivým proměnným
pca<-prcomp(mat[,-(1:4)],center=T,scale=T) #tady je pouzito centrovani a standardizace dat na smerodatnou odchylku

summary(pca) #zobrazí %vysvětlené variability
head(summary(pca)) #vypíše všechny výsledky z hoto se pak dá cokoliv kreslit, automaticky to umí několik balíků
plot(pca) #vykreslí množsví vysvětlené variability jednotlvými osami do garfu - není moc užitečné bez nejaké standardizace na celkovou variabilitu. 

plot(pca$x,pch=16,cex=1.8,col=as.character(tax.col), xlab="PCA 1 (24.7 %)",ylab="PCA 2 (15.0 %)")
text(pca$x,mat$combicislo, cex=0.6)
legend("bottomleft", levels(as.factor(mat$taxon)), col=levels(tax.col), pch = 16,pt.cex=1.8)

sc=10 #kolikrát budeme chtít natáhnout šipky znaků
arrows(0,0,sc*pca$rotation[,1],sc*pca$rotation[,2],col="gray50",length=0.1)
text(sc*pca$rotation[,1],sc*pca$rotation[,2],rownames(pca$rotation),col="gray50")


#_____Broken Stick model (MacArthur 1957) + distanční matice pro skóre významnách os
source ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/scripts/NumEcolR1/evplot.R') #https://www.davidzeleny.net/anadat-r/doku.php/en:numecolr
evplot(pca$sdev^2) #to sdev z výstupu prcomp je potřeba umocnit na druhou, abychom dostali eigenvalues (variabilitu)

#_____Další programy k zobrazení PCA

#__Vegan
library(vegan) #zakladn balík pro ordinační metody pro vegeťáky; víc k tomu třeba tady: https://www.davidzeleny.net/anadat-r/doku.php/en:ordiagrams_examples
ordiplot(pca)
ordiplot(pca,display="sites")

ordiplot(pca,display="sites",type="n")
orditorp(pca,display="sites",col=as.character(tax.col)) #přidá popisky vzorků

biplot (pca) #blbě se v tom ale označují skupiny

#__ggfortify
#další tipy ke kreslení tady: https://www.davidzeleny.net/anadat-r/doku.php/en:ordiagrams_examples
#Jak co zjistit a vykreslit se dá najít třeba tady: https://bioconductor.org/packages/release/bioc/vignettes/PCAtools/inst/doc/PCAtools.html,nebo tady: https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_pca.html, nebo se to dá udělat pomocí MorphoTools funkcí (návod v manuálu k MorphoTools).

library(ggfortify) # umí kreslit pekne pca objekty; co a jak najdete tady: 
browseVignettes("ggfortify")# jaka je k tomu dokunetace
browseVignettes("ggfortify")$ggfortify #jaky mají ty vignettes jména
vignette("plot_pca") #vyvolání konkerétní vignette

autoplot(pca, data=mat, colour="taxon", size=5, loadings=T, loadings.label=T, loadings.colour="black", loadings.label.colour="black") # nakreslí objekty i proměnné

autoplot(pca, data=mat, colour="taxon", size=5, label=T,label.colour="grey50", loadings=T, loadings.label=T, loadings.colour="black",loadings.label.colour="black") # tadz to samé s popisky vzorků

#_____PCA s pomocí MorphoTools + zobrazení dalších výsledků
# Ujistíme se, že máme v R nakopírované potřebné funkce ze souboru "MorphoTools1-1.R.txt" cely obsah souboru stačí zkopírovat ctrl+c a vložit do R (ctrl+V)

# Načteme námi upravená data pro použití s MorphoTools funkcemi
md<-read.morphodata("Festuca_matice_transformovana.txt") #načte data pro použití s MorphoTools funkcemi; v prvních třech sloupcích musí být id vzorku, ve druhém populace a ve třetím taxonu.
md<-md[,-4] #ploidii (čtvrtý sloupec v datech) si vyhodíme - to není nezávislý znak ale to hlavní rozlišovací kritérium (navíc stabilní ve skupinách) 

mpca<-pca.calc(md) # Takhle by se dělala PCA se surovými daty - MorphoTools proměnné vycentruje a standardizuje na směrodatnou odchylku.

#Vzhledem k tomu, už máme PCA spočítamou z vlastní maticí, mohli bychom využít funkce MorphoTools jen k zobrazení našich pca výsledků a ne nutně museli počítat pca znovu (místo mpca bychom do dalších příkazů dávali pouze pca)

pca.samples.score<-pca.scores(mpca,md) #extrahuje skóre pro jednotlivé vzorky, to se dá potom uložit buď pomocí write. table nebo pomocí morphoToolovské fce export.res, ktera to exportuje do clippboardu - pak to vlozite do Excelu.

export.res(pca.samples.score) #zkuste si to

pca.variables.correlations<-pca.cor(mpca) #korelace proměnných s prvními 4 ordinačními osami (pca.cor(Mpca,6) by třeba spočítalo pro 6 prvních os)
pca.variables.correlations

plot.scores(pca.samples.score) #tohle mi nějak nefunguje - bylo to vidět už při načítání této funkce do R, že je tam něco blbě v kódu této funkce, ale biploty kreslit jde - možná by se to podle toho dalo opravit (domácí úkol?)
#a nebo se dá spustit funkce plot.biplot, které se nedají data k proměnným - bude se vztekat, ale obrázek nakreslí.

plot.characters(pca.variables.correlations,col="black") #pro větší písmo popisek znaků je to asi potřeba upravit přímo v posledním řádku MorphoTools funkce

plot.biplot (pca.samples.score,pca.variables.correlations) 
plot.biplot (pca.samples.score,pca.variables.correlations,scale=6, labels=T) #s delšími šipkami a popisky

#_____Porovnání PCA s daty standardizovanými na rozsah a sd
# Pokud budeme chtít udělat standardizaci na rozsah - hodně zvýrazňuje rozdíly v prezenčně absenčních znacích: To není vždycky výhodné viz příklad dole, kde mají proměnné odění hrozně velkou váhu:

mat.range.scale<-mat
mat.range.scale[5:27]<-lapply(mat.range.scale[5:27],function(x) (x-min(x))/(max(x)-min(x))) #na rozsah od 0 do 1

pcar<-prcomp(mat.range.scale[,-(1:4)],center=T) #tady je pouzito centrovani a standardizace dat na smerodatnou odchylku
pca<-prcomp(mat[,-(1:4)],center=T,scale=T) #standardní pca centrovaná a standardizovaná na sd

require(gridExtra) #gg ploty nejde kreslit do jednoho okna pomocí mfrow=c(x,y) a je na to potřeba tahle knihovna
range.pca<-autoplot(pcar, data=mat, colour="taxon", size=3, loadings=T, loadings.label=T, loadings.colour="black", loadings.label.colour="black")
sd.pca<-autoplot(pca, data=mat, colour="taxon", size=3, loadings=T, loadings.label=T, loadings.colour="black", loadings.label.colour="black")
grid.arrange(range.pca, sd.pca, ncol=2)


#_____NMDS
Edist.scaled<-read.table("Euclidian_distance_sd_scale_transformed_data.xls",header=T,sep="\t")
Edist.scaled<-as.dist(Edist.scaled)

library(vegan)
nmds<- metaMDS(Edist.scaled, k=2)
nmds3<- metaMDS(Edist.scaled, k=3)
nmds4<- metaMDS(Edist.scaled, k=4)
nmds5<- metaMDS(Edist.scaled, k=5)

par(mfrow=c(2,2))
stressplot (nmds) #porovná vzdálenosti v originální matici a v NMDS diagramu
stressplot (nmds3)
stressplot (nmds4)
stressplot (nmds5)

par(mfrow=c(1,2))
plot(nmds3$points,pch=16,cex=1.8,col=as.character(tax.col),xlab="NMDS 1",ylab="NMDS 2")
text(nmds3$points,mat$combicislo, cex=0.6)
legend("topleft", ncol=2,levels(as.factor(mat$taxon)), col=levels(tax.col), pch = 16,pt.cex=1.8)

plot(nmds3$points[,c(3,2)],pch=16,cex=1.8,col=as.character(tax.col),xlab="NMDS 3",ylab="NMDS 2")
text(nmds3$points[,c(3,2)],mat$combicislo, cex=0.6)
legend("topleft", ncol=2,levels(as.factor(mat$taxon)), col=levels(tax.col), pch = 16,pt.cex=1.8)

# NMDS pro vzdalenost vzorku z PCA pro vyznamne PCA osy
pca4.dist<-dist(pca$x[,1:4])
nmdspca<- metaMDS(pca4.dist, k=2)
stressplot (nmdspca)

par(mfrow=c(1,2))
plot(nmdspca$points,pch=16,cex=1.8,col=as.character(tax.col),xlab="NMDS 1",ylab="NMDS 2")
text(nmdspca$points,mat$combicislo, cex=0.6)
legend("topleft", ncol=2,levels(as.factor(mat$taxon)), col=levels(tax.col), pch = 16,pt.cex=1.8)

plot(pca$x,pch=16,cex=1.8,col=as.character(tax.col),xlab="PCA 1",ylab="PCA 2")
text(pca$x,mat$combicislo, cex=0.6)
legend("topleft", ncol=2,levels(as.factor(mat$taxon)), col=levels(tax.col), pch = 16,pt.cex=1.8)




### CV6 ########################################################

setwd("G:/prednaskaMetody/Cvika/Cvika2021")

#_____Kanonická diskriminační analýza pomocí Morpho Tools (balíčku "vegan")
md<-read.morphodata("Festuca_matice_transformovana.txt") #načte data pro použití s MorphoTools funkcemi; v prvních třech sloupcích musí být id vzorku, ve druhém populace a ve třetím taxonu.
md<-md[,-4] #ploidii (čtvrtý sloupec v datech) si vyhodíme - to není nezávislý znak ale to hlavní rozlišovací kritérium (navíc stabilní ve skupinách)

da<-discr.calc(md) #spočítá DA
da
str(da)


da.sum<-discr.sum(da) #vypíše summary DA modelu - vysvětlenou variabilitu os a hlavně permutačně jejich významnost
da.sum

da.coef<-discr.coef(da) #vypíše koeficienty diskriminační fuknce proměnných
da.coef

da.taxa<-discr.taxa(da) #vypíše skóre centroidů taxonů
da.taxa

da.samples<-discr.scores(da,md) #vypíše skóre vzorků

da.variables<-discr.bip(da) #vypíše skóre proměnných


da.tests<-discr.test(da) # otestuje významnost (přínos) jednotlivých proměnných pro rozlišování taxonů dvěma metodami: 1-kdyby byl znak sám v modelu, a 2- a unikátní příspěvek - t.j. kolik vysvětlí, když se přídá jako poslední ke všem ostatním


# obrazek (podobně jako u PCA)
#základní
tax.col<-as.factor(md$Taxon)
levels(tax.col)<-c("tomato","gold","green3","dodgerblue","darkred","violet") #barvičky pro taxony

DA1.perc<-round(100*summary(da)$concont$importance[2,1],1) # % vysvetlene variability prvni DA osou
DA2.perc<-round(100*summary(da)$concont$importance[2,2],1) # % vysvetlene variability druhou DA osou

plot(da.samples[,c("CCA1","CCA2")],pch=16,cex=1.8,col=as.character(tax.col), xlab=paste("DA 1 (",DA1.perc," %)",sep=""), ylab=paste("DA 2 (",DA2.perc," %)",sep=""))
text(da.samples[,c("CCA1","CCA2")],labels=da.samples[,"ID"], cex=0.6)
legend("bottomleft", levels(as.factor(da.samples$Taxon)), col=levels(tax.col), pch = 16,pt.cex=1.8)

sc=6 #kolikrát budeme chtít natáhnout šipky znaků
arrows(0,0,sc*da.variables[,"CCA1"],sc*da.variables[,"CCA2"],col="gray50",length=0.1)
text(sc*da.variables[,"CCA1"],sc*da.variables[,"CCA2"],da.variables[,"Character"],col="gray50")


# zobrazení DA pomocí MorphoTools funkcí
plot.scores(da.samples) 
point.labels(da.samples,cex=0.5)

plot.characters(da.variables,col="black",cex.lab=2) #pro větší písmo popisek znaků je to asi potřeba upravit přímo v posledním řádku MorphoTools funkce

plot.biplot (da.samples,da.variables) 
plot.biplot (da.samples,da.variables,scale=6, labels=T) #s delšími šipkami a popisky


#___DA jen pro dvě skupiny
# jak je z těch grafů vidět, chtělo by to některé taxony sloučit (otestováno to bude níž v klasifikační DA)
# my si to z cvičných důvodů uděláme na dvě skupiny (psamofilní a saxikolní), abychom si vyzkoušeli použití histogramů skóre
md2<-md
levels(md2$Taxon)<-c("sax","psam","psam","sax","sax","sax")
da2<-discr.calc(md2) #spočítá DA

da2.sum<-discr.sum(da2)
da2.coef<-discr.coef(da2)
da2.taxa<-discr.taxa(da2)
da2.samples<-discr.scores(da2,md2)
da2.variables<-discr.bip(da2)
da2.tests<-discr.test(da2)

discr.2hist(da2.samples)

discr.hist(da2.samples, legend=T)

#___Stepwise DA 
sda<-discr.step(da,p=0.01) #forward selekce proměnných (dá se zvolit forward nebo backward, nebo obě, p hodnota pro zařazení proměnné - viz https://botanika.prf.jcu.cz/systematics/morphotools/MorphoTools%20manual%201-1.pdf str. 20

sda<-discr.step(da,perm=1000, dir="both",p=0.01) #při zadávání pravděpodobnosti stojí za to brát v úvahu nějakou korekci na opakované provádění testů

#pro zobrazení je potřeba udělat novou DA, ale jen s vybranými znaky - za domácí úkol



#___Klasifikační DA pomocí Morpho Tools (balíčku "MASS")


cda<-classif.da(md) # tohle krosvalidizuje jen jedince
#cdap<-classif.da(md,crossval="pop") # tohle crossvalidizuje (počítá s vyhozením) celé populace jedinců - s našimi daty to nemá smysl počítat

cda # ve výsledném objektu je uložena pravděpodobnost klasifikace do skupin (taxonů) pro jednotlivé vzorky

classif.barplot(cda, col=levels(tax.col)) # je grafické zobrazení předchozího (pravděpodobnosti klasifikace vzorků do taxonů)
par(xpd=T) #umožní kreslit legendu mimo plot
legend("top", inset=-0.1, levels(as.factor(da.samples$Taxon)), pt.bg=levels(tax.col), pch = 22, pt.cex=1.8,ncol=length(levels(as.factor(da.samples$Taxon)))); par(xpd=F)

classif.matrix(cda) #celkový úspech klasifikace do taxonů

#classif.pmatrix(cda) #úspěch klasifikace vzorků za jednotlivé populace - může ukázat divné populace; s našimi daty tohle počítat nepůjde

#___Domácí úkol: za domácí úkol si uděláte klasifikační DA s rozlišením na psamofilní diploidy (mm+vd), saxikolní diploidy (pa) a saxikolní tetraploidy (hu+sc+st)  a zobarzíte pro to obrázek skóre vzorků a proměnných z kanonické DA



#___klasifikace neznámého vzorku
cotoje<-read.morphodata("Festuca_unclear.txt")
cotoje<-cotoje[,-4] #ploidii (čtvrtý sloupec v datech) si vyhodíme
 
clas.da<-classif.samp(cotoje,md)
clas.da
export.res(clas.da) #exportuje do clippboadu - to nejlíp vložit do Excelu a obarvit podmíněným formátováním


#___KNN
knn.select(md) #hledá nejlepší k; dá se vybrat, jestli se bude krosvalidizovat defaultně na základě vzorků (crossval="indiv") nebo populací (crossval="pop"), což by ale v našem případě nedávalo smysl

knn<-knn.classif(md,16) # klasifikuje vzorky na základě daného k (množství sousedů), s našimi daty vyšlo 13
knn

classif.matrix(knn) # zobrazí klasifikaci v tabulkové podobě jako u DA
#classif.pmatrix(knn) # to samé pro případ, že bychom krosvalidizovali na základě populací

# klasifikace neznámých vzorků
# Původní funkce knn.samp v MorphoTools obsahuje chybu - data se dělí sd zvlášť pro soubor s tréninkovými data i vzorky na určení - což nedává smysl (a hází to pak někdy chyby). Použijte proto takto upravenou funkci, kde se standardizují datový i tréninkový soubor společně):

knn.samp<-function(SAMPLE,TRAINING,K){
  require(class)
  xsamp<-na.omit(SAMPLE)
  xtrain<-na.omit(TRAINING)
  scaled<-scale(rbind(xsamp[,-(1:3)],xtrain[,-(1:3)]))
  samp<-scaled[1:dim(xsamp)[1],]
  train<-scaled[(1+dim(xsamp)[1]):dim(scaled)[1],]
  classif<-xtrain[,3]
  knn.samp<-knn(train,samp,classif,k=K,prob=T,use.all=T)
  res<-data.frame(ID=xsamp$ID,Population=xsamp$Population,Taxon=xsamp$Taxon,
                  Classif=knn.samp,Prob=attr(knn.samp,"prob"))
  res$Correct<-as.character(res$Taxon)==as.character(res$Classif)
  return(res)}


clas.knn<-knn.samp(cotoje, md, 16) # klasifikace neznámých vzorků
clas.knn
classif.matrix(clas.knn)



# Dalsi metody DA treba tady: http://www.sthda.com/english/articles/36-classification-methods-essentials/146-discriminant-analysis-essentials-in-r/



### CV7 ########################################################
setwd("G:/prednaskaMetody/Cvika/Cvika2021")
mat<-read.table("Festuca_matice_transformovana.txt", header=T, sep = '\t')

##_____dendrogramy z PCA skóre
library(vegan)
pca<-prcomp(mat[,-(1:4)],center=T,scale=T)

source ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/scripts/NumEcolR1/evplot.R')
evplot(pca$sdev^2) # Broken Stick model

PCA.dist<-dist(pca$x[,1:4])
library(dendextend)
labels(PCA.dist)<-mat$combicislo


SL<-hclust(PCA.dist,method="single")
str(SL)

CL<-hclust(PCA.dist,method="complete")
UPGMA<-hclust(PCA.dist,method="average")
WARD<-hclust(PCA.dist,method="ward.D")


plot(CL)
plot(CL, hang=-1, cex=0.7, main="Complete linkage",xlab="",sub="") # to negativní hang protáhně koncové větve až k ose

rect.hclust(CL, k=6, border="red")
rect.hclust(CL, h=7, which=c(2,7), border=c("gold","dodgerblue"))


# graficky se s tím ale pracuje líp jako s objektem typu "dendrogram" - as.dendrogram(x) nebo jako s fylogenetickým stromem - as.phylo(x) což budeme dělat později
# pěkný povídání o zobrazování stromů je třeba tady: http://www.sthda.com/english/wiki/beautiful-dendrogram-visualizations-in-r-5-must-known-methods-unsupervised-machine-learning#plot.dendrogram-function
# další věci s dendrogramy jako jednoduše měnit jednotlivé parametry nebo porovnávat dendrogramy umí třeba balíček dendextended: https://cran.r-project.org/web/packages/dendextend/vignettes/Cluster_Analysis.html

plot(as.dendrogram(CL),horiz=T,label=15,main="Complete linkage",xlab="Distance")

nodePar <- list(lab.cex = 0.6, pch = c(16, 15), col = "blue") #pomocí nodePar a edgePar se dájí měnit vlastnosti uzlů a větví a dá se na ně kreslit; to NA je v tom pch proto, aby se to nekreslilo na uzel ale na konec větve

plot(as.dendrogram(CL), horiz=T, main="Complete linkage", xlab="Distance", nodePar=nodePar)
plot(as.dendrogram(CL), horiz=T, main="Complete linkage", xlab="Distance", nodePar=nodePar, ylim=c(1,20)) # pomocí xlim a ylim se dá koukat jen na část dendrogramu - tady prvních 20 vzorků

# minimum evolution stromy
#instal.package("ape")
library(ape)
NJ<-nj(PCA.dist) # výsledkem bude strom ve formátu phylo
plot(NJ) # strom není ultrametrický, ani zakořeněný

NJr<-root(NJ, node=2+length(NJ$tip.label),resolve.root=T) # zakořenění na uzel, který bude spojovat nejvzdálenější vzorky v matici
NJu<-chronos(NJr, lambda=0) # ultrametrizace
NJd<-as.hclust(NJu)
plot(as.dendrogram(NJd), horiz=T, main="Neighbour joining tree", xlab="Distance")

LS<-fastme.bal(PCA.dist)
#LS<-fastme.ols(PCA.dist)
#library(phytools)
#LS<-optim.phylo.ls(PCA.dist,NJ) # tohle je přesně podle Felsensteina 2004, ale je to pomalý, jako startovací strom použijeme NJ strom

LS$edge.length[LS$edge.length<0] <- 0 # musíme se zbavit případných neativních délek větví
LSr<-root(LS, node=2+length(LS$tip.label),resolve.root=T) # zakořenění na uzel, který bude spojovat nejvzdálenější vzorky v matici
LSu<-chronos(LSr, lambda=0) # ultrametrizace
LSd<-as.hclust(LSu)
plot(as.dendrogram(LSd), horiz=T, main="Least square tree", xlab="Distance")

# Změna parametrů dendrogramů pomocí balíčku dendextent; ten používá řetězení příkazů pomocí %>%, na příkazy,kterými se co mění je potřeba kouknout do vignette
#instal.package("dendextend")
library(dendextend)
browseVignettes("dendextend") #koukněte hlavně na introduction 

tax.col<-as.factor(mat$taxon)
levels(tax.col)<-c("tomato","gold","green3","dodgerblue","darkred","violet") #barvičky pro taxony
tax.col<-as.character(tax.col)
names(tax.col)<-mat$combicislo

dend.CL<-CL %>% as.dendrogram %>% set("labels_col",tax.col[CL$order]) # barvičky je potřeba dávat v pořadí, v jakém jsou vzorky v dendrogramu - tahle informace je uložena ve struktuře dendrogramu v $order
plot(dend.CL)

CL %>% as.dendrogram %>%
set("labels_col",tax.col[CL$order]) %>%
set("leaves_pch",15) %>%
set("leaves_col",tax.col[CL$order]) %>%
set("labels_cex",0.6) %>%
plot(horiz=T, main="Complete linkage", xlab="Distance")

legend("topleft", levels(as.factor(mat$taxon)), col=c("tomato","gold","green3","dodgerblue","darkred","violet"), pch = 15,pt.cex=1.8)



par(mfrow=c(2,3))
SL %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(SL))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(SL))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Single linkage", xlab="Distance")
CL %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(CL))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(CL))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Complete linkage", xlab="Distance")
WARD %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(WARD))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(WARD))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="WARD", xlab="Distance")
UPGMA %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(UPGMA))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(UPGMA))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="UPGMA", xlab="Distance")
NJd %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(NJd))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(NJd))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Neighbour joining tree", xlab="Distance")
LSd %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(LSd))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(LSd))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Least square tree", xlab="Distance")


#___dendrogramy z vlastní matice vzdáleností - je vidět, že tohle je esměshorší než práce s daty "vyčištěnými" v PCA
Edist.scaled<-read.table("Euclidian_distance_sd_scale_transformed_data.xls",header=T,sep="\t")
Edist.scaled<-as.dist(Edist.scaled)

#Edist.scaled<-read.table("Euclidian_distance_range_scaled_data.xls",header=T,sep="\t")
#Edist.scaled<-as.dist(Edist.scaled)


E.SL<-hclust(Edist.scaled,method="single")
E.CL<-hclust(Edist.scaled,method="complete")
E.UPGMA<-hclust(Edist.scaled,method="average")
E.WARD<-hclust(Edist.scaled,method="ward.D")

E.NJ<-nj(Edist.scaled) # výsledkem bude strom ve formátu phylo
E.NJr<-root(E.NJ, node=2+length(E.NJ$tip.label),resolve.root=T) # zakořenění na uzel, který bude spojovat nejvzdálenější vzorky v matici
E.NJu<-chronos(E.NJr, lambda=0) # ultrametrizace
E.NJd<-as.hclust(E.NJu)
plot(as.dendrogram(E.NJd), horiz=T, main="Neighbour joining tree", xlab="Distance")

E.LS<-fastme.bal(Edist.scaled)
E.LS$edge.length[E.LS$edge.length<0] <- 0 # musíme se zbavit případných neativních délek větví
E.LSr<-root(E.LS, node=2+length(E.LS$tip.label),resolve.root=T) # zakořenění na uzel, který bude spojovat nejvzdálenější vzorky v matici
E.LSu<-chronos(E.LSr, lambda=0) # ultrametrizace
E.LSd<-as.hclust(E.LSu)
plot(as.dendrogram(E.LSd), horiz=T, main="Least square tree", xlab="Distance")

par(mfrow=c(2,3))
E.SL %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.SL))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.SL))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Single linkage", xlab="Distance")
E.CL %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.CL))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.CL))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Complete linkage", xlab="Distance")
E.WARD %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.WARD))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.WARD))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="WARD", xlab="Distance")
E.UPGMA %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.UPGMA))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.UPGMA))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="UPGMA", xlab="Distance")
E.NJd %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.NJd))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.NJd))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Neighbour joining tree", xlab="Distance")
E.LSd %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.LSd))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.LSd))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Least square tree", xlab="Distance")



# Pro zajímavost klastrovka pomocí MorphoTools
md<-read.morphodata("Festuca_matice_transformovana.txt") #načte data pro použití s MorphoTools funkcemi; v prvních třech sloupcích musí být id vzorku, ve druhém populace a ve třetím taxonu.
md<-md[,-4] #ploidii (čtvrtý sloupec v datech) si vyhodíme - to není nezávislý znak ale to hlavní rozlišovací kritérium (navíc stabilní ve skupinách) 

SL.MT<-clust(md,method="single")
CL.MT<-clust(md,method="complete")
UPGMA.MT<-clust(md,method="average")
WARD.MT<-clust(md,method="ward")

plot(as.dendrogram(CL.MT), horiz=T, main="Complete linkage using MorphoTools", xlab="Distance")



###___Porovnání pozic vzorků
dends0<-dendlist(as.dendrogram(UPGMA),as.dendrogram(WARD))
tanglegram(dends0)
tanglegram(dends0, fast=T)

# do dendlist můžeme zadat i už jinak naformátované dendrogramy
SLd<-SL %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(SL))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(SL))]) %>% set("labels_cex",0.6)
CLd<-CL %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(CL))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(CL))]) %>% set("labels_cex",0.6)
WARDd<-WARD %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(WARD))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(WARD))]) %>% set("labels_cex",0.6)
UPGMAd<-UPGMA %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(UPGMA))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(UPGMA))]) %>% set("labels_cex",0.6) 
NJdd<-NJd %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(NJd))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(NJd))]) %>% set("labels_cex",0.6)
LSdd<-LSd %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(LSd))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(LSd))]) %>% set("labels_cex",0.6)

E.SLd<- E.SL %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.SL))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.SL))]) %>% set("labels_cex",0.6)
E.CLd<- E.CL %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.CL))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.CL))]) %>% set("labels_cex",0.6)
E.WARDd<- E.WARD %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.WARD))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.WARD))]) %>% set("labels_cex",0.6)
E.UPGMAd<- E.UPGMA %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.UPGMA))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.UPGMA))]) %>% set("labels_cex",0.6) 
E.NJdd<-E.NJd %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.NJd))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.NJd))]) %>% set("labels_cex",0.6)
E.LSdd<-E.LSd %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(E.LSd))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(E.LSd))]) %>% set("labels_cex",0.6)



dends<-dendlist(SLd,CLd,UPGMAd,WARDd,NJdd,LSdd,E.SLd,E.CLd,E.UPGMAd,E.WARDd,E.NJdd,E.LSdd)
names(dends)<-c("SL - PCA","CL - PCA","UPGMA - PCA","WARD - PCA","NJ - PCA","LS - PCA","SL - dist","CL - dist","UPGMA - dist","WARD - dist","NJ - dist","LS - dist")

# různá data
dends  %>% tanglegram(which=c(1,7), fast=T)
dends  %>% tanglegram(which=c(2,8), fast=T)
dends  %>% tanglegram(which=c(3,9), fast=T)
dends  %>% tanglegram(which=c(4,10), fast=T)
dends  %>% tanglegram(which=c(5,11), fast=T)
dends  %>% tanglegram(which=c(6,12), fast=T)

# různé klastrovací
dends  %>% tanglegram(which=c(3,4), fast=T)
dends  %>% tanglegram(which=c(3,5), fast=T)
dends  %>% tanglegram(which=c(3,6), fast=T)
dends  %>% tanglegram(which=c(5,6), fast=T)



###___korelace dendrogramů
dendlist.corP<-cor.dendlist(dends)
dendlist.corP
dendlist.corS<-cor.dendlist(dends, method_coef="spearman")
dendlist.corC<-cor.dendlist(dends, method="common_nodes")


#instal.package("corrplot")
library(corrplot)
par(mfrow=c(1,3))
corrplot(dendlist.corP, "pie", "lower")
corrplot(dendlist.corS, "pie", "lower")
corrplot(dendlist.corC, "pie", "lower")



###___Kofenetické korelace
SL.dist<-cophenetic(SL)
CL.dist<-cophenetic(CL)
WARD.dist<-cophenetic(WARD)
UPGMA.dist<-cophenetic(UPGMA)
NJ.dist<-cophenetic(NJd)
LS.dist<-cophenetic(LSd)

#PCA.dist<-dist(pca$x[,1:4]) # mám pocit, že se někde po tom použití dendextend něco uděje s tou maticí a pak vychází ty konefetický korelace plně nesmyslně - v tom případě je to lepší spočítat znovu načisto
#labels(PCA.dist)<-mat$combicislo

cor(PCA.dist,SL.dist)
cor(PCA.dist,SL.dist,method="spearman")
plot(PCA.dist,SL.dist, main=paste("SL, copehnetic cor= ",round(cor(PCA.dist,SL.dist),3)))


par(mfrow=c(2,3))
plot(PCA.dist,SL.dist, main=paste("SL, copehnetic cor= ",round(cor(PCA.dist,SL.dist),3)))
plot(PCA.dist,CL.dist, main=paste("CL, copehnetic cor= ",round(cor(PCA.dist,CL.dist),3)))
plot(PCA.dist,UPGMA.dist, main=paste("UPGMA, copehnetic cor= ",round(cor(PCA.dist,UPGMA.dist),3)))
plot(PCA.dist,WARD.dist, main=paste("WARD, copehnetic cor= ",round(cor(PCA.dist,WARD.dist),3)))
plot(PCA.dist,NJ.dist, main=paste("NJ tree, copehnetic cor= ",round(cor(PCA.dist,NJ.dist),3)))
plot(PCA.dist,LS.dist, main=paste("LS tree, copehnetic cor= ",round(cor(PCA.dist,LS.dist),3)))
par(mfrow=c(1,1))


###___Bootstrap
# na to budeme potřebovat originální data, abychom mohli vybírat (bootstrapovat) znaky
#mat<-read.table("Festuca_matice_transformovana.txt", header=T, sep = '\t')
mat2<-mat[,-(1:4)]
rownames(mat2)<-mat$combicislo
mat2<-scale(mat2) # standardizace na směrodatnou odchylku

#UPGMA.orig<-as.phylo(hclust(dist(mat2),method="average")) #uděláme potom samostatnou funkcí

fun<-function(x) as.phylo(hclust(dist(x), "average")) # funkce na dělání stromu (pomocí UPGMA)
Nb=100 # kolik bootstrap cyklů budeme počítat, na vážno samozřejmě víc
UPGMA.orig<-fun(mat2) # originální strom ze surových dat
UPGMA.bstrees <- boot.phylo(UPGMA.orig, mat2, fun, B=Nb, trees=T)$trees   # vygeneruje B bootstrapovaných stromů
#clad <- prop.clades(UPGMA.orig, UPGMA.bstrees, rooted = TRUE) # proporční zastoupení jednotlivých kládů:
boot <- prop.clades(UPGMA.orig, UPGMA.bstrees) ## get proportions of each bipartition:
boot[is.na(boot)] <- 0
clad[is.na(clad)] <- 0

par(mar=c(0,0,2,0))
plot(UPGMA.orig, cex=0.6, main="UPGMA Festuca",tip.color=tax.col[UPGMA.orig$tip.label])
#nodelabels(boot,frame="none",bg="transparent")
nodelabels(round(100*boot/Nb,0),frame="circle",bg=rgb((Nb-boot)/Nb,0,0),cex=0.6, col="white")

# tohle jde udělat i pro ty výsledky z PCA, ale do funkce by se muselo dát i udělání PCA:
mat3<-mat[,-(1:4)]
rownames(mat3)<-mat$combicislo

fun<-function(x) as.phylo(hclust(dist(prcomp(x,center=T,scale=T)$x[,1:4]), "average")) # funkce na dělání stromu (pomocí UPGMA metody)
Nb=327 # kolik bootstrap cyklů budeme počítat, na vážno samozřejmě víc
UPGMA.orig<-fun(mat3) # originální strom z PCA skóre prvních 4 os
UPGMA.bstrees <- boot.phylo(UPGMA.orig, mat3, fun, B=Nb, trees=T)$trees   # vygeneruje B bootstrapovaných stromů
#clad <- prop.clades(UPGMA.orig, UPGMA.bstrees, rooted = TRUE) # proporční zastoupení jednotlivých kládů:
boot <- prop.clades(UPGMA.orig, UPGMA.bstrees) ## get proportions of each bipartition:
boot[is.na(boot)] <- 0
clad[is.na(clad)] <- 0

par(mar=c(0,0,2,0))
plot(UPGMA.orig, cex=0.6, main="UPGMA PCA scores Festuca",tip.color=tax.col[UPGMA.orig$tip.label])
#nodelabels(boot,frame="none",bg="transparent")
nodelabels(round(100*boot/Nb,0),frame="circle",bg=rgb((Nb-boot)/Nb,0,0),cex=0.6, col="white")


###___Začlenění minimum spanning stromu do ordinace
library(vegan)
MS<-spantree(PCA.dist)

plot(MS,type="t",labels=mat$combicislo)
plot(as.hclust(MS))
plot(as.dendrogram(as.hclust(MS)),hang=-1, cex=0.7, main="Minimum_spanning",xlab="",sub="")

MSd<-as.hclust(MS)
MSd %>% as.dendrogram %>% set("labels_col",tax.col[labels(as.dendrogram(MSd))]) %>% set("leaves_pch",15) %>% set("leaves_col",tax.col[labels(as.dendrogram(MSd))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Minimum spannin tree", xlab="Distance")

plot(MS, pca, pch=16,cex=1.8,col=as.character(tax.col))
text(pca$x,mat$combicislo, cex=0.6)



###___K means clustering
KM<-kmeans(mat2, centers=2, nstart=50) # těch tsratů je potřeba dát víc - aspoň 25, aby se výsledky opakovaly, defaultně je jen 1.
str(KM)

# optimální počet klastrů
#install.packages("factoextra")
library(factoextra)
fviz_nbclust(mat2, kmeans, method = "silhouette")
fviz_nbclust(mat2, kmeans, method = "wss")

#install.packages("ClusterR")
library(ClusterR)
Optimal_Clusters_KMeans(mat2, max_clusters = 10, criterion = "distortion_fK", plot_clusters = T)

# obrázek pro surová, sd stadardizovaná data a pro PCA skóre
par(mfrow=c(1,3))
KM<-kmeans(mat[,-(1:4)], centers=3, nstart=50)
pca<-prcomp(mat[,-(1:4)],center=T,scale=T)
plot(pca$x,pch=16,cex=2.8,col=1+KM$cluster, xlab="PCA 1",ylab="PCA 2",main="K means surová data")
text(pca$x,mat$combicislo, cex=0.6)

KM.sd<-kmeans(mat2, centers=3, nstart=50)
plot(pca$x,pch=16,cex=2.8,col=1+KM.sd$cluster, xlab="PCA 1",ylab="PCA 2",main="K means sd standardizovaná data")
text(pca$x,mat$combicislo, cex=0.6)

KM.PCA<-kmeans(pca$x[,1:4], centers=3, nstart=50)
plot(pca$x,pch=16,cex=2.8,col=1+KM.PCA$cluster, xlab="PCA 1",ylab="PCA 2",main="K means první 4 PCA osy")
text(pca$x,mat$combicislo, cex=0.6)



###___ AFLP Cirsium
Jdist.aflp<-read.table("Jaccard_distance_Cirsium_AFLP_data.xls",header=T,sep="\t")
Jdist.aflp<-as.dist(Jdist.aflp)

aflp<-read.csv("Cirsium_AFLP.csv",header=F,row.names=1)
Cirsium.sp<-gsub('[0-9]+', '', rownames(aflp))
Cirsium.sp<-as.factor(Cirsium.sp)

#install.packages("viridis")
#library(viridis)
Cirs.col<-as.factor(Cirsium.sp)
levels(Cirs.col)<-  rainbow(length(levels(Cirs.col))) # viridis(length(levels(Cirs.col)))    
Cirs.col<-as.character(Cirs.col)
names(Cirs.col)<-rownames(aflp)


C.SL<-hclust(Jdist.aflp,method="single")
C.CL<-hclust(Jdist.aflp,method="complete")
C.UPGMA<-hclust(Jdist.aflp,method="average")
C.WARD<-hclust(Jdist.aflp,method="ward.D")


C.NJ<-nj(Jdist.aflp) # výsledkem bude strom ve formátu phylo
C.NJ$edge.length[C.NJ$edge.length<0] <- 0 # musíme se zbavit případných neativních délek větví
C.NJr<-root(C.NJ, "Sylibum48",resolve.root=T) # zakořenění na uzel, který bude spojovat nejvzdálenější vzorky v matici
C.NJu<-chronos(C.NJr, lambda=0) # ultrametrizace
C.NJd<-as.hclust(C.NJu)

C.LS<-fastme.bal(Jdist.aflp)
#C.LS<-fastme.ols(Jdist.aflp)
#library(phytools)
#C.LS<-optim.phylo.ls(Jdist.aflp,C.NJ) # tohle je přesně podle Felsensteina 2004, ale je to pomalý, jako startovací strom použijeme NJ strom
C.LS$edge.length[C.LS$edge.length<0] <- 0 # musíme se zbavit případných neativních délek větví
C.LSr<-root(C.LS, "Sylibum48",resolve.root=T) # zakořenění na uzel, který bude spojovat nejvzdálenější vzorky v matici
C.LSu<-chronos(C.LSr, lambda=0) # ultrametrizace
C.LSd<-as.hclust(C.LSu)

par(mfrow=c(2,3))
C.SL %>% as.dendrogram %>% set("labels_col",Cirs.col[labels(as.dendrogram(C.SL))]) %>% set("leaves_pch",15) %>% set("leaves_col",Cirs.col[labels(as.dendrogram(C.SL))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Single linkage", xlab="Distance")
C.CL %>% as.dendrogram %>% set("labels_col",Cirs.col[labels(as.dendrogram(C.CL))]) %>% set("leaves_pch",15) %>% set("leaves_col",Cirs.col[labels(as.dendrogram(C.CL))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Complete linkage", xlab="Distance")
C.WARD %>% as.dendrogram %>% set("labels_col",Cirs.col[labels(as.dendrogram(C.WARD))]) %>% set("leaves_pch",15) %>% set("leaves_col",Cirs.col[labels(as.dendrogram(C.WARD))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="WARD", xlab="Distance")
C.UPGMA %>% as.dendrogram %>% set("labels_col",Cirs.col[labels(as.dendrogram(C.UPGMA))]) %>% set("leaves_pch",15) %>% set("leaves_col",Cirs.col[labels(as.dendrogram(C.UPGMA))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="UPGMA", xlab="Distance")
C.NJd %>% as.dendrogram %>% set("labels_col",Cirs.col[labels(as.dendrogram(C.NJd))]) %>% set("leaves_pch",15) %>% set("leaves_col",Cirs.col[labels(as.dendrogram(C.NJd))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Neighbour joining tree", xlab="Distance")
C.LSd %>% as.dendrogram %>% set("labels_col",Cirs.col[labels(as.dendrogram(C.LSd))]) %>% set("leaves_pch",15) %>% set("leaves_col",Cirs.col[labels(as.dendrogram(C.LSd))]) %>% set("labels_cex",0.6) %>% plot(horiz=T, main="Least square tree", xlab="Distance")



###___ rbcL sekvenční data
Jdist.rbcl<-read.table("Jukes_Cantor_1969_distance_rbcL_data.xls",header=T,sep="\t")
Jdist.rbcl<-as.dist(Jdist.rbcl)

R.SL<-hclust(Jdist.rbcl,method="single")
R.CL<-hclust(Jdist.rbcl,method="complete")
R.UPGMA<-hclust(Jdist.rbcl,method="average")
R.WARD<-hclust(Jdist.rbcl,method="ward.D")


R.NJ<-nj(Jdist.rbcl) # výsledkem bude strom ve formátu phylo
R.NJ$edge.length[R.NJ$edge.length<0] <- 0 # musíme se zbavit případných neativních délek větví
R.NJr<-root(R.NJ, "Ginkgo_biloba",resolve.root=T) # zakořenění na uzel, který bude spojovat nejvzdálenější vzorky v matici
R.NJu<-chronos(R.NJr, lambda=0) # ultrametrizace
R.NJd<-as.hclust(R.NJu)

R.LS<-fastme.bal(Jdist.rbcl)
#R.LS<-fastme.ols(Jdist.rbcl)
#library(phytools)
#R.LS<-optim.phylo.ls(Jdist.rbcl,R.NJ) # tohle je přesně podle Felsensteina 2004, ale je to pomalý, jako startovací strom použijeme NJ strom
R.LS$edge.length[R.LS$edge.length<0] <- 0 # musíme se zbavit případných neativních délek větví
R.LSr<-root(R.LS, "Ginkgo_biloba",resolve.root=T) # zakořenění na uzel, který bude spojovat nejvzdálenější vzorky v matici
R.LSu<-chronos(R.LSr, lambda=0) # ultrametrizace
R.LSd<-as.hclust(R.LSu)

par(mfrow=c(2,3),no.margin=T)
R.SL %>% as.dendrogram %>% plot(horiz=T, main="Single linkage", xlab="Distance")
R.CL %>% as.dendrogram %>% plot(horiz=T, main="Complete linkage", xlab="Distance")
R.WARD %>% as.dendrogram %>% plot(horiz=T, main="WARD", xlab="Distance")
R.UPGMA %>% as.dendrogram %>% plot(horiz=T, main="UPGMA", xlab="Distance")
R.NJd %>% as.dendrogram %>% plot(horiz=T, main="Neighbour joining tree", xlab="Distance")
R.LSd %>% as.dendrogram %>% plot(horiz=T, main="Least square tree", xlab="Distance")

# Bootstrap rbcl - za domácí úkol (aspoň NJ a LS strom)



### CV8 ########################################################
setwd("G:/prednaskaMetody/Cvika/Cvika2021")

# Multiple alignment
# dá se dělat v mnoha balících jako DECIPHER, ips, msa (viz. https://stackoverflow.com/questions/4497747/how-to-perform-basic-multiple-sequence-alignments-in-r?utm_medium=organic&utm_source=google_rich_qa&utm_campaign=google_rich_qa; https://gtpb.github.io/MEVR16/align/exercise1.html). My zkusíme balík ips, který pracuje přímo s běžně používanými ne R balíky, které je ale potřeba si vesměs naistalovat


#_ClustalW a Muscle 
if (!requireNamespace("BiocManager", quietly = TRUE)) # balik msa neni jednoduse dostupny na beznych R cloudech / musi se instalovat ze stranek Bioconductor 
    install.packages("BiocManager")

BiocManager::install("msa")

library(msa)
#browseVignettes("msa")

rbcl.raw<-readDNAStringSet("Angiosperms_rbcL.FASTA")
names(rbcl.raw)<-c("Nymphaea_odorata", "Ranunculus_acris", "Geranium_divaricatum", "Poa_trivialis", "Ginkgo_biloba", "Cabomba_caroliniana", "Amborella_trichopoda", "Schisandra_sphenanthera", "Tasmannia_insipida", "Ceratophyllum_demersum", "Oryza_sativa", "Sagittaria_latifolia", "Lilium_superbum", "Platanus_occidentalis", "Nelumbo_lutea", "Spinacia oleracea", "Pisum sativum", "Aster_bakeranus")

rbcl.clustalW <- msa(rbcl.raw)
rbcl.clustalW
print(rbcl.clustalW, show="complete") #tohle by teď bylo potřeba ořezat - tím se zabývat nebudeme

rbcl.muscle <- msa(rbcl.raw, method="Muscle")
rbcl.muscle
print(rbcl.muscle, show="complete") #tohle by teď bylo potřeba ořezat - tím se zabývat nebudeme


#_ Mafft a PRANK se mi nepodařilo rozchodit - nejde mi nějak zadat cestu k patřičným programům, co to počítají mimo R

install.packages("ips")
library(ips) # interface pro R na ostatni fylogeneticke programy, nacte v ramci toho i ape pro manipulaci se stomy

# MAFFT - je potřeba stáhnout a naistalovat si program mafft a někam si poznačit cestu k programu - https://mafft.cbrc.jp/alignment/software/windows_without_cygwin.html

#rbcl.raw <- read.fas("Alignment_Angiosperms_rbcL.FAS") #alternativni nacteni dat pomoci funkce v baliku ips
rbcl.mafft <- mafft(rbcl.raw,path="G:/prednaskaMetody/Cvika/Cvika2021/mafft-win/usr/bin/mafft")
write.fas(rbcl.mafft,"rbcl_mafft.fas")

# PRANK - je potřeba stáhnout a naistalovat si program mafft a někam si poznačit cestu k programu - http://wasabiapp.org/download/prank/

rbcl.prank <- prank(rbcl.raw,path="G:/prednaskaMetody/Cvika/Cvika2021/prank.cygwin/prank/bin/prank")
write.fas(rbcl.prank,"rbcl_prank.fas")



### Fylogenetické stromy______
#něco umí ape a to, co už jsme si zkoušeli, nejlíp ale v balíku phangorn; ostatní populární programy RAxML, MrBayes si umí v R volat balík ips
# pro jednoduchost si načteme si už ořezaný alignment

library(ape)
rbcl.long<-read.dna("Alignment_Angiosperms_rbcL.FAS",format = "fasta")
rbcl<-rbcl.long
row.names(rbcl)<-c("Nymphaea_odorata", "Ranunculus_acris", "Geranium_divaricatum", "Poa_trivialis", "Ginkgo_biloba", "Cabomba_caroliniana", "Amborella_trichopoda", "Schisandra_sphenanthera", "Tasmannia_insipida", "Ceratophyllum_demersum", "Oryza_sativa", "Sagittaria_latifolia", "Lilium_superbum", "Platanus_occidentalis", "Nelumbo_lutea", "Spinacia oleracea", "Pisum sativum", "Aster_bakeranus")

#install.packages("phangorn")
library(phangorn) #asi hlavní čistě R balíček na konstrukci fylogenetických stromů
browseVignettes("phangorn")# jaka je k tomu dokunetace

rbcl<-as.phyDat(rbcl) #phangorn potřebuje ve vlastním phyDat formátu

# jednoduché distanční stromy, co budou sloužit jako základ pro složitější fylogenetické metody
rbcl.dist <- dist.ml(rbcl)

rbcl.UPGMA <- upgma(rbcl.dist)
rbcl.UPGMAr<-root(rbcl.UPGMA, "Ginkgo_biloba",resolve.root=T)
plot(rbcl.UPGMAr)
parsimony(rbcl.UPGMAr, rbcl)  # jak je strom dlouhý

rbcl.NJ <- NJ(rbcl.dist)

rbcl.NJr<-root(rbcl.NJ, "Ginkgo_biloba",resolve.root=T)
plot(rbcl.NJr)
parsimony(rbcl.NJr, rbcl) # tenhle je kratší, ten použijem dál


# Maximum parsimonie____
rbcl.bab<-bab(rbcl) #čistě jen branch and bound algoritmus - počítá se dlouho a má smysl jen do velmi malého počtu taxonů tak do 15-20
rbcl.babr<-root(rbcl.bab, "Ginkgo_biloba",resolve.root=T)
parsimony(rbcl.babr, rbcl)

rbcl.parsNJ <- optim.parsimony(rbcl.NJr, rbcl)
rbcl.parsUPGMA <- optim.parsimony(rbcl.UPGMAr, rbcl)

rbcl.ratchet <- pratchet(rbcl, trace=0)
rbcl.ratchetr<-root(rbcl.ratchet, "Ginkgo_biloba",resolve.root=T)

parsimony(c(rbcl.UPGMAr, rbcl.NJr, rbcl.parsNJ, rbcl.parsUPGMA, rbcl.ratchetr), rbcl)

plot(rbcl.ratchetr)

boot.rbcl.ratchetr <- bootstrap.phyDat(rbcl, FUN=function(x) pratchet(rbcl, trace=0), bs=100)
boot.rbcl.ratchetr.tree<-plotBS(rbcl.ratchetr, boot.rbcl.ratchetr, p=0, "phylogram")

boot.rbcl.ratchetr.tree; add.scale.bar()
write.tree(boot.rbcl.ratchetr.tree)



# Maximum likelihood____
fit<-pml(rbcl.NJ, rbcl) # vypočte likelihood konkerétního stromu
fit

rbcl.JC <- optim.pml(fit, TRUE) # tohle optimalizuje délky větví, defaultně pomocí Jukes-Cantor modelu
rbcl.JC   # kdybychom chtěli jen hodnotu log likelighood, stačí: logLik(rbcl.JC)


model<-modelTest(rbcl) # tady tím se dají spočítat orientačně (bez branch rearrangements) všechny možné modely, což může někdy chvíli trvat
model
model[order(model$AIC),] #seřazeno podle AIC

rbcl.GTR <- update(fit, k=4, inv=0.2) #update modelu z Jukes-Cantor na maximální model e GTR + ?(4) + I
rbcl.GTR <- optim.pml(rbcl.GTR, model="GTR", optInv=TRUE, optGamma=TRUE, rearrangement = "NNI", control = pml.control(trace = 0))
rbcl.GTR

rbcl.GTR2 <- optim.pml(rbcl.GTR, model="GTR", optInv=TRUE, optGamma=TRUE, rearrangement = "stochastic", control = pml.control(trace = 0)) # s jinou verzí rearrangement, která prý líp uniká z lokálních maxim
rbcl.GTR2

# ty stromy můžeme mezi sebou porovnat pomocí AIC a anova
AIC(rbcl.JC, rbcl.GTR, rbcl.GTR2) 
anova(rbcl.GTR, rbcl.GTR2) 

boot.rbcl.GTR2<-bootstrap.pml(rbcl.GTR2, bs=100, optNni=TRUE, control = pml.control(trace = 0))

plotBS(midpoint(rbcl.GTR2$tree), boot.rbcl.GTR2, p = 0, type="p")

net.boot.rbcl.GTR2<- consensusNet(boot.rbcl.GTR2, p=0.2)  #ty všechny bootstrapované stromy se schovávají a dá se na ně kouknout pomocí konsenzuální sítě - jetsli tam není nějaký zásadní konflikt
plot(net.boot.rbcl.GTR2, "2D", show.edge.label=TRUE)


### Maticove testy________
#__Data
matk<-read.table("Festuca_matice_kultivace.txt", header=T, sep = '\t')
matk<-matk[order(matk$mereno, matk$trs),]
str(matk)


logaritmovat<-c("maxlaty", "maxramen", "max1N", "medklas", "mednrklas","maxpochv", "maxlaty", "maxskldole", "medcevsv") 
for (q in logaritmovat)
   {matk[,q]<-log(matk[,q])} #ln tranformace promenych

# matk$X.prerus<-log((3*matk$X.prerus+1/2)/(3-3*matk$X.prerus+1/2)) #logit transformace / asi ted nebudem delat

cult<-matk[which(matk$mereno=="cultivation"),]
rownames(cult)<-cult$trs

wild<-matk[which(matk$mereno=="wild"),]
rownames(wild)<-wild$trs


#___PCA
pca.wild<-prcomp(wild[,-(1:3)],center=T,scale=T) #je pouzito centrovani a standardizace dat na smerodatnou odchylku
pca.cult<-prcomp(cult[,-(1:3)],center=T,scale=T)

par(mfrow=c(1,2))
plot(pca.cult$x,pch=16,cex=0.5, xlab="PCA 1",ylab="PCA 2", main="Cultivation")
text(pca.cult$x,cult$trs, cex=0.6)

plot(pca.wild$x,pch=16,cex=0.5, xlab="PCA 1",ylab="PCA 2", main="Wild")
text(pca.wild$x,wild$trs, cex=0.6)
par(mfrow=c(1,1))


#___Mantel test & Protest
library(vegan)
mantel(dist(pca.wild$x), dist(pca.cult$x), permutations=999)

proc<-procrustes(dist(pca.wild$x),dist(pca.cult$x))
protest(dist(pca.wild$x),dist(pca.cult$x), permutations = how(nperm = 999))

plot(proc)
summary(proc)

proc.resid<-residuals(proc)
proc.resid

sort(proc.resid, decreasing=T)

plot(hclust(dist(residuals(proc)),method="single"))
plot(ladderize(as.phylo(as.hclust(spantree(dist(residuals(proc)))))),cex=0.6)

#___Které znaky se mění
# párový t test
t.test(medklas ~ mereno, data = matk, paired = TRUE) # to samé musí mimochodem vyjít z mixed effekt modelu: summary(lme(medklas ~ mereno, data = matk, random=~1|

#pro vsecny znaky
res<-data.frame(rep(NA,16),rep(NA,16),row.names=names(matk)[4:19])
names(res)<-c("t.value","p.value")

for (i in names(matk)[4:19]){
ttt<-t.test(matk[,i] ~ matk[,"mereno"], paired = TRUE)
res[i,"t.value"]<-ttt$statistic
res[i,"p.value"]<-ttt$p.value}

res
res[order(res$p.value),]

plot(wild$maxstbla, cult$maxstbla, pch=16, cex=1.5)
text(wild$maxstbla, cult$maxstbla, wild$trs, col="grey",adj=-0.15)
abline(1,1)
