R for Phylogenetic

Contract the UPGMA tree

##############LOAD PACKAGES###########################
install.packages("ape")
install.packages("BiocManager")
BiocManager::install("msa")
library(msa)
library(ape)

#################READ THE FASTA FILE############################
dna=readDNAStringSet("sequence.fasta")

################ALIGN THE SEQUENCES#############################
clu_alin = msa(dna,method = "ClustalW")

#############CONVERT TO DNAbin FORMAT###########################
Clu_bin=as.DNAbin(Clu_alin)

#############MAKE THE DISTANCE METRIX##########################
clu_dis = dist.dna(Clu_bin)

#############CONSTACT THE UPGMA TREE###########################
clu_upgm = hclust(clu_dis, method = "average")

############ADDITIONAL-DISPLAY PROPERLY########################
clu_upgm_den=as.dendrogram(clu_upgm)

############PLOT THE TREE######################################
plot(clu_upgm_den,horiz=T)



 

Neighbor Joining tree

#############INSTALL PACKAGES###########
install.packages("phangorn")
library(phangorn)
install.packages("ape")
install.packages("BiocManager")
BiocManager::install("msa")
library(msa)
library(ape)

#########READ THE SEQUENCE DATA###########
dna=readDNAStringSet("sequences.fasta")

#########ALIGN THE DATA##############
Clu_alin=msa(dna,method="ClustalW")

#########CONVERT TO THE DNAbin format###########
Clu_bin=as.DNAbin(Clu_alin)

#########BIULD THE PAIRWISE DISTANCE##############
Clu_dist = dist.dna(Clu_bin, model = "K80")

#########BUILD THE NJ TREE##############
Clu_nj = nj(Clu_dist)

#########PLOT THE NJ TREE##############
plot(Clu_nj,horiz =T)

########CONVERT THE DISTANCE MATRIX TO CHARACTER ALIGNMENT########
Clu_data=as.phyDat(Clu_alin)

#######CAL THE PARSIMONY SCORE##############
parsimony(Clu_nj,Clu_data)

####### UPDATA NJ TREE WITH THE PASIMONY SCORE########
Clu_nj_opt = optim.parsimony(Clu_nj,Clu_data)

#######PLOT THE PARIMONY OPTIMIZED NJ TREES#########
plot(Clu_nj_opt,horiz=T)

Maximum Parsimony Tree & Maximum Likelihood Tree

#############INSTALL PACKAGES###########
install.packages("phangorn")
library(phangorn)
install.packages("ape")
install.packages("BiocManager")
BiocManager::install("msa")
library(msa)
library(ape)

###########READ THE FILE###########
fastaa = readAAStringSet("sequences.fasta", format = "fasta")

########ALIGN THE DATA########
mus_alin = msa(fastaa,method = "Muscle")

###########CONVERT BIN_FORMAT#################
mus_bin = as.AAbin(mus_alin)

##########DISTANCE METRIX##########
mus_dis = dist.ml(mus_bin, model = "mtArt")

########## BUILD NJ TREE###############
mus_tree = nj(mus_dis)

###########PLOT NJ  TREEE##############
plot(mus_tree)


bstrees <- boot.phylo(mus_tree, mus_bin, FUN = 
                        function(e) nj(dist.ml(e,model = "JTT")),
                      B = 1000, trees = T)$trees

######### MP TREE BIULDING############
phyd=as.phyDat(mus_alin)
ptree = pratchet(phyd,minit = 100,trace = 0)
plot(ptree)

########## DISCRIPTIVE STAT############
CI(ptree,phyd)
CI(mus_tree,phyd)
#RI(bstrees,phyd)

######the likelihood of a phylogenetic tree###################
mltree = pml(mus_tree,phyd, k=4)
plot(mltree)

####################### Optimize the ml tree#################
opt_pml = optim.pml(mltree, optNni = TRUE, optBf = TRUE, optQ = TRUE, optGamma = 
                      TRUE )

######################## outgroup rooting#################
rooted_tree = root(opt_pml$tree, "Populus trichocarpa" )
plot(rooted_tree)

##########################Bootstrapping ml trees###################
ts = bootstrap.pml(opt_pml,bs = 100, trees = T , multicore = T ) 
plotBS(root(opt_pml$tree, "Populus trichocarpa" ),ts,type = 
         "phylogram" )