# clean workspace

rm(list=ls())

# read in libraries

library(ape)
library(vegan)
library(adephylo)
library(picante)
library(ade4)
library(adiv)
library(phangorn)
library(car)
library(maps)
library(phytools)
library(sf)

# 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, '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 (Figure S1)
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,]

# generate a presence absence version of dataset
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 parameters #
##############################

# assign colors to the most common phyla
phyCols <- c('#f98e09', '#bc3754', '#fcffa4',  '#57106e', '#000004')
names(phyCols) <- c("Proteobacteria", "Firmicutes", "Planctomycetes", "Actinobacteria", "noPhy")

# assign colors to each location
locCols <- c('#fde725', '#7ad151', '#22a884', '#2a788e', '#414487', '#440154')
names(locCols) <- lg_locs2

# label locations by 
locLabels <- c('Michigan', 'Georgia', 'Nicaragua highlands', 'Nicaragua dry forest', 'Nicaragua wet forest', 'Peru wet forest')

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

# save host species sampled per location for later plotting
locLabels2 <- paste0(locLabels, ' (', host_sp, ')')

# make a table of all pairwise comparisons 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)))

# ASV names in each most common phylum
namesProt <- names(asvPhyVec)[asvPhyVec == "Proteobacteria"]

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

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

# randomly sample equal numbers of hosts and ASVs, then save by writing out a dataset of six host x ASV matrices that can be used consistently throughout analyses

# find minimum number hosts

lapply(lg_list3, nrow)

lg_list4 <- lapply(lg_list3, function(x){x[sample(1:nrow(x), 15),]})

# this line removes per-location singleton ASVs

lg_list5 <- list()

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

# find minimum number ASVs

lapply(lg_list5, ncol)

# code to write out is commented out so script can be run without writing out new test files every time 

# lg_list6 <- lapply(lg_list5, function(x){x[,sample(1:ncol(x), 140)]})

# for(i in 1:6){	write.table(lg_list6[[i]], paste0(myPath, lg_locs2[i], '_hostxotu'), quote=F)}

# read files back in - not commented out so these files can be consistent throughout different runs of the script

lg_list7 <- list()

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

####################################################
# phylum breakdown #################################
####################################################

# get per ASV phylum breakdown for the subsampled dataset and plot

par(mar=c(1, 3, 3, 1))
par(mfrow=c(6, 1))
par(xpd=NA)

for(i in 1:6){

	t1 <- lg_list7[[i]]

	# get phyla of each ASV and make a table
	t2 <- apply(t1, 1, function(x){x1 <- colnames(t1)[x > 0]; table(asvPhyVec[names(asvPhyVec) %in% x1])})
	
	# plot the table as polygons that represent the proportion of phyla in each major group
	plot(c(0, nrow(t1)), c(0, 1), type='n', axes=F, xlab='', ylab='')
	
	for(j in 1:length(t2)){
		
		tmpTab <- t2[[j]]
		
		tt1 <- unlist(lapply(names(phyCols), function(x){x1 <- tmpTab[names(tmpTab) == x]; if(length(x1) > 0){return(x1)}else{return(0)}}))
		
		tt1[length(tt1)] <- tt1[length(tt1)] + sum(tmpTab)  - sum(tt1)
		
		tt2 <- c(0, cumsum(tt1/sum(tt1)))
	
		for(k in 1:(length(tt2) -1)){
			
			polygon(c(j-1, j-1, j, j), c(tt2[k], tt2[k+1], tt2[k+1], tt2[k]), col= phyCols[k], border=NA)
		}		
	}
	
	text(-1, 0.5, locLabels[i], srt=90)
}	


legend(0, 11.5, ncol=5, fill= phyCols, legend=names(phyCols), bty='n')



##################################################################
# code for range map figure - takes a long time to run ###########
##################################################################

# download IUCN scaled reptiles files: https://www.iucnredlist.org/resources/spatial-data-download

shp_path <- '~/path/to/SCALED_REPTILES/'

rep1 <- st_read(paste0(shp_path, "SCALED_REPTILES_PART1.shp"))

rep2 <- st_read(paste0(shp_path, "SCALED_REPTILES_PART2.shp"))

hostCex <- host_sp/max(host_sp)*5

# plot base map

quartz(height=7, width=6)
par(mar=c(1,0,1,0))
par(oma=c(0,0,0,0))

map("usa", xlim=c(-147, -35), ylim=c(-26,50), mar=c(0,0,0,0))
map("world","Mexico", add=T)
map("world","Guatemala", add=T)
map("world","Belize", add=T)
map("world","El Salvador", add=T)
map("world","Honduras", add=T)
map("world","Nicaragua", add=T)
map("world","Costa Rica", add=T)
map("world","Panama", add=T)
map("world","Colombia", add=T)
map("world","Ecuador", add=T)
map("world","Peru", add=T)
map("world","Venezuela", add=T)
map("world","Guyana", add=T)
map("world","Surinam", add=T)
map("world","French Guiana", add=T)
map("world","Brazil", add=T)
map("world","Bolivia", add=T)
map("world","Chile", add=T)
map("world","Cuba", add=T)
map("world","Jamaica", add=T)
map("world","Haiti", add=T)
map("world","Puerto Rico", add=T)
map("world","Dominican Republic", add=T)

for(i in 1:length(rep1$geometry)){
	plot(rep1$geometry[i], add=T, col=adjustcolor('dark blue', alpha.f=0.01), lty=0)
	print(i)
}

for(i in 1:length(rep2$geometry)){
	plot(rep2$geometry[i], add=T, col=adjustcolor('dark blue', alpha.f=0.01), lty=0)
	print(i)
}

# points for each location, colored by location colors and sized by the number of hosts sampled at each location - note that these locations are not exact, they are selected for plotting clarity

# Michigan
points(-83.9, 41.9, pch=21, bg='#4575b4', cex= hostCex[1])
# Georgia
points(-82.5, 31.4, pch=21, bg='#91bfdb', cex= hostCex[2])
# Nica wet forest
points(-84.3, 11.12, pch=21, bg='#fc8d59', cex= hostCex[5])
# Nica montane
points(-86, 12.8, pch=21, bg='#e0f3f8', cex= hostCex[3])
# Nica dry forest
points(-85, 11.5, pch=21, bg='#fee090', cex= hostCex[4])
# Peru wet forest
points(-71.8, -12.62, pch=21, bg='#d73027', cex= hostCex[6])

legend('bottomleft', ncol=1, legend = locLabels2, fill= locCols, bty='n')

# get the 
# number host species per location from the range maps

# get the exact centroid locations of each sampling area

MI <- c(-83.663, 42.302)
GA <- colMeans(lg_hosts[lg_hosts$sample.name %in% rownames(lg_list3[[2]]),9:10])
NIC_high <- colMeans(lg_hosts[lg_hosts$sample.name %in% rownames(lg_list3[[3]]),9:10])
NIC_dry <- apply(lg_hosts[lg_hosts$sample.name %in% rownames(lg_list3[[4]]),9:10], 2, function(x){mean(x[!is.na(x)])})
NIC_wet <- apply(lg_hosts[lg_hosts$sample.name %in% rownames(lg_list3[[5]]),9:10], 2, function(x){mean(x[!is.na(x)])})
PR <- apply(lg_hosts[lg_hosts$sample.name %in% rownames(lg_list3[[6]]),9:10], 2, function(x){mean(x[!is.na(x)])})

sampling_pts <- st_sfc(st_point(MI), st_point(GA), st_point(NIC_high),st_point(NIC_dry),st_point(NIC_wet),st_point(PR), crs=st_crs(rep1))

# this throws errors at some range maps - usually the ones with many polygons

# loop over each species polygon, see which overlap sampling locations

intersect_list <- list()

for(i in 1:length(rep1$geometry)){
	
	tmp_union <- st_make_valid(rep1$geometry[i])
	
	intersect_list[[i]] <- st_intersects(sampling_pts, tmp_union, sparse=F)
}


names(intersect_list) <- rep1$sci_name

intersect_list1.1 <- intersect_list[unlist(lapply(intersect_list, length)) == 6]

rep1_mat <- do.call(cbind, intersect_list1.1)

rowSums(rep1_mat)

# these are the layers that throw errors
rep1$sci_name[c(1325, 2504, 4054, 4260, 4668, 5788)]

# repeat for second set of polygons

intersect_list2 <- list()

for(i in 1:length(rep2$geometry)){
	
	tmp_union <- st_make_valid(rep2$geometry[i])
	
	intersect_list2[[i]] <- st_intersects(sampling_pts, tmp_union, sparse=F)
}

rep2$sci_name[1293]

names(intersect_list2) <- rep2$sci_name

intersect_list2.1 <- intersect_list2[unlist(lapply(intersect_list2, length)) == 6]

rep2_mat <- do.call(cbind, intersect_list2.1)

rowSums(rep2_mat)

############################################################
# get number of species, genera, and families per location #
############################################################

samp_sp <- unlist(lapply(lg_list3, function(x){xnm <- rownames(x); x1 <-  lg_hosts[lg_hosts$sample.name %in% xnm, 7]; length(unique(x1))}))

samp_gen <- unlist(lapply(lg_list3, function(x){xnm <- rownames(x); x1 <-  lg_hosts[lg_hosts$sample.name %in% xnm, 6]; length(unique(x1))}))

samp_fam <- unlist(lapply(lg_list3, function(x){xnm <- rownames(x); x1 <-  lg_hosts[lg_hosts$sample.name %in% xnm, 5]; length(unique(x1))}))


###########################################################################
# create a function to find phylogenetic clustering in the host community #
###########################################################################


prep_ses <- function(host_mat, host_tree){
	
	# read in a host x ASV matrix
	hostTmp <- colnames(t(host_mat))
		
	# get phylogenetic tree of hosts in the matrix
	# get genus and species of each host
	spTmp <- unlist(lapply(hostTmp, function(x){x1 <- lg_hosts[lg_hosts$sample.name == x,]; return(paste0(x1[6], '_', x1[7]))}));	
	# match to tree				
	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(host_mat)[, order(spTonini)]
	
	# make a vector of matched host names 
	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
	
	# get tree for samples
	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 get effect size
	ses <- ses.pd(hostTmp2, host_tree_clip)
	
	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)	
}


####################################################
# figure 1: specialization vs. prevalence vs. locs #
####################################################

asv_names <- unique(unlist(lapply(lg_list7, function(x){colnames(x)})))

# for each location, number of hosts in which each ASV occurs
abundance <- lapply(lg_list7, function(x){apply(x, 2, function(y){sum(y != 0)})})

# get host phylogenetic clustering relative to other host in the dataset for each ASV
host.spec <- lapply(lg_list7, function(x){prep_ses(x, lg_host_tree)})

# get mean prevalence of ASVs for reference in figure
meanAbAll <- mean(unlist(abundance))

# plot prevalence and specialization for each ASV, with points sized by number of locations in which they occur
quartz(width=8, height=5)
par(mfrow=c(1, 2))

plot(c(1, max(unlist(abundance))), c(-5, 2), type='n', xlab='mean prevalence', ylab='mean host specialization', main='Firmicutes')

lines(c(meanAbAll, meanAbAll), c(-5, 2), lwd=3, col='grey') 
lines(c(1, max(unlist(abundance))), c(0, 0), lwd=3, col='grey') 

for(i in 1:length(asv_names)){
	
	t1 <- asv_names[i]
	
	meanAb <- mean(unlist(lapply(abundance, function(x){x[names(x) == t1]})))
	
	meanSES <- mean(unlist(lapply(host.spec, function(x){x[names(x) == t1]})))
	
	tmpCex <- length(unlist(lapply(abundance, function(x){x[names(x) == t1]})))
	
	points(meanAb, meanSES, cex= 0.7, pch=4)

}


for(i in 1:length(asv_names)){
	
	t1 <- asv_names[i]
	
	meanAb <- mean(unlist(lapply(abundance, function(x){x[names(x) == t1]})))
	
	meanSES <- mean(unlist(lapply(host.spec, function(x){x[names(x) == t1]})))
	
	tmpCex <- length(unlist(lapply(abundance, function(x){x[names(x) == t1]})))
	
	tmpPhy <- unique(asvPhyVec[names(asvPhyVec) == t1])
	
	if(t1 %in% namesFirm){
		points(meanAb, meanSES, cex= tmpCex, pch=21, bg = phyCols[names(phyCols) == 'Firmicutes'])
	}
	

}


plot(c(1, max(unlist(abundance))), c(-5, 2), type='n', xlab='mean prevalence', ylab='mean host specialization', main='Proteobacteria')

lines(c(meanAbAll, meanAbAll), c(-5, 2), lwd=3, col='grey') 
lines(c(1, max(unlist(abundance))), c(0, 0), lwd=3, col='grey')

for(i in 1:length(asv_names)){
	
	t1 <- asv_names[i]
	
	meanAb <- mean(unlist(lapply(abundance, function(x){x[names(x) == t1]})))
	
	meanSES <- mean(unlist(lapply(host.spec, function(x){x[names(x) == t1]})))
	
	tmpCex <- length(unlist(lapply(abundance, function(x){x[names(x) == t1]})))
	
	points(meanAb, meanSES, cex= 0.7, pch=4)

}


for(i in 1:length(asv_names)){
	
	t1 <- asv_names[i]
	
	meanAb <- mean(unlist(lapply(abundance, function(x){x[names(x) == t1]})))
	
	meanSES <- mean(unlist(lapply(host.spec, function(x){x[names(x) == t1]})))
	
	tmpCex <- length(unlist(lapply(abundance, function(x){x[names(x) == t1]})))
	
	tmpPhy <- unique(asvPhyVec[names(asvPhyVec) == t1])
	
	if(t1 %in% namesProt){
		points(meanAb, meanSES, cex= tmpCex, pch=21, bg = phyCols[names(phyCols) == 'Proteobacteria'])
	}
	

}

#####################################################
# figure S4: specialization vs. prevalence per site #
#####################################################

# for each location, plot prevalence vs. abundance

abundance <- lapply(lg_list7, function(x){apply(x, 2, function(y){sum(y != 0)})})
abFirm <- lapply(lg_list7, function(x){x1 <- x[, colnames(x) %in% namesFirm]; apply(x1, 2, function(y){sum(y != 0)})})
abProt <- lapply(lg_list7, function(x){x1 <- x[, colnames(x) %in% namesProt]; apply(x1, 2, function(y){sum(y != 0)})})

medAb <- mean(unlist(abundance))
maxAb <-  max(unlist(abundance))

host.spec <- lapply(lg_list7, function(x){prep_ses(x, lg_host_tree)})

host.specFirm <- list()
for(i in 1:length(host.spec)){host.specFirm[[i]] <- host.spec[[i]][names(host.spec[[i]]) %in% namesFirm]}

host.specProt <- list()
for(i in 1:length(host.spec)){host.specProt[[i]] <- host.spec[[i]][names(host.spec[[i]]) %in% namesProt]}


par(mfrow=c(2, 3))
par(xpd=NA)

for(i in 1:6){
	
	plot(c(1.5, maxAb), c(-5, 2.5), xlab='', ylab='', type='n', main= locLabels2[i], axes=F, cex.axis=1.5)
	
	lines(c(1.5, maxAb), c(0, 0), lwd=3, col='black')
	text(0, 0, '1', pch=1.5)
	text(17, 0, maxAb, pch=1.5)
	lines(c(3.5, 3.5), c(-5, 2.5), lwd=3, col='black')
	text(2.5, 2.5, '2.5', pch=1.5)
	text(2.5, -5, '-5', pch=1.5)
	
	text(-2.25, 0, 'host specialization', srt=90, cex=1.2)
	
	text(8, -6, 'prevalence', cex=1.2)
	
	
	points(abundance[[i]], host.spec[[i]], pch=4, cex=1.5)
	
	points(abundance[[i]][names(abundance[[i]]) %in% namesFirm], host.spec[[i]][names(abundance[[i]]) %in% namesFirm], pch=21, bg=phyCols[names(phyCols) == 'Firmicutes'], cex=2)
	
	points(abundance[[i]][names(abundance[[i]]) %in% namesProt], host.spec[[i]][names(abundance[[i]]) %in% namesProt], pch=23, bg=phyCols[names(phyCols) == 'Proteobacteria'], cex=2)

	
}

legend(-38, 5.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)


##########################################
# PERMANOVA + Figure 2 ###################
##########################################

# test whether microbial communities cluster by location or host family

lg_hosts2 <- lg_hosts[lg_hosts$sample.name %in% rownames(lg_reads_rare),]

# get families with more than seven samples
lg_fams <- names(table(lg_hosts2$family))[table(lg_hosts2$family) > 7]

lg_hosts3 <- lg_hosts2[lg_hosts2$family %in% lg_fams,]

# make host x ASV matrix with only the hosts in the focal families
lg_perm <- lg_reads_rare[rownames(lg_reads_rare) %in% lg_hosts3$sample.name,]
lg_perm2 <- lg_perm[,colSums(lg_perm) > 0]

# make a vector of host family and location for each host
lg_fam_fac <- as.factor(unlist(lapply(rownames(lg_perm2), function(x){lg_hosts3$family[lg_hosts3$sample.name == x]})))
lg_loc_fac <- as.factor(unlist(lapply(rownames(lg_perm2), function(x){lg_hosts3$locality[lg_hosts3$sample.name == x]})))

# run permanova
both_perm <- adonis2(lg_perm2 ~ lg_loc_fac + lg_fam_fac + lg_loc_fac:lg_fam_fac, permutations=100000)

# Repeat for only Firmicutes ASVs

lg_perm_firm <- lg_reads_rare[rownames(lg_reads_rare) %in% lg_hosts3$sample.name, colnames(lg_reads_rare) %in% namesFirm]
lg_perm2_firm <- lg_perm_firm[,colSums(lg_perm_firm) > 0]
lg_perm3_firm <- lg_perm2_firm[rowSums(lg_perm2_firm) > 0,]
lg_loc_fac_firm <- lg_loc_fac[rowSums(lg_perm2_firm) > 0]
lg_fam_fac_firm <- lg_loc_fac[rowSums(lg_perm2_firm) > 0]

perm_firm <- adonis2(lg_perm3_firm ~ lg_fam_fac_firm + lg_loc_fac_firm + lg_fam_fac_firm:lg_loc_fac_firm, permutations=100000)

# Repeat for only Proteobacteria ASVs ASVs

lg_perm_prot <- lg_reads_rare[rownames(lg_reads_rare) %in% lg_hosts3$sample.name, colnames(lg_reads_rare) %in% namesProt]
lg_perm2_prot <- lg_perm_prot[,colSums(lg_perm_prot) > 0]
lg_perm3_prot <- lg_perm2_prot[rowSums(lg_perm2_prot) > 0,]
lg_loc_fac_prot <- lg_loc_fac[rowSums(lg_perm2_prot) > 0]
lg_fam_fac_prot <- lg_loc_fac[rowSums(lg_perm2_prot) > 0]

perm_prot <- adonis2(lg_perm3_prot ~ lg_loc_fac_prot + lg_fam_fac_prot + lg_loc_fac_prot:lg_fam_fac_prot, permutations=100000)

# plot distance between groups on XY plane using nmds

# assign point character by family (grouping snakes together) and color points by location

lg_fam_fac <- as.character(lg_fam_fac)
lg_fam_fac[lg_fam_fac == 'Colubridae'] <- 24
lg_fam_fac[lg_fam_fac == 'Dactyloidae'] <- 21
lg_fam_fac[lg_fam_fac == 'Dipsadinae'] <- 24
lg_fam_fac[lg_fam_fac == 'Gekkonidae'] <- 23
lg_fam_fac[lg_fam_fac == 'Gymnopthalmidae'] <- 22
lg_fam_fac[lg_fam_fac == 'Natricinae'] <- 24
lg_fam_fac[lg_fam_fac == 'Teiidae'] <- 22
lg_fam_fac[lg_fam_fac == 'Tropiduridae'] <- 25
lg_fam_fac[lg_fam_fac == 'Vipiridae'] <- 24

lg_loc_fac <- as.character(lg_loc_fac)
lg_loc_fac[lg_loc_fac == 'Michigan'] <- locCols[1]
lg_loc_fac[lg_loc_fac == 'Georgia'] <- locCols[2]
lg_loc_fac[lg_loc_fac == 'Nicaragua highlands'] <- locCols[3]
lg_loc_fac[lg_loc_fac == 'Nicaragua dry forest'] <- locCols[4]
lg_loc_fac[lg_loc_fac == 'Nicaragua wet forest'] <- locCols[5]
lg_loc_fac[lg_loc_fac == 'Peru wet forest'] <- locCols[6]


metaAll <- plot(metaMDS(lg_perm2))

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

quartz(width=6, height=6)
par(xpd=NA)
par(mar=c(5, 5, 8, 3))

plot(metaAll$sites[,1], metaAll$sites[,2], type='n', xlab='NMDS 1', ylab='NMDS 2', axes=F)
axis(2, seq(-0.1, 0.2, by=0.1), seq(-0.1, 0.2, by=0.1), las=2)
axis(1, seq(-0.15, 0.15, by=0.1), seq(-0.15, 0.15, by=0.1))

points(metaAll$sites[,1], metaAll$sites[,2], pch=as.numeric(lg_fam_fac), bg= lg_loc_fac, cex=1.5)

legend(-0.2, 0.27, legend= c('Michigan', 'Georgia', 'Nicaragua highlands', 'Nicaragua dry forest', 'Nicaragua wet forest', 'Peru wet forest'), fill= locCols, bty='n', ncol=2)

legend(0.1, 0.2, legend= c('Dactyloidae', 'Teiidae', 'Gekkonidae', 'snakes', 'Tropiduridae'), pch=c(21, 22, 23, 24, 25), bty='n', ncol=1)


###############################################
# species accumulation curve + Figure 3 #######
###############################################

pool_firm <- lapply(lg_list7, function(x){x1 <- x[,colnames(x) %in% namesFirm]; poolaccum(x1)})

pool_prot <- lapply(lg_list7, function(x){x1 <- x[,colnames(x) %in% namesProt]; poolaccum(x1)})

poolAcc_7 <- lapply(lg_list7, poolaccum)


quartz(width=6, height=3)
par(mfrow=c(1, 3))
par(xpd=NA)
par(las=2)
par(mar=c(5, 5, 8, 2))

plot(c(1, 16), c(0, 250), type='n', xlab='hosts', ylab='bootstrap number species', axes=F, main='All', cex=1.5)
axis(1, seq(0, 16, by=4), seq(0, 16, by=4), las=1)
axis(2, seq(0, 250, by=50), seq(0, 250, by=50), las=2)

text(-1, 300, 'A)', cex=1.5)

for(i in 1:6){	
	
	lines(summary(poolAcc_7[[i]])$boot[,2], lwd=3, col=locCols[i])
	lines(summary(poolAcc_7[[i]])$boot[,3], lwd=1, col=locCols[i])
	lines(summary(poolAcc_7[[i]])$boot[,4], lwd=1, col=locCols[i])
}

legend(-2, 480, legend= locLabels2, fill= locCols, bty='n', ncol=3, cex=1.2)


plot(c(1, 16), c(0, 100), type='n', xlab='hosts', ylab='bootstrap number species', axes=F, main='Firmicutes')
axis(1, seq(0, 16, by=4), seq(0, 16, by=4), las=1)
axis(2, seq(0, 100, by=25), seq(0, 100, by=25))

text(-1, 120, 'B)', cex=1.5)

for(i in 1:6){	
	
	lines(summary(pool_firm[[i]])$boot[,2], lwd=3, col=locCols[i])
	lines(summary(pool_firm[[i]])$boot[,3], lwd=1, col=locCols[i])
	lines(summary(pool_firm[[i]])$boot[,4], lwd=1, col=locCols[i])
}


plot(c(1, 16), c(0, 100), type='n', xlab='hosts', ylab='bootstrap number species', axes=F, main='Proteobacteria')
axis(1, seq(0, 16, by=4), seq(0, 16, by=4), las=1)
axis(2, seq(0, 100, by=25), seq(0, 100, by=25))

text(-1, 120, 'C)', cex=1.5)

for(i in 1:6){	
	
	lines(summary(pool_prot[[i]])$boot[,2], lwd=3, col=locCols[i])
	lines(summary(pool_prot[[i]])$boot[,3], lwd=1, col=locCols[i])
	lines(summary(pool_prot[[i]])$boot[,4], lwd=1, col=locCols[i])
}


#################################################
# Figure 4: richness and specializtion boxplots #
#################################################

# get shannon index for full dataset, Firmicutes only, and Proteobacteria only

shannon <- lapply(lg_list7, function(x){vegan::diversity(x)})
shannonFirm <- lapply(lg_list7, function(x){x1 <- x[, colnames(x) %in% namesFirm]; vegan::diversity(x1)})
shannonProt <- lapply(lg_list7, function(x){x1 <- x[, colnames(x) %in% namesProt]; vegan::diversity(x1)})

# prepare shannon values to write out in supplemental table
tabS2_1 <- unlist(lapply(shannon, mean))
tabS2_2 <- unlist(lapply(shannonFirm, mean))
tabS2_3 <- unlist(lapply(shannonProt, mean))

# run ANOVA between sites to test for differences 
shannonMat <- data.frame(cbind(unlist(shannon), rep(lg_locs2, each= nrow(lg_list7[[1]]))))
colnames(shannonMat) <- c('shannon', 'loc')

shan_aov <- aov(shannon ~ loc, shannonMat)
summary(shan_aov)

# run a Tukey test to find signifcant pairwise comparisons, then prepare to write out in supplmental table 
TukeyHSD(shan_aov)
shan_aov_tukey <- TukeyHSD(shan_aov)$loc


tabS3_0 <- apply(pairwiseNames, 1, function(x){shan_aov_tukey[grepl(x[1], rownames(shan_aov_tukey)) + grepl(x[2], rownames(shan_aov_tukey)) == 2,]})
colnames(tabS3_0) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_0) <- paste0('All Shannon ', rownames(tabS3_0))


# Firmicutes

shannonMatFirm <- data.frame(cbind(unlist(shannonFirm), rep(lg_locs2, each=nrow(lg_list7[[1]]))))
colnames(shannonMatFirm) <- c('shannon', 'loc')

firm_aov <- aov(shannon ~ loc, shannonMatFirm)
summary(firm_aov)
TukeyHSD(firm_aov)
firm_aov_tukey <- TukeyHSD(firm_aov)$loc

tabS3_1 <- apply(pairwiseNames, 1, function(x){firm_aov_tukey[grepl(x[1], rownames(firm_aov_tukey)) + grepl(x[2], rownames(firm_aov_tukey)) == 2,]})
colnames(tabS3_1) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_1) <- paste0('Firmicutes Shannon ', rownames(tabS3_1))

# Proteobacteria

shannonMatProt <- data.frame(cbind(unlist(shannonProt), rep(lg_locs2, each=nrow(lg_list7[[1]]))))
colnames(shannonMatProt) <- c('shannon', 'loc')

prot_aov <- aov(shannon ~ loc, shannonMatProt)
summary(prot_aov)
TukeyHSD(prot_aov)
prot_aov_tukey <- TukeyHSD(prot_aov)$loc

tabS3_2 <- apply(pairwiseNames, 1, function(x){prot_aov_tukey[grepl(x[1], rownames(prot_aov_tukey)) + grepl(x[2], rownames(prot_aov_tukey)) == 2,]})
colnames(tabS3_2) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_2) <- paste0('Proteobacteria Shannon ', rownames(tabS3_2))


# Repeat with Faith's PD - phylogenetically informed diversity metric

faiths <- lapply(lg_list7, function(x){pd(x, lg_tree, include.root=F)})
names(faiths) <- lg_locs2	
faith_pd <- lapply(faiths, function(x){x$PD})

faithsFirm <- lapply(lg_list7, function(x){x1 <- x[, colnames(x) %in% namesFirm]; pd(x1, lg_tree, include.root=F)})
names(faithsFirm) <- lg_locs2	
faith_pdFirm <- list()
for(i in 1:6){x1 <- as.numeric(faithsFirm[[i]]$PD); x1[is.na(x1)] <- 0; faith_pdFirm[[i]] <- x1}

faithsProt <- lapply(lg_list7, function(x){x1 <- x[, colnames(x) %in% namesProt]; pd(x1, lg_tree, include.root=F)})
names(faithsProt) <- lg_locs2	
faith_pdProt <- list()
for(i in 1:6){x1 <- as.numeric(faithsProt[[i]]$PD); x1[is.na(x1)] <- 0; faith_pdProt[[i]] <- x1}

tabS2_4 <- unlist(lapply(faith_pd, mean))
tabS2_5 <- unlist(lapply(faith_pdFirm, mean))
tabS2_6 <- unlist(lapply(faith_pdProt, mean))

# all 
faithMat <- data.frame(cbind(unlist(faith_pd), rep(lg_locs2, each=16)))
colnames(faithMat) <- c('faith', 'loc')

all_aov_faith <- aov(faith ~ loc, faithMat)
summary(all_aov_faith)
TukeyHSD(all_aov_faith)

all_faith_tukey <- TukeyHSD(all_aov_faith)$loc

tabS3_3a <- apply(pairwiseNames, 1, function(x){all_faith_tukey[grepl(x[1], rownames(all_faith_tukey)) + grepl(x[2], rownames(all_faith_tukey)) == 2,]})
colnames(tabS3_3a) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_3a) <- paste0('Firmicutes Faith ', rownames(tabS3_3a))

# Firmicutes

faithMatFirm <- data.frame(cbind(unlist(faithsFirm), rep(lg_locs2, each=16)))
colnames(faithMatFirm) <- c('faith', 'loc')

lapply(faith_pdFirm, mean)

firm_aov_faith <- aov(faith ~ loc, faithMatFirm)
summary(firm_aov_faith)
TukeyHSD(firm_aov_faith)

firm_faith_tukey <- TukeyHSD(firm_aov_faith)$loc

tabS3_3 <- apply(pairwiseNames, 1, function(x){firm_faith_tukey[grepl(x[1], rownames(firm_faith_tukey)) + grepl(x[2], rownames(firm_faith_tukey)) == 2,]})
colnames(tabS3_3) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_3) <- paste0('Firmicutes Faith ', rownames(tabS3_3))

# Proteobacteria

faithMatProt <- data.frame(cbind(unlist(faithsProt), rep(lg_locs2, each=16)))
colnames(faithMatProt) <- c('faith', 'loc')

lapply(faith_pdProt, mean)

prot_aov_faith <- aov(faith ~ loc, faithMatProt)
summary(prot_aov_faith)
TukeyHSD(prot_aov_faith)

prot_faith_tukey <- TukeyHSD(prot_aov_faith)$loc

tabS3_4 <- apply(pairwiseNames, 1, function(x){prot_faith_tukey[grepl(x[1], rownames(prot_faith_tukey)) + grepl(x[2], rownames(prot_faith_tukey)) == 2,]})
colnames(tabS3_4) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_4) <- paste0('Proteobacteria Faith ', rownames(tabS3_4))


# host specificity: ses.pd for each location

host.spec <- lapply(lg_list7, function(x){prep_ses(x, lg_host_tree)})

host.specFirm <- list()
for(i in 1:length(host.spec)){host.specFirm[[i]] <- host.spec[[i]][names(host.spec[[i]]) %in% namesFirm]}

host.specProt <- list()
for(i in 1:length(host.spec)){host.specProt[[i]] <- host.spec[[i]][names(host.spec[[i]]) %in% namesProt]}

tabS2_7 <- unlist(lapply(host.spec, mean))
tabS2_8 <- unlist(lapply(host.specFirm, mean))
tabS2_9 <- unlist(lapply(host.specProt, mean))

# all 
host.specMat <- data.frame(cbind(as.numeric(unlist(host.spec)), rep(lg_locs2, each=unlist(lapply(host.spec, length)))))
colnames(host.specMat) <- c('host.spec', 'loc')

aov_host.spec <- aov(host.spec ~ loc, host.specMat)
summary(aov_host.spec)
TukeyHSD(aov_host.spec)

all_host.spec_tukey <- TukeyHSD(aov_host.spec)$loc

tabS3_5 <- apply(pairwiseNames, 1, function(x){all_host.spec_tukey[grepl(x[1], rownames(all_host.spec_tukey)) + grepl(x[2], rownames(all_host.spec_tukey)) == 2,]})
colnames(tabS3_5) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_5) <- paste0('all host.spec ', rownames(tabS3_5))

# Firmicutes

host.specMatFirm <-  data.frame(cbind(unlist(host.specFirm), rep(lg_locs2, unlist(lapply(host.specFirm, length)))))
colnames(host.specMatFirm) <- c('host.spec', 'loc')

firm_aov_host.spec <- aov(host.spec ~ loc, host.specMatFirm)
summary(firm_aov_host.spec)
TukeyHSD(firm_aov_host.spec)

firm_host.spec_tukey <- TukeyHSD(firm_aov_host.spec)$loc

tabS3_6 <- apply(pairwiseNames, 1, function(x){firm_host.spec_tukey[grepl(x[1], rownames(firm_host.spec_tukey)) + grepl(x[2], rownames(firm_host.spec_tukey)) == 2,]})
colnames(tabS3_6) <- paste0(pairwiseNames[,1], '-', pairwiseNames[,2])
rownames(tabS3_6) <- paste0('Firm host.spec ', rownames(tabS3_6))

# Proteobacteria

host.specMatProt <- data.frame(cbind(unlist(host.specProt), rep(lg_locs2, unlist(lapply(host.specProt, length)))))
colnames(host.specMatProt) <- c('host.spec', 'loc')


prot_aov_host.spec <- aov(host.spec ~ loc, host.specMatProt)
summary(prot_aov_host.spec)
TukeyHSD(prot_aov_host.spec)

prot_host.spec_tukey <- TukeyHSD(prot_aov_host.spec)$loc

#################################################################################
# Figure 4: boxplots of Faith's diversity and host specificity across locations #
#################################################################################

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


plot(c(0, 24), c(0, 12), type='n', axes=F, xlab='', ylab='Faiths Phylogenetic Diversity')
axis(2, seq(0, 12, by=3), seq(0, 12, by=3))
text(seq(2, 23, by=4), rep(-1.5, 6), labels= locLabels2, srt=90, pos=2)
text(-4, 14, 'A)')

boxplot(faith_pd, add=T, at=seq(1, 22, by=4), col= 'grey',  border = 'black', axes=F)
boxplot(faith_pdFirm, add=T, at=seq(2, 23, by=4), col= phyCols[names(phyCols) == 'Firmicutes'],  axes=F)
boxplot(faith_pdProt, add=T, at=seq(3, 24, by=4), col = phyCols[names(phyCols) == 'Proteobacteria'],  axes=F)


plot(c(0, 24), c(-7, 2), type='n', axes=F, xlab='', ylab='Host specialization index')
axis(2, -7:2, -7:2)
text(seq(2, 23, by=4), rep(-8, 6), labels= locLabels2, srt=90, pos=2)
text(-4, 3.7, 'B)')

boxplot(host.spec, add=T, at=seq(1, 22, by=4), col= 'grey',  border = 'black', axes=F)
boxplot(host.specFirm, add=T, at=seq(2, 23, by=4), col= phyCols[names(phyCols) == 'Firmicutes'],  axes=F)
boxplot(host.specProt, add=T, at=seq(3, 24, by=4), col=  phyCols[names(phyCols) == 'Proteobacteria'],  axes=F)


legend(-40, 5.5, legend=c('all', 'Firmicutes', 'Proteobacteria'), fill = c('grey', phyCols[names(phyCols) == 'Firmicutes'], phyCols[names(phyCols) == 'Proteobacteria']), ncol=3, bty='n')


################################
# write out supplemental files #
################################


tabS2_script2 <- list(tabS2_1, tabS2_2, tabS2_3, tabS2_4, tabS2_5, tabS2_6, tabS2_7, tabS2_8, tabS2_9)

tabS2_script2_mat <- do.call(rbind, tabS2_script2)
rownames(tabS2_script2_mat) <- c('Shannon all', 'Shannon Firm', 'Shannon Prot', 'Faith all', 'Faith Firm', 'Faith Prot', 'Host spec all', 'Host spec Firm', 'Host spec Prot')


write.csv(tabS2_script2_mat, paste0(myPath, '/tabS2_script1.csv'))


tabS3_script2 <- rbind(tabS3_0, tabS3_1, tabS3_2, tabS3_3, tabS3_3a, tabS3_4, tabS3_5, tabS3_6)
write.csv(tabS3_script2, paste0(myPath, 'tabS3_script1.csv'))
















