Codiversification_tests.Rmd

---
title: "Codiversification tests"
author: "Taichi Suzuki"
date: "10/1/2021"
output: html_document
---

##Introduction

This document describes the codes that are used to generate the codiversification statistics using three tests (Phytools, Paco, and Parafit) on two types of trees (majority-rule consensus phylogenies and best maximum likelihood phylogenies).

1.Phytools
2.Paco
3.Parafit


##Setup
```{r setup, include=FALSE}
#Set root r
knitr::opts_knit$set(root.dir = "/ebio/abt3_projects/Bifido_Coevolution/publication_codes/")
```


##Install packages
```{r}
library(ape)
library(phylotools)
library(phytools)
library(ggtree)
library(dynamicTreeCut)
library(stringr)
library(cluster)
```


##Data example
```{r}
#Name of the bacteria to test
Bac_name = "s__Collinsella_aerofaciens" 

#Best maximum likelihood tree out of 100 trees (output of StrainPhlAn)
Bac_tree_best = read.tree("/ebio/abt3_projects/Bifido_Coevolution/publication_codes/Example_data/s__Collinsella_aerofaciens_best.tree") 

#Best maximum likelihood tree out of 100 trees (output of SNPhylo)
Host_tree_best = read.tree("/ebio/abt3_projects/Bifido_Coevolution/publication_codes/Example_data/Host_tree_best.tree")

#Association file to match host tip and bacterial tip. 
assoc_f = read.table("/ebio/abt3_projects/Bifido_Coevolution/publication_codes/Example_data/HostTip_BacTip_assoc.txt", header = TRUE)
```
#Create majority-rule consensus trees
```{r}
#Majority-rule consensus bacterial tree by collapsing branches with low bootstrap support <50
Bac_tree_cons = as.polytomy(Bac_tree_best, feature = 'node.label', fun=function(x) as.numeric(x) < 0.5) 
#Majority-rule consensus host tree by collapsing braches with low bootstrap support <50
Host_tree_cons = as.polytomy(Host_tree_best, feature = 'node.label', fun=function(x) as.numeric(x) < 0.5)
```



##1. Phytools function
```{r}
#This function takes 6 inputs and calculate the mean and standard deviation of test statistics and p-values using N permutations
#1. Bac_name: name of bacteria
#2. Bac_tree: bacterial tree (best tree or consensus tree)
#3. Host_tree: host tree (best tree or consensus tree)
#4. N_perm: the number of permutations for a given cospeciation function
#5. assoc_f: association file to match the host tips and bacterial tips
#6. output_dir: output directory of the resulting file 

phytools.function = function(Bac_name, Bac_tree, Host_tree, N_perm, assoc_f, output_dir){
  
  #Get bac tips
  Bac_tip = Bac_tree$tip.label
  Bac_tip_df = as.data.frame(Bac_tip)
  
  #Filter association file
  Host_Bac_filtered = merge(x=Bac_tip_df, y=assoc_f)
  Host_Bac_filtered_m = as.matrix(Host_Bac_filtered)
  
  #Run cospeciation function with N_perm permutation 
  result = cospeciation(Bac_tree, Host_tree, distance = c("RF"), method = c("permutation"), assoc=Host_Bac_filtered_m, nsim = N_perm)

  #Change matrix to list, extract row, and then to vector
  result2 = unlist(result, recursive = FALSE, use.names= TRUE)
  obs_RFdist = as.numeric(result2$d)
  result3 = split(unlist(result2, use.names = FALSE), rep(names(result2), lengths(result2)))
  pvalues = as.numeric(result3$P.val)
  null = as.numeric(result3$d.null)
  
  #Write output
  sink(paste0(output_dir, Bac_name, "_n", N_perm, "_phytools_output.txt"), append=TRUE)
  print(Bac_name)
  print(result)
  sink()
}
```

#1.1. Phytools Example
```{r}
#output directory
phytools_example_output = "/ebio/abt3_projects/Bifido_Coevolution/publication_codes/Example_data/"
  
phytools.function(Bac_name, Bac_tree_best, Host_tree_best, 999, assoc_f, phytools_example_output)
```



##2. Paco function
```{r}
#This function takes 5 inputs and calculate the Paco test statistics and p-values using N permutations
#1. Bac_name: name of bacteria
#2. Bac_tree: bacterial tree (best tree or consensus tree)
#3. Host_tree: host tree (best tree or consensus tree)
#4. N_perm: the number of permutations for a given cospeciation function
#5. output_dir: output directory of the resulting file 

PACO.cutree.function = function(Bac_name, Bac_tree, Host_tree, N_perm, output_dir){

##Required subfunctions:
#PACo
  PACo <- function (H.dist, P.dist, HP.bin)
{ 
HP.bin <- which(HP.bin > 0, arr.in=TRUE)
H.PCo <- pcoa(H.dist, correction="cailliez")$vectors #Performs PCo of Host distances 
P.PCo <- pcoa(P.dist, correction="cailliez")$vectors #Performs PCo of Parasite distances
H.PCo <- H.PCo[HP.bin[,1],] #adjust Host PCo vectors 
P.PCo <- P.PCo[HP.bin[,2],]  ##adjust Parasite PCo vectors
list (H.PCo = H.PCo, P.PCo = P.PCo)
}

#RemoveDups to select unique bacterial strain to create cuttree 
RemoveDups <- function(df, column) {
  inds = sample(1:nrow(df))  
  df   = df[inds, ]

  dups = duplicated(df[, column])
  df   = df[!dups, ]
  inds = inds[!dups]

  df[sort(inds, index=T)$ix, ]
}


##Main function  
##1. Create Bac tree using cutree
#Get tip labels
Bac_tip_label = Bac_tree$tip.label

#Convert to distance matrix
Bac_dist = cophenetic(Bac_tree) 
  
#Convert to hclust
Bac_hclust = diana(Bac_dist)
Bac_hclust_1 = as.hclust(Bac_hclust)
 
#Cutree
Bac_cutree = cutreeDynamic(Bac_hclust_1, cutHeight = 0.99, method = "tree", minClusterSize = 2) 
Bac_cutree2 = cbind(Bac_tip_label,Bac_cutree)
Bac_cutree3 = as.data.frame(Bac_cutree2)
Bac_cutree3$Bac_cutree = sub("^", "X", Bac_cutree3$Bac_cutree)

#Filter cutree IDs
Bac_cutree3_unique = RemoveDups(Bac_cutree3, "Bac_cutree")
Bac_cutree3_unique_ID = Bac_cutree3_unique$Bac_tip_label

#Filter bac tree
BacTree_cutree = keep.tip(Bac_tree, Bac_cutree3_unique_ID)

#rename bac tree
BacTree_cutree_renamed = sub.taxa.label(BacTree_cutree, Bac_cutree3)
Bac_dist_cutree = cophenetic(BacTree_cutree_renamed) 


##2. Edit host tree
#Host tree with IDs that exist in Bac tree
Host_tree_filter = keep.tip(Host_tree, Bac_tip_label)
Host_tree_filter_m = cophenetic(Host_tree_filter)


##3. Create HP file
Bac_cutree3_filtered = dcast(Bac_cutree3, Bac_tip_label~Bac_cutree, length, value.var = "Bac_cutree")
Bac_cutree3_filtered2 = data.frame(Bac_cutree3_filtered[,-1], row.names = Bac_cutree3_filtered[,1])
HP_file = as.matrix(Bac_cutree3_filtered2)


#4. Rename variables to input to PACO function
Bac_dist = Bac_dist_cutree #bac tree that are representing one strain per cutree groups
host_dist = Host_tree_filter_m #host tree filtered by bac tree IDs
HP = HP_file


#5. PACO function
PACo.fit <- PACo(host_dist, Bac_dist, HP)
HP.proc <- procrustes(PACo.fit$H.PCo, PACo.fit$P.PCo) #Procrustes Ordination 
NLinks = sum(HP)

m2.obs <- HP.proc$ss #observed sum of squares
N.perm = N_perm #set number of permutations for testing
P.value = 0
seed <-.Random.seed[trunc(runif(1,1,626))]
set.seed(seed)
    #set.seed(5) ### use this option to obtain reproducible randomizations


for (n in c(1:N.perm))
	{ 
	if (NLinks <= nrow(HP) | NLinks <= ncol(HP)) 	#control statement to avoid all parasites beig associated to a single host 
		{	flag2 <- TRUE 
			while (flag2 == TRUE)	{ 
		HP.perm <- t(apply(HP,1,sample))
		if(any(colSums(HP.perm) == NLinks)) flag2 <- TRUE else flag2 <- FALSE
									}  
		} else { HP.perm <- t(apply(HP,1,sample))} #permutes each HP row independently
		PACo.perm <- PACo(host_dist, Bac_dist, HP.perm)
		m2.perm <- procrustes(PACo.perm$H.PCo, PACo.perm$P.PCo)$ss #randomized sum of squares
		#write (m2.perm, file = paste0(output_dir,Bac_name,"_n",N_perm,"_m2_null.txt"), sep ="\t", append =TRUE) #option to save m2 from each permutation
		#write (m2.obs, file = paste0(output_dir,Bac_name,"_n",N_perm,"_m2_obs.txt"), sep ="\t", append =FALSE) #option to save m2 observed
		if (m2.perm <= m2.obs)
		{P.value = P.value + 1} 
	}
P.value <- P.value/N.perm


#Calculate m2_null mean and sd
m2_null = read.table(file = paste0(output_dir,Bac_name,"_n",N_perm,"_m2_null.txt"), sep ="\t", header = FALSE)
mean_null = mean(m2_null$V1)
sd_null = sd(m2_null$V1)
m2_obs = read.table(file = paste0(output_dir,Bac_name,"_n",N_perm,"_m2_obs.txt"), sep ="\t", header = FALSE)
  
#print output
sink(paste0(output_dir, Bac_name, "_n", N_perm, "_PACO_output.txt"), append=TRUE)
  print(Bac_name)
  print(paste0("p-value = ",P.value))
  print(paste0("m2_obs = ",m2_obs))
  print(paste0("Mean m2_null = ",mean_null))
  print(paste0("SD m2_null = ",sd_null))
sink()
}

```

#2.1. PACo Example
```{r}
#output directory
paco_example_output = "/ebio/abt3_projects/Bifido_Coevolution/publication_codes/Example_data/"
  
PACO.cutree.function(Bac_name, Bac_tree_best, Host_tree_best, 999, paco_example_output)
```



##3. Parafit function
```{r}
#This function takes 5 inputs and calculate the mean and standard deviation of Parfit test statistics and p-values using N permutations
#1. Bac_name: name of bacteria
#2. Bac_tree: bacterial tree (best tree or consensus tree)
#3. Host_tree: host tree (best tree or consensus tree)
#4. N_perm: the number of permutations for a given cospeciation function
#5. output_dir: output directory of the resulting file 

parafit.cutree.function = function(Bac_name, Bac_tree, Host_tree, N_perm, output_dir){
  
##1. Create Bac tree using cutree
#Get tip labels
Bac_tip_label = Bac_tree$tip.label

#Convert to distance matrix
Bac_dist = cophenetic(Bac_tree) 
  
#Convert to hclust
Bac_hclust = diana(Bac_dist)
Bac_hclust_1 = as.hclust(Bac_hclust)
 
#Cutree
Bac_cutree = cutreeDynamic(Bac_hclust_1, cutHeight = 0.99, method = "tree", minClusterSize = 2) 
Bac_cutree2 = cbind(Bac_tip_label,Bac_cutree)
Bac_cutree3 = as.data.frame(Bac_cutree2)
Bac_cutree3$Bac_cutree = sub("^", "X", Bac_cutree3$Bac_cutree)


#RemoveDups to select unique bacterial strain to create cuttree 
RemoveDups <- function(df, column) {
  inds = sample(1:nrow(df))  
  df   = df[inds, ]

  dups = duplicated(df[, column])
  df   = df[!dups, ]
  inds = inds[!dups]

  df[sort(inds, index=T)$ix, ]
}

#Filter cutree IDs
Bac_cutree3_unique = RemoveDups(Bac_cutree3, "Bac_cutree")
Bac_cutree3_unique_ID = Bac_cutree3_unique$Bac_tip_label

#Filter bac tree
BacTree_cutree = keep.tip(Bac_tree, Bac_cutree3_unique_ID)


#rename bac tree
BacTree_cutree_renamed = sub.taxa.label(BacTree_cutree, Bac_cutree3)
Bac_dist_cutree = cophenetic(BacTree_cutree_renamed) 


##2. Edit host tree
#Host tree with IDs that exist in Bac tree
Host_tree_filter = keep.tip(Host_tree, Bac_tip_label)
Host_tree_filter_m = cophenetic(Host_tree_filter)

##3. Create HP file
Bac_cutree3_filtered = dcast(Bac_cutree3, Bac_tip_label~Bac_cutree, length, value.var = "Bac_cutree")
Bac_cutree3_filtered2 = data.frame(Bac_cutree3_filtered[,-1], row.names = Bac_cutree3_filtered[,1])
HP_file = as.matrix(Bac_cutree3_filtered2)


##Parafit function
##############Pierre Legendre, May 2009#########################
'parafit_Legendre' <-
	function(host.D, para.D, HP, nperm=999, test.links=FALSE, seed=NULL, correction="none", silent=FALSE)
# Test of host-parasite coevolution
# host.D = host distance or patristic matrix (class dist or matrix)
# para.D = parasite distance or patristic matrix (class dist or matrix)
# HP = host-parasite link matrix (n.host, n.para)

{
epsilon <- sqrt(.Machine$double.eps)
if(is.null(seed)) {
	runif(1)
	seed <- .Random.seed[trunc(runif(1,1,626))]
	}
HP <- as.matrix(HP)
 
host.D <- as.matrix(host.D)
host.pc <- pcoa(host.D, correction=correction)
if(host.pc$correction[2] == 1) {
	if(min(host.pc$values[,2]) < -epsilon) stop('Host D matrix has negative eigenvalues. Rerun with correction="lingoes" or correction="cailliez"')
	sum.host.values.sq <- sum(host.pc$values[,1]^2)
	host.vectors <- host.pc$vectors
	} else {
	sum.host.values.sq <- sum(host.pc$values[,2]^2)
	host.vectors <- host.pc$vectors.cor
	}
n.host <- nrow(host.D)
 
para.D <- as.matrix(para.D)
para.pc <- pcoa(para.D, correction=correction)
if(para.pc$correction[2] == 1) {
	if(min(para.pc$values[,2]) < -epsilon) stop('Parasite D matrix has negative eigenvalues. Rerun with correction="lingoes" or correction="cailliez"')
	sum.para.values.sq <- sum(para.pc$values[,1]^2)
	para.vectors <- para.pc$vectors
	} else {
	sum.para.values.sq <- sum(para.pc$values[,2]^2)
	para.vectors <- para.pc$vectors.cor
	}
n.para <- nrow(para.D)
 
if(!silent) cat("n.hosts =", n.host, ", n.parasites =", n.para,'\n')
 
a <- system.time({
tracemax <- max(sum.host.values.sq, sum.para.values.sq)
 
if(n.host == n.para) {
	if(!silent) cat("The function cannot check if matrix HP has been entered in the right way.",'\n')
	if(!silent) cat("It will assume that the rows of HP are the hosts.",'\n')
	} else {
	temp <- dim(HP)
	if(temp[1] == n.host) {
		if(temp[2] != n.para) stop("Matrices host.D, para.D and HP not comformable")
		} else if(temp[2] == n.host) {
			if(temp[1] != n.para) stop("Matrices host.D, para.D and HP not comformable")
			HP <- t(HP)
			if(!silent) cat("Matrix HP has been transposed for comformity with host.D and para.D.",'\n')
		} else {
		stop("Matrices host.D, para.D and HP not comformable")
		}
	}
p.per.h <- apply(HP, 1, sum)
h.per.p <- apply(HP, 2, sum)
#
# Compute and test the global statistics
mat.4 <- t(host.vectors) %*% HP %*% para.vectors
global <- sum(mat.4^2)
if(nperm > 0) {
	set.seed(seed)
	nGT <- 1
	global.perm <- NA
	for(i in 1:nperm) {
		HP.perm <- apply(HP, 2, sample)
		mat.4.perm <- t(host.vectors) %*% HP.perm %*% para.vectors
		global.perm <- c(global.perm, sum(mat.4.perm^2))
		if(global.perm[i+1] >= global) nGT <- nGT+1
		}
	global.perm <- global.perm[-1]
	p.global <- nGT/(nperm+1)
	} else { p.global <- NA }
 
#
# Test individual H-P links
if(test.links) {
	# 1. Create the list of H-P pairs
	list.hp <- which( t(cbind(HP,rep(0,n.host))) > 0)
	HP.list <- cbind((list.hp %/% (n.para+1))+1, list.hp %% (n.para+1))
	colnames(HP.list) <- c("Host","Parasite")
	n.links <- length(list.hp)
 
	stat1 <- NA
	stat2 <- NA
	p.stat1 <- NA
	p.stat2 <- NA
	for(k in 1:n.links) {
		#
		# 2. Compute reference values of link statistics
		HP.k <- HP
		HP.k[HP.list[k,1], HP.list[k,2]] <- 0
		mat.4.k <- t(host.vectors) %*% HP.k %*% para.vectors
		trace.k <- sum(mat.4.k^2)
		stat1 <- c(stat1, (global-trace.k))
		den <- tracemax-global
		if(den > epsilon) { 
			stat2 <- c(stat2, stat1[k+1]/den) 
			} else { 
			stat2 <- c(stat2, NA) 
			}
		#
		# 3. Test link statistics by permutations
		if(nperm > 0) {
			set.seed(seed)
			nGT1 <- 1
			nGT2 <- 1
			nperm2 <- nperm
			#
			for(i in 1:nperm) {
				HP.k.perm <- apply(HP.k, 2, sample)
				mat.4.k.perm <- t(host.vectors) %*% HP.k.perm %*% para.vectors
				trace.k.perm <- sum(mat.4.k.perm^2)
				stat1.perm <- global.perm[i]-trace.k.perm
				if(stat1.perm >= stat1[k+1]) nGT1 <- nGT1+1
				#
				if(!is.na(stat2[k+1])) {
					den <- tracemax-global.perm[i]
					if(den > epsilon) {
						stat2.perm <- stat1.perm/den 
						if(stat2.perm >= stat2[k+1]) nGT2 <- nGT2+1
						} else {
						nperm2 <- nperm2-1
						# if(!silent) cat("In permutation #",i,"den < epsilon",'\n')
						}
					}
				}
			p.stat1 <- c(p.stat1, nGT1/(nperm+1))
			if(!is.na(stat2[k+1])) {
				p.stat2 <- c(p.stat2, nGT2/(nperm2+1))
				} else {
				p.stat2 <- c(p.stat2, NA) ### Error in previous version, corrected here
				}
			} else { 
			p.stat1 <- c(p.stat1, NA)     ### Error in previous version, corrected here
			p.stat2 <- c(p.stat2, NA)     ### Error in previous version, corrected here
			}
		}
	#
	link.table <- cbind(HP.list, stat1[-1], p.stat1[-1], stat2[-1], p.stat2[-1])
	colnames(link.table) = c("Null","Host","Parasite","F1.stat","p.F1","F2.stat","p.F2")
	out <-list(Null=global.perm,ParaFitGlobal=global, p.global=p.global, link.table=link.table, para.per.host=p.per.h, host.per.para=h.per.p, nperm=nperm)
	} else {
	if(!silent) cat("Rerun the program with option 'test.links=TRUE' to test the individual H-P links",'\n')
	out <-list(Null=global.perm,ParaFitGlobal=global, p.global=p.global, para.per.host=p.per.h, host.per.para=h.per.p, nperm=nperm)
	}
#
})
a[3] <- sprintf("%2f",a[3])
if(!silent) cat("Computation time =",a[3]," sec",'\n')
#
class(out) <- "parafit"
out
}
##########################

#Run parafit with N_perm permutation
parafit_result = parafit_Legendre(Host_tree_filter_m, Bac_dist_cutree, HP_file, nperm = N_perm, correction="cailliez")

#Write output
  sink(paste0(output_dir, Bac_name, "_n", N_perm, "_parafit_output.txt"), append=TRUE)
  print(Bac_name)
  print(parafit_result)
  sink()
}
```


#3.1. Parafit Example
```{r}
#output directory
parafit_example_output = "/ebio/abt3_projects/Bifido_Coevolution/publication_codes/Example_data/"
  
parafit.cutree.function(Bac_name, Bac_tree_best, Host_tree_best, 999, parafit_example_output)
```