install.packages('ape')
library('ape')
install.packages('TreeDist')
library('TreeDist')
library(slendr)
init_env()
eate populations
Chimp <- population("Chimp", N = 5000, time = 7000000)
Chimp
Africa <- population("Africa", N = 1500, time = 1000000, parent = Chimp)
Africa
Africa <- population("Africa", N = 1500, time = 6000000, parent = Chimp)
Neanderthal <- population("Neanderthal", N = 1000, time = 550000, parent = Africa)
Europe <- population("Europe", N = 6000, time = 60000, parent= Africa)
model <- compile_model(
list(Chimp, Africa, Neanderthal, Europe),
generation_time = 30
)
model
plot_model(model, sizes=TRUE, log=FALSE)
plot_model(model, sizes=FALSE, log=FALSE)
plot_model(model, sizes=TRUE, log=FALSE)
plot_model(model, sizes=TRUE, log=TRUE)
Africa[[1]][1]
Africa[[1]][2]
Africa[[2]][1]
ts <- msprime(model, sequence_length = 100e6, recombination_rate = 1e-8)
ts
ts_mutations <- ts_mutate(ts,mutation_rate=1e-8)
ts_mutations
ts_genotypes(ts)
ts_genotypes(ts_mutations)
ts_samples(ts)
ts_samples(ts)$pop
length(ts_samples(ts)$pop)
samples = ts_samples(ts)
samples
samples = split(samples, samples$pop)
samples
samples = lapply(samples,pull,"name")
samples = lapply(samples,"name")
samples <- ts_samples(ts) %>%
split(., .$pop) %>%
lapply(pull, "name")
samples = ts_samples(ts)
samples = split(samples, samples$pop)
samples$Africa
head(samples$Europe, 3)
Europe_samples=samples$Europe
Europe_samples
Chimp_samples=samples["Chimp"]
Chimp_samples
Africa_samples=samples$Africa
Neanderthal_samples=samples$Neanderthal
ts_diversity(ts, sample_sets = list(samples))
ts_diversity(ts, sample_sets = list(Chimp_samples))
ts_diversity(ts, sample_sets = Chimp_samples)
ts_diversity(ts_mutations, sample_sets = Chimp_samples)
ts_diversity(ts_mutations, sample_sets = Neanderthal_samples)
ts_diversity(ts_mutations, sample_sets = Neanderthal_samples)
samples
samples = ts_samples(ts_mutations)
samples = split(samples, samples$pop)
### Get all european samples
Europe_samples=samples$Europe
### Do the same for other populations
Chimp_samples=samples$Chimp
## or
Chimp_samples=samples["Chimp"]
Africa_samples=samples$Africa
Neanderthal_samples=samples$Neanderthal
ts_diversity(ts_mutations, sample_sets = Neanderthal_samples)
View(Neanderthal_samples)
View(Neanderthal)
View(Neanderthal_samples)
Neanderthal_samples
Neanderthal_samples$name
ts_diversity(ts_mutations, sample_sets = Neanderthal_samples$name)
GENES=c("AHSG","ALB","AMBN","AMELX","AMELY","AMTN","COL1A1","COL17A1","ENAM","MMP20","ODAM")
for (GENE in GENES){print(paste(GENE,'_lol'))}
GENES=c("AHSG","ALB","AMBN","AMELX","AMELY","AMTN","COL1A1","COL17A1","ENAM","MMP20","ODAM")
for (GENE in GENES){print(paste(GENE,'_lol',sep=""))}
GENES=c("AHSG","ALB","AMBN","AMELX","AMELY","AMTN","COL1A1","COL17A1","ENAM","MMP20","ODAM")
for (GENE in GENES){print(paste0(GENE,'_lol',sep=""))}
GENES=c("AHSG","ALB","AMBN","AMELX","AMELY","AMTN","COL1A1","COL17A1","ENAM","MMP20","ODAM")
for (GENE in GENES){print(paste0(GENE,'_lol',"JK",sep=""))}
setwd("C:/Users/rjt939/Desktop/PhD/ILS - Hominids/Enamel_Hominds_Tree_Distances")
library('TreeDist')
GENES=c("AHSG","ALB","AMBN","AMELX","AMELY","AMTN","COL1A1","COL1A2","COL17A1","ENAM","MMP20","ODAM")
GENE="AMELX"
INFERED_TREE_PATH=paste0(GENE,'/',GENE,"_aln_e.phy_phyml_tree.txt",sep="")
INFERED_TREE_PATH
OUTPUT_FILES=paste0(GENE,'/',GENE,"_Tree_Distances.txt",sep="")
Tree2=ape::read.tree(INFERED_TREE_PATH) #### Infered Tree
Tips=length(Tree2$tip.label)
counter=0
for (T in 1:Tips){
TIP_HERE = Tree2$tip.label[T]
TIP_HERE = unlist(strsplit(TIP_HERE, split=paste('_',GENE,sep=''), fixed=TRUE))[1]
Tree2$tip.label[T]=TIP_HERE
counter=counter+1
}
Tree2
### names of Genes
GENES=c("AHSG","ALB","AMBN","AMELX","AMELY","AMTN","COL1A1","COL1A2","COL17A1","ENAM","MMP20","ODAM")
for (GENE in GENES){
INITIAL_TREE_PATH="Hominid_Tree.txt"
INFERED_TREE_PATH=paste0(GENE,'/',GENE,"_aln_e.phy_phyml_tree.txt",sep="")
OUTPUT_FILES=paste0(GENE,'/',GENE,"_Tree_Distances.txt",sep="")
###################### Import Trees ####################
Tree1=ape::read.tree(INITIAL_TREE_PATH) #### Input Tree
Tree2=ape::read.tree(INFERED_TREE_PATH) #### Infered Tree
Tips=length(Tree2$tip.label)
counter=0
for (T in 1:Tips){
TIP_HERE = Tree2$tip.label[T]
TIP_HERE = unlist(strsplit(TIP_HERE, split=paste('_',GENE,sep=''), fixed=TRUE))[1]
Tree2$tip.label[T]=TIP_HERE
counter=counter+1
}
print('Trees loaded')
##################### Calculate Distances #################
##### Normal Tree Distances
distance_Generic <- TreeDistance(Tree1, Tree2)
distance_Nye <- NyeSimilarity(Tree1, Tree2, normalize = TRUE)
distance_JRF <- JaccardRobinsonFoulds(Tree1,Tree2)
distance_Bod_Giaro <- MatchingSplitDistance(Tree1,Tree2)
distance_KC <- KendallColijn(Tree1,Tree2)
##### For comparisons of tree with not - equal number of leaves
Clustering_Entropy1=ClusteringEntropy(Tree1)
Clustering_Entropy2=ClusteringEntropy(Tree2)
Mutual_Clustering_Info=MutualClusteringInfo(Tree1, Tree2)
#### Output Distances into file
METRICS=c(distance_Generic,distance_Nye,distance_JRF,distance_Bod_Giaro,distance_KC,Mutual_Clustering_Info)
file.create(OUTPUT_FILES)
write(paste(METRICS,collapse='\t'), file = OUTPUT_FILES,append=FALSE)
}
