# clean workspace

rm(list=ls())

# read in libraries

library(ape)
library(vegan)
library(adephylo)
library(bipartite)
library(picante)
library(ade4)
library(adiv)
library(phangorn)
library(cooccur)
library(car)
library(maps)
library(phytools)
library(fossil)

# clean the data

myPath <- '~/path/to/my/files/'

# read in host by ASV matrix
lg_db <- read.table(paste0(myPath, 'infiles/lg_tab.txt'))

# read in ASV tree
lg_tree <- read.tree(paste0(myPath, 'infiles/lg_asv.tre'))

# check that all ASVs in the matrix are represented in the tree

sum(!colnames(lg_db) %in% lg_tree$tip.label)

# read in host metadata
lg_hosts <- read.csv(paste0(myPath, 'infiles/Table_S1.csv'))

# select localities of interest
lg_locs <- unique(lg_hosts$locality)

# read in host tree
lg_host_tree <- read.tree(paste0(myPath, 'infiles/lg_host.tre'))

# read in database to match host species calls to host tree
tax_match_tonini <- read.csv(paste0(myPath, 'infiles/lat.grad.tax.match.csv'))

# match host species in metadata to host species in tree
lg_db_sp <- unlist(lapply(rownames(lg_db), function(x){x1 <- lg_hosts[lg_hosts$sample.name == x,]; return(paste0(x1[6], '_', x1[7]))}))

# remove hosts that do not have an entry on the host tree
lg_db_match <- lg_db[lg_db_sp %in% tax_match_tonini$my.species,]

# remove ASVs with fewer than 5 reads
lg_db_match[lg_db_match < 6] <- 0

# plot histogram of the number of reads per host
hist(rowSums(lg_db_match), xlab='read depth', main='Histogram of read depths prior to truncation')

# remove hosts with fewer than 25000 reads
lg_db_match2 <- lg_db_match[rowSums(lg_db_match) > 25000,]

# plot histogram of the number of ASVs per host
lg_db_pa <- lg_db_match2
lg_db_pa[lg_db_pa > 0] <- 1

# clean up empty columns
lg_db_match4 <- lg_db_match2[, colSums(lg_db_match2) > 0]

# rarefy data to 25,000 ASVs per host and reformat
lg_reads_rare <- apply(lg_db_match4, 1, function(x){x2 <- x[x>0]; x3 <- sample(rep(names(x2), x2), 25000);x4 <- table(x3);fullVec <- c();for(i in 1:ncol(lg_db_match4)){if(colnames(lg_db_match4)[i] %in% names(x4)){fullVec[i] <- x4[names(x4) == colnames(lg_db_match4)[i]]}else{fullVec[i] <- 0}};return(fullVec)})

# reformat
lg_reads_rare <- t(lg_reads_rare)
colnames(lg_reads_rare) <- colnames(lg_db_match4)

# clean up empty columns
lg_reads_rare <- lg_reads_rare[, colSums(lg_reads_rare) > 0]

# convert matrix into a per-locality list of host x ASV matrices

lg_rare_locs <- unlist(lapply(rownames(lg_reads_rare), function(x){lg_hosts$locality[lg_hosts$sample.name == x]}))

# manually order localities north to south
lg_locs2 <- c('Michigan', 'Georgia', 'Nicaragua highlands', 'Nicaragua dry forest',  'Nicaragua wet forest', 'Peru wet forest')

lg_list3 <- lapply(lg_locs2, function(x){tmpNames <- rownames(lg_reads_rare)[lg_rare_locs == x]; tMat <- lg_reads_rare[rownames(lg_reads_rare) %in% tmpNames,]; tMat2 <- tMat[,colSums(tMat) > 0]; return(tMat2)})

# get phylum ID of each ASV from the tree

# find tree node labels starting with 'p__', which indicates a phylum identifier
phyNodes <- cbind((1:length(lg_tree$node.label))[grepl('p__', lg_tree$node.label)], lg_tree$node.label[grepl('p__', lg_tree$node.label)])

# use the 'Descendants' function to get the tips that descend from each node
asvPhyNo <- lapply(phyNodes[,1], function(x){Descendants(lg_tree, as.numeric(x) + length(lg_tree$tip.label), type='tips')})

# use the numerical values for the tips to get the named tip label 
asvPhy <- lapply(asvPhyNo, function(x){lg_tree$tip.label[unlist(x)]})

# match phylum names to list of tipes
names(asvPhy) <- lg_tree$node.label[grepl('p__', lg_tree$node.label)]

# convert list into vector
asvPhyVec <- rep(names(asvPhy), unlist(lapply(asvPhy, length)))
names(asvPhyVec) <- unlist(asvPhy)

asvPhyVec <- substr(asvPhyVec, 5, nchar(asvPhyVec) - 1)

asvPhyVec <- unlist(lapply(strsplit(asvPhyVec, '-'), function(x){x[1]}))

# check for ASVs with no phylum ID from the tree and append to vector

allASV <- unlist(lapply(lg_list3, colnames))

noPhyNames <- allASV[!allASV %in% names(asvPhyVec)]

noPhy <- rep('noPhy', length=length(noPhyNames))
names(noPhy) <- noPhyNames

asvPhyVec <- c(asvPhyVec, noPhy)

# print out number of ASVs per phylum
table(asvPhyVec)[rev(order(table(asvPhyVec)))]


############################
# set up plotting features #
############################

# vector of colors for common phyla
phyCols <- c('#f98e09', '#bc3754', '#fcffa4',  '#57106e', '#000004')
names(phyCols) <- c("Proteobacteria", "Firmicutes", "Planctomycetes", "Actinobacteria", "noPhy")

# vector of colors for locations
locCols <- c('#fde725', '#7ad151', '#22a884', '#2a788e', '#414487', '#440154')
names(locCols) <- lg_locs2

# location labels for plotting
locLabels <- c('Michigan', 'Georgia', 'Nicaragua highlands', 'Nicaragua dry forest', 'Nicaragua wet forest', 'Peru wet forest')

# find the number of sampled hosts per location
host_sp <- unlist(lapply(lg_list3, function(x){x1 <- rownames(x); length(unique(lg_hosts$species[lg_hosts$sample.name %in% x1]))}))

# add the number of host species sampled to location labesl for plotting
locLabels2 <- paste0(locLabels, ' (', host_sp, ')')

# make a two-column matrix of every unique pair of locations
pairwiseRow2 <- list()

for(i in 2:6){
	
	pairwiseRow2[[i-1]] <- lg_locs2[i:6]
	
}

pairwiseNames <- cbind(c(rep('Michigan', 5), rep('Georgia', 4), rep('Nicaragua highlands', 3), rep('Nicaragua dry forest', 2), 'Nicaragua wet forest'), c(unlist(pairwiseRow2)))

# for later use, find names of ASVs that are in Proteobacteria and Firmicutes 
namesProt <- names(asvPhyVec)[asvPhyVec == "Proteobacteria"]

namesFirm <- names(asvPhyVec)[asvPhyVec == "Firmicutes"]

# read in pre-processed dataset (see diversity script)

lg_list7 <- list()

for(i in 1:6){
	
	lg_list7[[i]] <- read.table(paste0(myPath, 'infiles/', lg_locs2[i], '_hostxotu'))
	
}

#################

# create a function to process our host x otu matrix and host phylogenetic tree to find host specificity of a set of ASVs
# formats data for the ses.pd function from the package 'picante'

# the function takes a vector of all hosts at a site, a list of host in which an ASV occurs, and a phylogenetic tree of the hosts

# this is equivalent to prep_ses in the diversity script, but formatted for the output object from the bipartite package 'metaComputeModules' command

# host_tot is a list of all hosts in a community, host_list is the hosts in which an OTU (or module) occurs, can be multiple vectors of hosts, host_tree is the phylogenetic tree of hosts

prep_ses_2 <- function(host_tot, host_list, host_tree){
	
	# take a list of vectors of hosts that are associated with a module and process them into a presence/absence matrix relative to all available hosts for a location
	host_list_tmp <- list()
	
	# format hosts into presence-absence vectors relative to full community
	for(i in 1:length(host_list)){
		
		tmp_h <- host_list[[i]]
		
		host_list_tmp[[i]] <- unlist(lapply(host_tot, function(x){if(x %in% tmp_h){return(1)}else(return(0))}))
		
	}
	
	# make into a matrix
	hostMat <- do.call(cbind, host_list_tmp)
	
	rownames(hostMat) <- host_tot
	colnames(hostMat) <- paste0('mod_', 1:ncol(hostMat))
	
	hostTmp <- colnames(t(hostMat))
		
	# take host names and find taxa that match the host tree	
	spTmp <- unlist(lapply(hostTmp, function(x){x1 <- lg_hosts[lg_hosts$sample.name == x,]; return(paste0(x1[6], '_', x1[7]))}));					
	spTonini1 <- lapply(spTmp, function(x){tax_match_tonini[tax_match_tonini $my.species %in% x, 2]})					
	spTonini <- unlist(spTonini1)
	spTonini_tab <- table(spTonini)
	
	spTonini_order <- spTonini[order(spTonini)]
	
	hostTmp2 <- t(hostMat)[, order(spTonini)]
	
	names_hostTmp2 <- list()
	
	for(i in 1:length(spTonini_tab)){
		names_hostTmp2[[i]] <- paste0(names(spTonini_tab)[i], '_', 1: spTonini_tab[i])
	}

	names_hostTmp2_vec <- unlist(names_hostTmp2)

	colnames(hostTmp2) <- names_hostTmp2_vec
	
	# take the full host tree (all locations) and clip down only to 	the hosts present in the full sample
	host_tree_clip <- keep.tip(host_tree, names(spTonini_tab))
			
	for(i in 1:ncol(hostTmp2)){
		
		host_tree_clip <- bind.tip(host_tree_clip, names_hostTmp2_vec[i], edge.length=0, where=which(host_tree_clip$tip.label == spTonini_order[i]))		
	}
	
	host_tree_clip <- drop.tip(host_tree_clip, names(spTonini_tab))
	
	# run ses.pd and return observed phylogenetic specificity
	ses <- ses.pd(hostTmp2, host_tree_clip)
	
	# get p value of clusters (ie, is the group significantly phylogenetically clustered?)
	pd_z <- ses$pd.obs.z
	
	# NaNs are returned when an ASV is in all hosts (sd is zero and dividing by zero gives NaN). Replace with max value of any other ASV in the dataset
	pd_z[pd_z == 'NaN'] <- max(pd_z[!is.na(pd_z)])
	
	names(pd_z) <- rownames(ses)
	
	return(pd_z)	
}

####################################

# resample the dataset and re-run cooccurrence analysis to account for sampling effects

# subsampled dataset
lg_list_cooc <- list()

# pairs with greater than expected overlap
lg_cooc_pos <- list()

# pairs with lower than expected overlap
lg_cooc_neg <- list()

# randomized host x ASV matrix
lg_list_rand <- list()

# pairs with greater than expected overlap for randomized matrix
lg_rand_pos <- list()

# pairs with lower than expected overlap for randomized matrix
lg_rand_neg <- list()


for(i in 1:50){
	
	# sample 15 hosts per location from full list
	lg_t1 <- lapply(lg_list3, function(x){x[sample(1:nrow(x), 15),]})
	
	# remove ASVs with no reads from each lcoation
	lg_t2 <- list()

	for(j in 1:6){
		
		t1 <- lg_t1[[j]]
		
		lg_t2[[j]] <- t1[, apply(t1, 2, function(x){sum(x == 0)}) < 14]
		
	}

	# subsample ASVs down to 120
	lg_t3 <- lapply(lg_t2, function(x){x[,sample(1:ncol(x), 120)]})
	
	lg_list_cooc[[i]] <- lg_t3
	
	# run cooccurrence analysis
	cooc <- lapply(lg_t3, function(x){x[x > 0] <- 1; cooccur(t(x))})

	# get positively and negatively cooccurring ASVs
	lg_cooc_pos[[i]] <- lapply(cooc, function(x){x$results[x$results$p_gt < 0.01,]})

	lg_cooc_neg[[i]] <- lapply(cooc, function(x){x$results[x$results$p_lt < 0.01,]})
	
	########################
	
	# randomize column order to create null dataset
	lg_null <- lapply(lg_t3, function(x){apply(x, 2, function(y){y[sample(1:length(y), length(y))]})})
	
	lg_list_rand[[i]] <- lg_null
	
	# rerun cooccurrence analysis on null dataset
	coocNull <- lapply(lg_null, function(x){x[x > 0] <- 1; cooccur(t(x))})

	lg_rand_pos[[i]] <- lapply(coocNull, function(x){x$results[x$results$p_gt < 0.01,]})

	lg_rand_neg[[i]] <- lapply(coocNull, function(x){x$results[x$results$p_lt < 0.01,]})
	
	# track iteration
	print(i)
}

length(lg_rand_neg)

# sort into locations first, then runs within locations
# report a single metric per run
# for matched depth comparisons, single null pair per real pair

# number of pairs (positively associated)
n_pair_pos <- list()
n_pair_null_pos <- list()

# mean prevalence of ASVs in positive pairs
mean_size_pos <- list()
mean_size_null_pos <- list()

# number of pairs (negatively associated)
n_pair_neg <- list()
n_pair_null_neg <- list()

#mean prevalence of ASVs in negative pairs
mean_size_neg <- list()
mean_size_null_neg <- list()

# patristic distance between ASVs in positive pairs (real dataset, null dataset, and random pairs of ASVs matched by prevalence)
pair_phylo <- list()
pair_phylo_null <- list()
pair_phylo_sub <- list()

# host clustering value ASVs in positive pairs (real dataset, null dataset, and random pairs of ASVs matched by prevalence)
pair_host <- list()
pair_host_null <- list()
pair_host_sub <- list()

# the structure of these loops reverses the order of the co-occurrence analysis result list from run -> location to location -> run

for(loc in 1:6){
	
	n_pair_pos_1 <- c()
	n_pair_null_pos_1 <- c()	
	n_pair_neg_1 <- c()
	n_pair_null_neg_1 <- c()
	
	mean_size_pos_1 <- c()
	mean_size_null_pos_1 <- c()
	mean_size_neg_1 <- c()
	mean_size_null_neg_1 <- c()
	
	pair_phylo_1 <- c()
	pair_phylo_null_1 <- c()
	pair_phylo_sub_1 <- c()
		
	pair_host_1 <- c()
	pair_host_null_1 <- c()
	pair_host_sub_1 <- c()

	
	for(run in 1:length(lg_cooc_pos)){
		
		# reorder results into loction, then run
		real_pos <- lg_cooc_pos[[run]][[loc]]
		null_pos <- lg_rand_pos[[run]][[loc]]
		real_neg <- lg_cooc_neg[[run]][[loc]]
		null_neg <- lg_rand_neg[[run]][[loc]]
		
		# get the number of positively and negativley cooccurring pairs in the real and null datasets
		n_pair_pos_1[run] <- nrow(real_pos)
		n_pair_null_pos_1[run] <- nrow(null_pos)
		n_pair_neg_1[run] <- nrow(real_neg)
		n_pair_null_neg_1[run] <- nrow(null_neg)
		
		# get the mean prevalence of the pairs
		mean_size_pos_1[run] <- mean(rowMeans(real_pos[,3:4]))
		mean_size_null_pos_1[run] <- mean(rowMeans(null_pos[,3:4]))
		mean_size_neg_1[run] <- mean(rowMeans(real_neg[,3:4]))
		mean_size_null_neg_1[run] <- mean(rowMeans(null_neg[,3:4]))
		
		######################################
		
		# get the patristic distance between pairs of interacting ASVs 
		asv_names <- colnames(lg_list_cooc[[run]][[loc]])
		
		# account for empty objects or results with only one pair
		if(nrow(real_pos) > 0){pair_phylo_1[run] <-mean(apply(real_pos, 1, function(x){sum(keep.tip(lg_tree, asv_names[c(unlist(x[1]), unlist(x[2]))])$edge.length)}))}else{pair_phylo_1[run] <- NaN}
		
		if(nrow(null_pos) > 0){pair_phylo_null_1[run] <- mean(apply(null_pos, 1, function(x){sum(keep.tip(lg_tree, asv_names[c(unlist(x[1]), unlist(x[2]))])$edge.length)}))}else{pair_phylo_null_1[run] <- NaN}
		
		asv_depths <- apply(lg_list_cooc[[run]][[loc]], 2, function(x){sum(x != 0)})
		
		# get patristic distance between randomly drawn ASVs with matching prevalence
		if(nrow(real_pos) > 0){pair_phylo_sub_1[run] <- mean(apply(real_pos, 1, function(x){x1 <- unlist(x[3]); x2 <- unlist(x[4]); x1.1 <- names(asv_depths)[asv_depths == x1]; x2.1 <- names(asv_depths)[asv_depths == x2]; x1.2 <- x1.1[sample(1:length(x1.1), 1)]; x2.2 <- x2.1[sample(1:length(x2.1), 1)]; return(sum(keep.tip(lg_tree, c(x1.2, x2.2))$edge.length))}))}else{pair_phylo_sub_1[run] <- NaN}
		
		######################################
		
		# get the phylogenetic similarity value (using ses.pd from picante) of the hosts in which both ASVs of a pair cooccur
		
		real_mat <- lg_list_cooc[[run]][[loc]]		
		null_mat <- lg_list_rand[[run]][[loc]]
		
		# account for empty and singleton pairs
		if(nrow(real_pos) > 0){	
			pair_host_tmp <- c()		
			for(j in 1:nrow(real_pos)){		
				tmpHostNames <- rownames(real_mat)[apply(real_mat[,unlist(real_pos[1,1:2])], 1, function(x){x[1] & x[2] > 0})]		
				pair_host_tmp[j] <- mean(prep_ses_2(rownames(real_mat), tmpHostNames, lg_host_tree))
				
			}			
			pair_host_1[run] <- mean(pair_host_tmp)
		}else{pair_host_1[run] <- NaN}
		
		############	
		
		# phylogenetic value of null overlap hosts
		if(nrow(null_pos) == 0){pair_host_null_1[run] <- NaN}else{
			
			pair_host_tmp_null <- c()
			
			for(j in 1:nrow(null_pos)){
		
				tmpHostNames <- rownames(real_mat)[apply(null_mat[,unlist(null_pos[1,1:2])], 1, function(x){x[1] & x[2] > 0})]
		
				pair_host_tmp_null[j] <- mean(prep_ses_2(rownames(real_mat), tmpHostNames, lg_host_tree))
			
			}			
			
			pair_host_null_1[run] <- mean(pair_host_tmp_null)
		}
 	
 		#############
 		
 		# overlap hosts from abundance matched ASVs 
 		
 		matching_host <- apply(real_pos, 1, function(x){x1 <- unlist(x[3]); x2 <- unlist(x[4]); x1.1 <- names(asv_depths)[asv_depths == x1]; x2.1 <- names(asv_depths)[asv_depths == x2]; x1.2 <- x1.1[sample(1:length(x1.1), 1)]; x2.2 <- x2.1[sample(1:length(x2.1), 1)]; return(c(x1.2,x2.2))})
 		
 		if(nrow(real_pos) > 0){
	 		pair_host_tmp_sub <- c()		
	 		for(j in 1:ncol(matching_host)){
	 			
	 			host_overlap <- cbind(real_mat[,colnames(real_mat) %in% matching_host[1,j]], real_mat[,colnames(real_mat) %in% matching_host[2,j]])
	 			host_overlap[host_overlap > 0] <- 1
	 			
	 			ho2 <- rowSums(host_overlap)
	 			
	 			tmpHostNames <- names(ho2)[ho2 > 1]
	 			
	 			if(length(tmpHostNames) > 1){
	 				pair_host_tmp_sub[j] <- mean(prep_ses_2(rownames(real_mat), tmpHostNames, lg_host_tree))
	 			}else{pair_host_tmp_sub[j] <- 0}		
	 		}
	 		 		
	 		pair_host_sub_1[run] <- mean(pair_host_tmp_sub)	
 		}else{pair_host_sub_1[run] <- NaN}
 		
 		
 		##########
 		
 		}
	
	mean_size_pos_1 <- mean_size_pos_1[mean_size_pos_1 != 'NaN']
	mean_size_null_pos_1 <- mean_size_null_pos_1[mean_size_null_pos_1 != 'NaN']
	mean_size_neg_1 <- mean_size_neg_1[mean_size_neg_1 != 'NaN']
	mean_size_null_neg_1 <- mean_size_null_neg_1[mean_size_null_neg_1 != 'NaN']
	
	pair_phylo_1 <- pair_phylo_1[pair_phylo_1 != 'NaN']
	pair_phylo_null_1 <- pair_phylo_null_1[pair_phylo_null_1 != 'NaN']
	pair_phylo_sub_1 <- pair_phylo_sub_1[pair_phylo_sub_1 != 'NaN']
	
	pair_host_null_1 <- as.numeric(pair_host_null_1[pair_host_null_1 != 'NaN'])
	pair_host_1 <- as.numeric(pair_host_1[pair_host_1 != 'NaN'])
	pair_host_sub_1 <- as.numeric(pair_host_sub_1[pair_host_sub_1 != 'NaN'])
	
	n_pair_pos[[loc]] <- n_pair_pos_1
	n_pair_null_pos[[loc]] <- n_pair_null_pos_1
	n_pair_neg[[loc]] <- n_pair_neg_1
	n_pair_null_neg[[loc]] <- n_pair_null_neg_1
	
	mean_size_pos[[loc]] <- mean_size_pos_1
	mean_size_null_pos[[loc]] <- mean_size_null_pos_1 
	mean_size_neg[[loc]] <- mean_size_neg_1 
	mean_size_null_neg[[loc]] <- mean_size_null_neg_1

	pair_phylo[[loc]] <- pair_phylo_1
	pair_phylo_null[[loc]] <- pair_phylo_null_1
	pair_phylo_sub[[loc]] <- pair_phylo_sub_1

	pair_host[[loc]] <- pair_host_1
	pair_host_null[[loc]] <- pair_host_null_1
	pair_host_sub[[loc]] <- pair_host_sub_1

}

 
# test outputs for statistical significance and process outputs for export 


# number of positively coocurring pairs

# format matrix for input into 'aov' function (base R anova command)
locLength <- list()
for(i in 1:length(lg_locs2)){locLength[[i]] <- rep(lg_locs2[i], length(n_pair_pos[[i]]))}

npairMat <- data.frame(cbind(unlist(n_pair_pos), unlist(locLength)))
colnames(npairMat) <- c('npair', 'loc')

# write out mean number of pairs per location
lapply(n_pair_pos, mean)

# perform anova, Tukey test between locations
npair_aov <- aov(npair ~ loc, npairMat)
summary(npair_aov)
TukeyHSD(npair_aov)
npair_tukey <- TukeyHSD(npair_aov)$loc

# format into table to write out
tabS3_1 <- apply(pairwiseNames, 1, function(x){npair_tukey[grepl(x[1], rownames(npair_tukey)) + grepl(x[2], rownames(npair_tukey)) == 2,]})
colnames(tabS3_1) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_1) <- paste0('number positive pairs ', rownames(tabS3_1))


# negatively cooccurring pairs 
locLength <- list()
for(i in 1:length(lg_locs2)){locLength[[i]] <- rep(lg_locs2[i], length(n_pair_neg[[i]]))}


npairMatNeg <- data.frame(cbind(unlist(n_pair_neg), unlist(locLength)))
colnames(npairMatNeg) <- c('npair', 'loc')

lapply(n_pair_neg, mean)

npairNeg_aov <- aov(npair ~ loc, npairMatNeg)
summary(npairNeg_aov)
TukeyHSD(npairNeg_aov)

npair_neg_tukey <- TukeyHSD(npair_aov)$loc


tabS3_2 <- apply(pairwiseNames, 1, function(x){npair_neg_tukey[grepl(x[1], rownames(npair_neg_tukey)) + grepl(x[2], rownames(npair_neg_tukey)) == 2,]})
colnames(tabS3_2) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_2) <- paste0('number negative pairs ', rownames(tabS3_2))


# use t-tests to compare numbers of pairs in each location to the matched null results

t.pos <- list()

for(loc in 1:6){t.pos[[loc]] <- t.test(n_pair_pos[[loc]], n_pair_null_pos[[loc]])}

# format to write out for supplemental table
tabS2_1 <- unlist(lapply(t.pos, function(x){x$p.value}))

t.neg <- list()

for(loc in 1:6){t.neg[[loc]] <- t.test(n_pair_neg[[loc]], n_pair_null_neg[[loc]])}

tabS2_2 <- unlist(lapply(t.neg, function(x){x$p.value}))

t.neg2pos <- list()

for(loc in 1:6){t.neg2pos[[loc]] <- t.test(n_pair_neg[[loc]], n_pair_pos[[loc]])}

tabS2_3 <- unlist(lapply(t.neg2pos, function(x){x$p.value}))

# Repeat for prevalence of ASVs in cooccurring pairs

locLength <- list()
for(i in 1:length(lg_locs2)){locLength[[i]] <- rep(lg_locs2[i], length(mean_size_pos[[i]]))}


sizeMat <- data.frame(cbind(unlist(mean_size_pos), unlist(locLength)))
colnames(sizeMat) <- c('size', 'loc')

size_aov <- aov(size ~ loc, sizeMat)
summary(size_aov)
TukeyHSD(size_aov)

pair_size_tukey <- TukeyHSD(size_aov)$loc

tabS3_3 <- apply(pairwiseNames, 1, function(x){pair_size_tukey[grepl(x[1], rownames(pair_size_tukey)) + grepl(x[2], rownames(pair_size_tukey)) == 2,]})
colnames(tabS3_3) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_3) <- paste0('mean abundance positive pairs ', rownames(tabS3_3))


t.size <- list()

for(loc in 1:6){t.size[[loc]] <- t.test(mean_size_pos[[loc]], mean_size_null_pos[[loc]])}

tabS2_4 <- unlist(lapply(t.size, function(x){x$p.value}))

# repeat for phylogenetic clustering of shared hosts

locLength <- list()
for(i in 1:length(lg_locs2)){locLength[[i]] <- rep(lg_locs2[i], length(pair_host[[i]]))}

hostPhyMat <- data.frame(cbind(unlist(pair_host), unlist(locLength)))
colnames(hostPhyMat) <- c('host', 'loc')

host_aov <- aov(host ~ loc, hostPhyMat)
summary(host_aov)
TukeyHSD(host_aov)

pair_host_tukey <- TukeyHSD(host_aov)$loc

tabS3_5 <- apply(pairwiseNames, 1, function(x){pair_host_tukey[grepl(x[1], rownames(pair_host_tukey)) + grepl(x[2], rownames(pair_host_tukey)) == 2,]})
colnames(tabS3_5) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_5) <- paste0('specialization shared hosts ', rownames(tabS3_5))

t.host <- list()

for(loc in 1:6){t.host[[loc]] <- t.test(pair_host[[loc]], pair_host_null[[loc]])}

tabS2_7 <- unlist(lapply(t.host, function(x){x$p.value}))


t.host.match <- list()

for(loc in 1:6){t.host.match[[loc]] <- t.test(pair_host[[loc]], pair_host_sub[[loc]])}


tabS2_8 <- unlist(lapply(t.host.match, function(x){x$p.value}))


# repeat for patristic distance between 

locLength <- list()
for(i in 1:length(lg_locs2)){locLength[[i]] <- rep(lg_locs2[i], length(pair_phylo[[i]]))}


asvPhyMat <- data.frame(cbind(unlist(pair_phylo), unlist(locLength)))
colnames(asvPhyMat) <- c('asv', 'loc')

asv_aov <- aov(asv ~ loc, asvPhyMat)
summary(asv_aov)
TukeyHSD(asv_aov)

pair_phylo_tukey <- TukeyHSD(asv_aov)$loc

tabS3_4 <- apply(pairwiseNames, 1, function(x){pair_phylo_tukey[grepl(x[1], rownames(pair_phylo_tukey)) + grepl(x[2], rownames(pair_phylo_tukey)) == 2,]})
colnames(tabS3_4) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_4) <- paste0('mean phylo distance positive pairs ', rownames(tabS3_4))

t.phy <- list()

for(loc in 1:6){t.phy[[loc]] <- t.test(pair_phylo[[loc]], pair_phylo_null[[loc]])}

tabS2_5 <- unlist(lapply(t.phy, function(x){x$p.value}))

t.phy.match <- list()

for(loc in 1:6){t.phy.match[[loc]] <- t.test(pair_phylo[[loc]], pair_phylo_sub[[loc]])}

tabS2_6 <- unlist(lapply(t.phy.match, function(x){x$p.value}))


###########################################

# get patristic distance between nearest neighbors and mean patristic distance within each location (standard dataset)

tipDist <- c()

meanDist <- c()

for(i in 1:6){
	
	# get local tree
	lg_sub <- keep.tip(lg_tree, colnames(lg_list7[[i]]))
	
	# get mean patristic distance
	meanDist[i] <- mean(as.dist(cophenetic(lg_sub)))
	
	# get all subtree tips, keep subtrees with two tips (nearest neighbors)
	lg_subtree <- listTips(lg_sub)
	
	lg_subtree2 <- lg_subtree[unlist(lapply(lg_subtree, length)) == 2]
	
	lg_subtree3 <- lapply(lg_subtree2, function(x){keep.tip(lg_tree, x)})
	
	# get distance
	tipDist[i] <- mean(unlist(lapply(lg_subtree3, distTips)))
}

#########################################
# code for figure 5 #####################
#########################################


quartz(width=6, height=7)
par(mfrow=c(2, 1))
par(xpd=NA)
par(mar=c(7, 4, 4, 1))

# retain bonferroni corrected signficant comparisons for number of ASVs 

npair_tukey_sig <- npair_tukey[npair_tukey[,4] < 0.0001,]

sig_comparisons <- list()

for(i in 1:length(locCols)){
	
	t1 <- rownames(npair_tukey_sig)[grepl(names(locCols)[i], rownames(npair_tukey_sig))]
	
	t2 <- unlist(lapply(strsplit(t1, '-'), function(x){x[x != names(locCols)[i]]}))
	
	t3 <- letters[1:6][names(locCols) %in% t2]
	
	sig_comparisons[[i]] <- paste0(t3, collapse='')
}

# boxplots of number of pairs for positive and negative interactions

plot(c(0, 32), c(0, 60), type='n', axes=F, xlab='', ylab='number of pairs')
axis(2, seq(0, 60, by=20), seq(0, 60, by=20), las=2, cex.axis=1)
text(seq(3.5, 28.5, by=5), rep(-5, 6), labels= locLabels, srt=90, pos=2, cex=0.6)
text(-5, 70, 'A)')

boxplot(n_pair_pos, add=T, at=seq(1, 26, by=5), col= locCols, lwd=1.5, border = 'black', axes=F)
boxplot(n_pair_null_pos, add=T, at=seq(2, 27, by=5), col= locCols, lwd=1.5, axes=F, border='grey40')

boxplot(n_pair_neg, add=T, at=seq(3, 28, by=5), col= adjustcolor(locCols, alpha.f=0.5), lwd=1.5, border = 'black', axes=F)
boxplot(n_pair_null_neg, add=T, at=seq(4, 29, by=5), col= adjustcolor(locCols, alpha.f=0.5), lwd=1.5,  border = 'grey40', axes=F)

# text(seq(1, 26, by=5), rep(70, 6), unlist(sig_comparisons))

################

# comparison and boxplots for mean prevalence of ASVs in pairs

size_tukey_sig <- pair_size_tukey[pair_size_tukey[,4] < 0.0001,]

sig_comparisons_size <- list()

for(i in 1:length(locCols)){
	
	t1 <- rownames(size_tukey_sig)[grepl(names(locCols)[i], rownames(size_tukey_sig))]
	
	t2 <- unlist(lapply(strsplit(t1, '-'), function(x){x[x != names(locCols)[i]]}))
	
	t3 <- letters[1:6][names(locCols) %in% t2]
	
	sig_comparisons_size[[i]] <- paste0(t3, collapse='')
}


plot(c(0, 32), c(3, 12), type='n', axes=F, xlab='', ylab='mean prevalence')
axis(2, seq(3, 12, by=3), seq(3, 12, by=3), las=2, cex.axis=1)
text(seq(3.5, 28.5, by=5), rep(2, 6), labels= locLabels, srt=90, pos=2, cex=0.6)
text(-5, 15, 'B)')

boxplot(mean_size_pos, add=T, at=seq(1, 26, by=5), col= locCols, lwd=1.5, border = 'black', axes=F)
boxplot(mean_size_null_pos, add=T, at=seq(2, 27, by=5), col= locCols, lwd=1.5, axes=F, border='grey40')

boxplot(mean_size_neg, add=T, at=seq(3, 28, by=5), col= adjustcolor(locCols, alpha.f=0.5), lwd=1.5, border = 'black', axes=F)
boxplot(mean_size_null_neg, add=T, at=seq(4, 29, by=5), col= adjustcolor(locCols, alpha.f=0.5), lwd=1.5,  border = 'grey40', axes=F)

legend(3, 20, legend=c('real data', 'randomized'), fill=c('black', 'grey60'), bty='n', ncol=2)

# text(seq(1, 26, by=5), rep(12.5, 6), unlist(sig_comparisons_size))

#########################################
# code for figure 6 #####################
#########################################

# boxplots and significance comparisons for ASV and host phylogeny

quartz(width=7, height=7)
par(mfrow=c(2, 1))
par(xpd=NA)
par(mar=c(8, 4, 2, 1))


phylo_tukey_sig <- pair_phylo_tukey[pair_phylo_tukey[,4] < 0.0001,]

sig_comparisons_phylo <- list()

for(i in 1:length(locCols)){
	
	t1 <- rownames(phylo_tukey_sig)[grepl(names(locCols)[i], rownames(phylo_tukey_sig))]	
	t2 <- unlist(lapply(strsplit(t1, '-'), function(x){x[x != names(locCols)[i]]}))	
	t3 <- letters[1:6][names(locCols) %in% t2]	
	sig_comparisons_phylo[[i]] <- paste0(t3, collapse='')
}

plot(c(0, 24), c(0, 1.5), type='n', axes=F, xlab='', ylab='bacteria phylo. dist')
axis(2, seq(0, 1.5, by=0.5), seq(0, 1.5, by=0.5), las=2)
text(seq(3, 26, by=4), rep(-0.1, 6), labels= locLabels, srt=90, pos=2, cex=0.7)
text(-3, 1.7, 'A)')

boxplot(pair_phylo, add=T, at=seq(1, 22, by=4), col= locCols, lwd=1.5, border = 'black', axes=F)
boxplot(pair_phylo_sub, add=T, at=seq(2, 23, by=4), col= locCols, lwd=1.5, border='grey40',   axes=F)
boxplot(pair_phylo_null, add=T, at=seq(3, 24, by=4), col= locCols, lwd=1.5, border = 'grey60', axes=F)

lines(c(1, 3), c(tipDist[[1]], tipDist[[1]]), lwd=2)
lines(c(5, 8), c(tipDist[[2]], tipDist[[2]]), lwd=2)
lines(c(9, 12), c(tipDist[[3]], tipDist[[3]]), lwd=2)
lines(c(13, 16), c(tipDist[[4]], tipDist[[4]]), lwd=2)
lines(c(17, 20), c(tipDist[[5]], tipDist[[5]]), lwd=2)
lines(c(21, 24), c(tipDist[[6]], tipDist[[6]]), lwd=2)

lines(c(1, 3), c(meanDist[[1]], meanDist[[1]]), lwd=2)
lines(c(5, 8), c(meanDist[[2]], meanDist[[2]]), lwd=2)
lines(c(9, 12), c(meanDist[[3]], meanDist[[3]]), lwd=2)
lines(c(13, 16), c(meanDist[[4]], meanDist[[4]]), lwd=2)
lines(c(17, 20), c(meanDist[[5]], meanDist[[5]]), lwd=2)
lines(c(21, 24), c(meanDist[[6]], meanDist[[6]]), lwd=2)

# text(seq(1, 22, by=4), rep(1.6, 6), unlist(sig_comparisons_phylo))

#################

plot(c(0, 24), c(-2, 2), type='n', axes=F, xlab='', ylab='host phylo. diversity')
axis(2, seq(-2, 2, by=1), seq(-2, 2, by=1), las=2)
text(seq(3, 26, by=4), rep(-2, 6), labels= locLabels, srt=90, pos=2, cex=0.7)
text(-3, 2.5, 'B)')

host_tukey_sig <- pair_host_tukey[pair_host_tukey[,4] < 0.0001,]

sig_comparisons_host <- list()

for(i in 1:length(locCols)){
	
	t1 <- rownames(host_tukey_sig)[grepl(names(locCols)[i], rownames(host_tukey_sig))]	
	t2 <- unlist(lapply(strsplit(t1, '-'), function(x){x[x != names(locCols)[i]]}))	
	t3 <- letters[1:6][names(locCols) %in% t2]	
	sig_comparisons_host[[i]] <- paste0(t3, collapse='')
}


boxplot(pair_host, add=T, at=seq(1, 22, by=4), col= locCols, lwd=1.5, border = 'black', axes=F)
boxplot(pair_host_sub, add=T, at=seq(2, 23, by=4), col= locCols, lwd=1.5, border='grey40', axes=F)
boxplot(pair_host_null, add=T, at=seq(3, 24, by=4), col= locCols, lwd=1.5, border = 'grey60', axes=F)

legend(-3, 4, legend=c('real data', 'abundance matched ASVs', 'randomized'), fill=c('black', 'grey40', 'grey60'), bty='n', ncol=3)

# text(seq(1, 22, by=4), rep(2.2, 6), unlist(sig_comparisons_host))
# text(seq(1, 22, by=4), rep(2.2, 6), c('', '', '', '', 'f', 'e'))

###########################################################################################
# get phyla of interacting pairs, check for excess of Proteobacteria-Proteobacteria pairs #
###########################################################################################

tax_pos <- list()
tax_null_pos <- list()
tax_phylo_sub <- list()

for(loc in 1:6){
	
	tax_pos_1 <- list()
	tax_null_pos_1 <- list()
	tax_phylo_sub_1 <- list()
	
	for(run in 1:50){
		
		real_pos <- lg_cooc_pos[[run]][[loc]]
		null_pos <- lg_rand_pos[[run]][[loc]]
		
		######################################
		
		asv_names <- colnames(lg_list_cooc[[run]][[loc]])
		asv_depths <- apply(lg_list_cooc[[run]][[loc]], 2, function(x){sum(x != 0)})
		
		tax_pos_tmp <- cbind(unlist(lapply(asv_names[real_pos$sp1], function(x){asvPhyVec[names(asvPhyVec) == x][1]})), unlist(lapply(asv_names[real_pos$sp2], function(x){asvPhyVec[names(asvPhyVec) == x][1]})))
		tax_pos_tmp[!tax_pos_tmp %in% c('Proteobacteria', 'Firmicutes')] <- 'other'
		tax_pos_1[[run]] <- apply(matrix(tax_pos_tmp, ncol=2), 1, function(x){x1 <- x[order(x)]; paste0(x1, collapse='-')})
		
		tax_null_pos_tmp <- cbind(unlist(lapply(asv_names[null_pos $sp1], function(x){asvPhyVec[names(asvPhyVec) == x][1]})), unlist(lapply(asv_names[null_pos $sp2], function(x){asvPhyVec[names(asvPhyVec) == x][1]})))	
		tax_null_pos_tmp[! tax_null_pos_tmp %in% c('Proteobacteria', 'Firmicutes')] <- 'other'
		tax_null_pos_1[[run]] <- apply(matrix(tax_null_pos_tmp, ncol=2), 1, function(x){x1 <- x[order(x)]; paste0(x1, collapse='-')})
		
		tmpNullSp1 <- apply(real_pos, 1, function(x){x1 <- names(asv_depths)[asv_depths == x[3]]; return(x1[sample(1:length(x1), 1)])})
		tmpNullSp2 <- apply(real_pos, 1, function(x){x1 <- names(asv_depths)[asv_depths == x[4]]; return(x1[sample(1:length(x1), 1)])})
		
		tax_phylo_sub_tmp <- cbind(unlist(lapply(tmpNullSp1, function(x){asvPhyVec[names(asvPhyVec) == x][1]})), unlist(lapply(tmpNullSp2, function(x){asvPhyVec[names(asvPhyVec) == x][1]})))
		tax_phylo_sub_tmp[! tax_phylo_sub_tmp %in% c('Proteobacteria', 'Firmicutes')] <- 'other'
		tax_phylo_sub_1[[run]] <- apply(matrix(tax_phylo_sub_tmp, ncol=2), 1, function(x){x1 <- x[order(x)]; paste0(x1, collapse='-')})
				
	}
	
	tax_pos[[loc]] <- tax_pos_1 
	tax_null_pos[[loc]] <- tax_null_pos_1 
	tax_phylo_sub[[loc]] <- tax_phylo_sub_1 
}


#######################################################
# find proportion of pairs in each taxonomic category #
#######################################################

protNames <- names(asvPhyVec)[asvPhyVec == "Proteobacteria"]

phyPairs <- c('Firmicutes-Firmicutes', 'Firmicutes-other', 'Firmicutes-Proteobacteria', 'Proteobacteria-Proteobacteria', 'other-Proteobacteria', 'other-other')

tax_pos_prop <- list()
tax_null_prop <- list()
tab_phy_prop <- list()
tab_prot_exp <- list()
tab_prot_real <- list()

for(loc in 1:6){
	
	tmpProp <- list()
	tmpNull <- list()
	tmpPhy <- list()
	tmpExp <- list()
	
	for(run in 1:50){
		
		tp1 <- tax_pos[[loc]][[run]]
		tn1 <- tax_null_pos[[loc]][[run]]
		tphy <- tax_phylo_sub[[loc]][[run]]
		asv_names <- colnames(lg_list_cooc[[run]][[loc]])
		
		propProt <- length(asv_names[asv_names %in% protNames])/length(asv_names)
		
		tmpProp[[run]] <- unlist(lapply(phyPairs, function(x){sum(tp1 %in% x)}))/length(tp1)
		
		tmpNull[[run]] <- unlist(lapply(phyPairs, function(x){sum(tn1 %in% x)}))/length(tn1)
		
		tmpPhy[[run]] <- unlist(lapply(phyPairs, function(x){sum(tphy %in% x)}))/length(tphy)
		
		tmpExp[[run]] <- propProt* propProt
		
	}
	
	tax_pos_prop[[loc]] <- colMeans(do.call(rbind, tmpProp)[!apply(do.call(rbind, tmpProp), 1, function(x){'NaN' %in% x}),])
	tax_null_prop[[loc]] <- colMeans(do.call(rbind, tmpNull))
	tab_phy_prop[[loc]] <- colMeans(do.call(rbind, tmpPhy))
	tab_prot_exp[[loc]] <- unlist(tmpExp)
	tab_prot_real[[loc]] <- do.call(rbind, tmpProp)[,4][do.call(rbind, tmpProp)[,4] != 'NaN']
}


# t-test for number of Proteobacteria-Proteobacteria pairs vs null expectations

no_prot_prot_pairs <- list()

for(i in 1:6){
	
	no_prot_prot_pairs[[i]] <- t.test(tab_prot_real[[i]], tab_prot_exp[[i]])
	
}

tabS2_pairs <- unlist(lapply(no_prot_prot_pairs, function(x){x$p.value}))


#############################################################################################
# Figure 7: pairwise 'heatmap' of points representing proportions of pairs in each category #
#############################################################################################

phyPairs <- c('Firmicutes-Firmicutes', 'Firmicutes-other', 'Firmicutes-Proteobacteria', 'Proteobacteria-Proteobacteria', 'other-Proteobacteria', 'other-other')

par(xpd=NA)

par(mfrow=c(2, 3))

par(mar=c(7, 3, 4, 2))

for(i in 1:6){

	plot(c(0, 3), c(0, 3), type='n', axes=F, xlab='', ylab='', main= locLabels2[i])
	text(c(-0.25, -0.25, -0.25), c(0.5, 1.5, 2.5), c('P', 'F', 'O'))
	text(c(0.5, 1.5, 2.5), c(-0.25, -0.25, -0.25),  c('P', 'F', 'O'))

	cexTmp <- tax_pos_prop[[i]]*20
	
	points(0.5, 0.5, bg= phyCols[names(phyCols) == 'Proteobacteria'],  cex= cexTmp[4], pch=21)
	points(0.5, 1.5, bg='black', cex = cexTmp[3], pch=21)
	points(0.5, 2.5, bg='black', cex = cexTmp[5], pch=21)
	points(1.5, 0.5, bg='black', cex = cexTmp[3], pch=21)
	points(1.5, 1.5, bg= phyCols[names(phyCols) == 'Firmicutes'], pch=21, cex= cexTmp[1])
	points(1.5, 2.5, bg='black', cex = cexTmp[2], pch=21)	
	points(2.5, 0.5, bg='black', cex = cexTmp[5], pch=21)
	points(2.5, 1.5, bg='black', cex = cexTmp[2], pch=21)
	points(2.5, 2.5, bg='black', cex = cexTmp[6], pch=21)

}

legend(-8, 4.5, legend = c('Proteobacteria', 'Firmicutes', 'Other ASV'), fill=c(phyCols[names(phyCols) == 'Proteobacteria'], phyCols[names(phyCols) == 'Firmicutes'], phyCols[names(phyCols) == 'noPhy']), bty='n', ncol=3, cex=1.5)



###################
# modules #########
###################

# calculate modularity and find modules in bipartite host-bacteria network using metaComputeModule command from R package 'bipartite'

# single modularity statistic
mod_rare <- list()
mod_rare_rand <- list()

# output list object that includes hosts and bacteria
modules_info_rare <- list()	
modules_info_rare_rand <- list()	

# take 50 random subsamples of the full communities (15 hosts, 120 ASVs), calculate statistic and compare to randomized 

for(j in 1:50){
	
	lg_t1 <- lapply(lg_list3, function(x){x[sample(1:nrow(x), 15),]})

	lg_t2 <- list()

	for(i in 1:6){
		
		t1 <- lg_t1[[i]]
		
		lg_t2[[i]] <- t1[, apply(t1, 2, function(x){sum(x == 0)}) < 14]
		
	}

	lg_list_tmp <- lapply(lg_t2, function(x){x[,sample(1:ncol(x), 120)]})
	
	lg_rand_tmp <- lapply(lg_list_tmp, function(x){apply(x, 2, function(y){y[sample(1:length(y), length(y))]})})
	
	##########
	# loop over subsamples, find modularity and hosts/ASVs in modules
	modList <- list()
	webList <- list()

	modRand <- list()
	webRand <- list()

	for(i in 1:length(lg_list_tmp)){		
		
		lgTmp <- lg_list_tmp[[i]]	
		lgTmp[lgTmp > 0] <- 1		
		tmp <- DIRT_LPA_wb_plus(lg_list_tmp[[i]])		
		tmpLPA <- metaComputeModules(lgTmp, N=5)		
		webList[[i]] <- tmpLPA	
		modList[[i]] <- tmp$modularity
		
		lgTmpRand <- lg_rand_tmp[[i]]	
		lgTmpRand[lgTmpRand > 0] <- 1		
		tmpRand <- DIRT_LPA_wb_plus(lg_rand_tmp[[i]])	
		tmpLPARand <- metaComputeModules(lgTmpRand, N=5)		
		webRand[[i]] <- tmpLPARand	
		modRand[[i]] <- tmpRand$modularity
	}
		
	mod_rare[[j]] <- unlist(modList)
	mod_rare_rand[[j]] <- unlist(modRand)
	
	# output - within and between module connectivity - not validated for bipartite, may not make sense
	
	modInfo <- lapply(webList, function(x){listModuleInformation(x)})	
	modules_info_rare[[j]] <- modInfo	
	
	modInfoRand <- lapply(webRand, function(x){listModuleInformation(x)})	
	modules_info_rare_rand[[j]] <- modInfoRand	
	
	print(j)
	
}

# compare modularity

# visualize modularity stat ranges across runs

par(mfrow=c(1, 2))

boxplot(do.call(rbind, mod_rare))

boxplot(do.call(rbind, mod_rare_rand))

# format for t-test comparing distributions of modularity statistic within locations between real and null data

modMat <- do.call(rbind, mod_rare)
modRareMat <- do.call(rbind, mod_rare_rand)

mod_t <- list()

for(i in 1:6){
	
	mod_t[[i]] <- t.test(modMat[,i], modRareMat[,i])
	
}

tabS2_9 <- unlist(lapply(mod_t, function(x){x$p.value}))

# check for differences in modularity between locations 

modMat2 <- data.frame(cbind(as.vector(modMat), rep(lg_locs2, each=50)))
colnames(modMat2) <- c('mod', 'loc')

modularity_aov <- aov(mod ~ loc, modMat2)
summary(modularity_aov)
TukeyHSD(modularity_aov)

modularity_tukey <- TukeyHSD(modularity_aov)$loc

tabS3_6 <- apply(pairwiseNames, 1, function(x){modularity_tukey[grepl(x[1], rownames(modularity_tukey)) + grepl(x[2], rownames(modularity_tukey)) == 2,]})
colnames(tabS3_6) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_6) <- paste0('mean modularity ', rownames(tabS3_6))


#############

# compare number and size of modules

# metaComputeModules output is a list of length two: first entry has all hosts and ASVs, second entry is a list of the hosts-ASVs in each module. This code pulls out the number of modules per run, the mean number of ASVs in each module, and the mean number of hosts in each module

# real data
nMod <- list()
sizeMod <- list()
nHost <- list()

# randomized data
nModRand <- list()
sizeModRand <- list()
nHostRand <- list()

for(i in 1:6){
		
	nModtmp <- c()
	sizeModtmp <- c()
	nHostTmp <- c()
	
	for(j in 1:length(modules_info_rare)){
		
		t1 <- modules_info_rare[[j]][[i]][[2]]
		
		t2 <- t1[unlist(lapply(t1, function(x){length(x[[1]])})) > 1]
		
		nModtmp[j] <- length(t2)
		
		sizeModtmp[j] <- mean(unlist( lapply(t2, function(x){length(x[[2]])}) ) )
		
		nHostTmp[j] <- mean(unlist(lapply(t2, function(x){length(x[[1]])})))
		
	}
		
	nMod[[i]] <- nModtmp
	sizeMod[[i]] <- sizeModtmp
	nHost[[i]] <- nHostTmp
	
	nModtmp_rand <- c()
	sizeModtmp_rand <- c()
	nHostTmp_rand <- c()
	
	for(j in 1:length(modules_info_rare)){
		
		t1 <- modules_info_rare_rand[[j]][[i]][[2]]
		
		t2 <- t1[unlist(lapply(t1, function(x){length(x[[1]])})) > 1]
		
		nModtmp_rand[j] <- length(t2)
		
		sizeModtmp_rand[j] <- mean(unlist(lapply(t2, function(x){length(x[[2]])})))
		
		nHostTmp_rand[j] <- mean(unlist(lapply(t2, function(x){length(x[[1]])})))
		
	}
		
	nModRand[[i]] <- nModtmp_rand
	sizeModRand[[i]] <- sizeModtmp_rand
	nHostRand[[i]] <- nHostTmp_rand
	
}


nMod 
sizeMod 
nHost 

nModRand
sizeModRand 
nHostRand 

# compare number modules to random

n.mod_t <- list()

for(i in 1:6){
	
	n.mod_t[[i]] <- t.test(nMod[[i]], nModRand[[i]])
	
}

tabS2_10 <- unlist(lapply(n.mod_t, function(x){x$p.value}))

# compare number modules across sites in real data

n.modMat2 <- data.frame(cbind(unlist(nMod), rep(lg_locs2, each=50)))
colnames(n.modMat2) <- c('n.mod', 'loc')

nmod_aov <- aov(n.mod ~ loc, n.modMat2)
summary(nmod_aov)
TukeyHSD(nmod_aov)

nmod_tukey <- TukeyHSD(nmod_aov)$loc

tabS3_7 <- apply(pairwiseNames, 1, function(x){nmod_tukey[grepl(x[1], rownames(nmod_tukey)) + grepl(x[2], rownames(nmod_tukey)) == 2,]})
colnames(tabS3_7) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_7) <- paste0('mean number modules ', rownames(tabS3_7))


# compare number ASVs in modules to random

size.mod_t <- list()

for(i in 1:6){
	
	size.mod_t[[i]] <- t.test(sizeMod[[i]], sizeModRand[[i]])
	
}

tabS2_11 <- unlist(lapply(size.mod_t, function(x){x$p.value}))

# compare number ASVs in modules across locations

size.mod <- data.frame(cbind(unlist(sizeMod), rep(lg_locs2, each=50)))
colnames(size.mod) <- c('size.mod', 'loc')

sizemod_aov <- aov(size.mod ~ loc, size.mod)
summary(sizemod_aov)

sizemod_tukey <- TukeyHSD(sizemod_aov)$loc

tabS3_8 <- apply(pairwiseNames, 1, function(x){sizemod_tukey[grepl(x[1], rownames(sizemod_tukey)) + grepl(x[2], rownames(sizemod_tukey)) == 2,]})
colnames(tabS3_8) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_8) <- paste0('mean size modules ', rownames(tabS3_8))

# compare number hosts in modules to random

host.mod_t <- list()

for(i in 1:6){
	
	host.mod_t[[i]] <- t.test(nHost[[i]], nHostRand[[i]])
	
}

boxplot(nHost)

tabS2_12 <- unlist(lapply(host.mod_t, function(x){x$p.value}))

# compare number hosts across sites

n.host <- data.frame(cbind(unlist(nHost), rep(lg_locs2, each=50)))
colnames(n.host) <- c('nhost', 'loc')

nhost_aov <- aov(nhost ~ loc, n.host)
summary(nhost_aov)

hostmod_tukey <- TukeyHSD(nhost_aov)$loc

tabS3_9 <- apply(pairwiseNames, 1, function(x){hostmod_tukey[grepl(x[1], rownames(hostmod_tukey)) + grepl(x[2], rownames(hostmod_tukey)) == 2,]})
colnames(tabS3_9) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_9) <- paste0('mean host number per module ', rownames(tabS3_9))


# get mean patristic distance between ASVs within module and mean phylogenetic clustering of hosts within modules

phyASVmod <- list()
hostASVmod <- list()

phyASVmod_rand <- list()
hostASVmod_rand <- list()

for(i in 1:6){
	
	phyASVmod_tmp <- c()
	hostASVmod_tmp <- c()
	
	for(j in 1:length(modules_info_rare)){
		
		tot_asv_dist <- mean(distTips(keep.tip(lg_tree, modules_info_rare[[j]][[i]][[1]][[1]][[2]])))
		
		tot_host_vec <- modules_info_rare[[j]][[i]][[1]][[1]][[1]]
				
		t1 <- modules_info_rare[[j]][[i]][[2]]
		
		t2 <- t1[unlist(lapply(t1, function(x){length(x[[1]])})) > 1]
		
		t2_asv <- lapply(t2, function(x){x[[2]]})
		
		t2_host <- lapply(t2, function(x){x[[1]]})
		
		phyASVmod_tmp[j] <- mean(unlist(lapply(t2_asv, function(x){mean(distTips(keep.tip(lg_tree, x)))/tot_asv_dist})))
		
		hostASVmod_t1 <- prep_ses_2(tot_host_vec, t2_host, lg_host_tree)
		
		hostASVmod_tmp[j] <- mean(hostASVmod_t1)
		
	}
		
	phyASVmodRand_tmp <- c()
	hostASVmodRand_tmp <- c()
	
	for(j in 1:length(modules_info_rare)){
		
		tot_asv_dist <- mean(distTips(keep.tip(lg_tree, modules_info_rare_rand[[j]][[i]][[1]][[1]][[2]])))
		
		tot_host_vec <- modules_info_rare[[j]][[i]][[1]][[1]][[1]][as.numeric(modules_info_rare_rand[[j]][[i]][[1]][[1]][[1]])]
				
		t1 <- modules_info_rare_rand[[j]][[i]][[2]]
		
		t2 <- t1[unlist(lapply(t1, function(x){length(x[[1]])})) > 1]
		
		t2_asv <- lapply(t2, function(x){x[[2]]})
		
		t2_host <- lapply(t2, function(x){modules_info_rare[[j]][[i]][[1]][[1]][[1]][as.numeric(x[[1]])]})
		
		phyASVmodRand_tmp[j] <- mean(unlist(lapply(t2_asv, function(x){mean(distTips(keep.tip(lg_tree, x)))/tot_asv_dist})))
		
		hostASVmod_t1_rand <- prep_ses_2(tot_host_vec, t2_host, lg_host_tree)
		
		hostASVmodRand_tmp[j] <- mean(hostASVmod_t1_rand)
		
	}
				
	phyASVmod[[i]] <- phyASVmod_tmp
	hostASVmod[[i]] <- hostASVmod_tmp

	phyASVmod_rand[[i]] <- phyASVmodRand_tmp
	hostASVmod_rand[[i]] <- hostASVmodRand_tmp
	
	print(i)
	
}

# compare mean patristic distance of ASVs within modules to random

boxplot(phyASVmod)

asv_phy_t <- list()

for(i in 1:6){
	
	asv_phy_t[[i]] <- t.test(phyASVmod[[i]], phyASVmod_rand[[i]])
	
}

tabS2_13 <- unlist(lapply(asv_phy_t, function(x){x$p.value}))

# compare mean patristic distance across sites

asv_phy <- data.frame(cbind(unlist(phyASVmod), rep(lg_locs2, each=50)))
colnames(asv_phy) <- c('asvPhy', 'loc')

asv_phy_aov <- aov(asvPhy ~ loc, asv_phy)
summary(asv_phy_aov)
TukeyHSD(asv_phy_aov)

asv_phymod_tukey <- TukeyHSD(asv_phy_aov)$loc

tabS3_10 <- apply(pairwiseNames, 1, function(x){asv_phymod_tukey[grepl(x[1], rownames(asv_phymod_tukey)) + grepl(x[2], rownames(asv_phymod_tukey)) == 2,]})
colnames(tabS3_10) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_10) <- paste0('ASV phylo. clustering in modules ', rownames(tabS3_10))


# compare host phylogenetic clustering to random

boxplot(hostASVmod)

host_phy_t <- list()

for(i in 1:6){
	
	host_phy_t[[i]] <- t.test(hostASVmod[[i]], hostASVmod_rand[[i]])
	
}

tabS2_14 <- unlist(lapply(host_phy_t, function(x){x$p.value}))

# compare host phylogenetic clustering across sites

host_phy <- data.frame(cbind(unlist(hostASVmod), rep(lg_locs2, each=50)))
colnames(host_phy) <- c('hostPhy', 'loc')

host_phy_aov <- aov(hostPhy ~ loc, host_phy)
summary(host_phy_aov)
TukeyHSD(host_phy_aov)

host_phymod_tukey <- TukeyHSD(host_phy_aov)$loc

tabS3_11 <- apply(pairwiseNames, 1, function(x){host_phymod_tukey[grepl(x[1], rownames(host_phymod_tukey)) + grepl(x[2], rownames(host_phymod_tukey)) == 2,]})
colnames(tabS3_11) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_11) <- paste0('host phylo. clustering in modules ', rownames(tabS3_11))


# compile outputs into supplemental tables

tabS3_script2 <- rbind(tabS3_1, tabS3_2, tabS3_3, tabS3_4, tabS3_5, tabS3_6, tabS3_7, tabS3_8, tabS3_9, tabS3_10, tabS3_11)
write.csv(tabS3_script2, paste0(myPath, 'tabS3_script2.csv'))


tabS2_script2 <- list(tabS2_1, tabS2_2, tabS2_3, tabS2_4, tabS2_5, tabS2_6, tabS2_7, tabS2_8, tabS2_9, tabS2_10, tabS2_11, tabS2_12, tabS2_13, tabS2_14)

tabS2_script2_mat <- do.call(rbind, tabS2_script2)

rownames(tabS2_script2_mat) <- c('number pos vs. null', 'number neg vs. null', 'number pos vs. number neg', 'abundance pos pairs vs. null', 'phylo. relatedness of pos pairs vs. null', 'phylo. relatedness of pos pairs vs. size matched ASVs', 'phylo. relatedness of hosts vs. null', 'phylo. relatedness of hosts vs. size matched ASVs', 'modularity vs. null', 'number modules vs. null', 'size modules vs. null', 'mean number host per module vs. null', 'module ASV relatedness vs. null', 'module host relatedness vs. null')

colnames(tabS2_script2_mat) <- lg_locs2

write.csv(tabS2_script2_mat, paste0(myPath, 'tabS2_script2.csv'))
























