library(vegan)
arch = read.csv("archaea_abund.csv", row.names = "X.OTU.ID")
arch_m = read.csv("archaea_abund_match.csv", row.names = "X.OTU.ID")
vir_m = read.csv("pop_archvir_coverage_norm_log_raised_match.csv")
vir_m_b = vegdist(vir_m, method = "bray", binary = FALSE) #bray curtis
arch_m_b = vegdist(arch_m, method = "bray", binary = FALSE)#bray curtis
manarch_vir = mantel(arch_m_b, vir_m_b, method = "pearson", permutations = 999) #mantel
manarch_vir
arch_6 = read.csv("arch_6", row.names = "X.OTU.ID")
arch_2 = read.csv("arch_2", row.names = "X.OTU.ID")
arch_2_b = vegdist(arch_2, method = "bray", binary = FALSE)
arch_6_b = vegdist(arch_6, method = "bray", binary = FALSE)
manarch = mantel(arch_2_b, arch_6_b, method = "pearson", permutations = 999)
manarch
vir_6 = read.csv("vir_6")
vir_2 = read.csv("vir_2")
vir_6_b = vegdist(vir_6, method = "bray", binary = FALSE)
vir_2_b = vegdist(vir_2, method = "bray", binary = FALSE)
man_2 = mantel(arch_2_b, vir_2_b, method = "pearson", permutations = 999)
man_2
man_6 = mantel(arch_6_b, vir_6_b, method = "pearson", permutations = 999)
man_6
man_vir = mantel(vir_2_b, vir_6_b, method = "pearson", permutations = 999)
man_vir
env = read.csv("env_arch_all.csv", row.names = "X")

NMDSarch=metaMDS(arch, distance = "bray", k = 2, trymax = 1000, autotransform =TRUE, noshare = 0.1, expand = TRUE, trace = 1, plot = FALSE) #NMDS
attach(env)
plot(NMDSarch, type = "t", display = 'sites')
fit = envfit(NMDSarch~Temperature+Salinity+Oxygen+Chla+Depth+Feature, permutations=9999) #envfit
fit
NMDSarch
respo = ordisurf(NMDSarch~Oxygen, plot=FALSE) #responce surface
respt = ordisurf(NMDSarch~Temperature, plot=FALSE)
summary(respo)
summary(respt)
plot(respo)
plot(respt)
vir = read.csv("pop_archvir_coverage_norm_log_raised.csv")
NMDSvir=metaMDS(vir, distance = "bray", k = 2, trymax = 1000, autotransform =TRUE, noshare = 0.1, expand = TRUE, trace = 1, plot = FALSE) 
NMDSvir
plot(NMDSvir, type = "t", display = 'sites')
env_vir = read.csv("ETNP_ENV_VARxy.csv")
attach(env_vir)
fit_vir = envfit(NMDSvir~Temperature+Salinity+Oxygen+Chla+Depth+Feature, permutations=9999)
fit_vir
respo_vir = ordisurf(NMDSvir~Oxygen, plot=FALSE)
summary(respo_vir)
plot(NMDSvir, type = "t", display = 'sites')
plot(respo_vir)
respt_vir = ordisurf(NMDSvir~Temperature, plot=FALSE)
plot(respt_vir)
summary(respt_vir)
NMDSvir_m=metaMDS(vir_m, distance = "bray", k = 2, trymax = 1000, autotransform =TRUE, noshare = 0.1, expand = TRUE, trace = 1, plot = FALSE) 
NMDSvir_m
arch_abc = read.csv("archaea_abund_abc.csv", row.names = "X.OTU.ID")
attach(arch_abc)
fit_vir_host = envfit(NMDSvir_m~a+b+c+d+e+f+g+h+I+j+k+l+m+n+o+p+q+r, permutations = 9999)
plot(fit_vir_host)
fit_vir_host

cor_per = cor(arch_m, vir_m, use = 'everything', method = c('pearson')) #pearson correlation
cor_per
write.csv(cor_per, file = 'cor_per')
cor_spr = cor(arch_m, vir_m, use = 'everything', method = c('spearman')) #spearman correlation
cor_spr
write.csv(cor_spr, file = 'cor_spr')

richness<-specnumber(ab_raw) #calculate species richness
sumR <- summary(richness)
J <- H/log(specnumber(ab_raw)) #calculates peilou’s J
sumJ <- summary(J)
H<- diversity(ab_raw, index = "shannon") #calculates shannon’s index
sumH <- summary(H) 

library(RColorBrewer)
library(pvclust)
library(heatmap3)
PC <- read.csv("your_csv.csv", sep=",", header=TRUE) 
attach(PC)
PC_matrix <- data.matrix(PC)
Colv <- pvclust(PC_matrix, method.dist="max", method.hclust="average", nboot=100) #pvclust
Rowv <- pvclust((t(PC_matrix)), method.dist="max", method.hclust="average", nboot=100)
Col_dend <-as.dendrogram(Colv$hclust)
Row_dend <-as.dendrogram(Rowv$hclust)
hmcol<-colorRampPalette(c("#FFFFFF","#fff7bc","#ff0000"))(1000)
RowSideColors1<-rep("#FFFFFF",nrow(PC_matrix))
RowSideColors1[c(1,6,7)]<-"#DEEBF7"
RowSideColors1[c(2,3,8)]<-"#6BAED6"
RowSideColors1[c(4,9)]<-"#2171B5"
RowSideColors1[c(5,10)]<-"#08306B"
ColSideColors1<-rep("#FFFFFF",ncol(PC_matrix))
ColSideColors1[c(1,2,3,4,5,27,28,29,30,31,32,33,34,35,36,37,38,39,40,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)]<-"#DEEBF7"
ColSideColors1[c(40,41,42,43,22,23,24,25)]<-"#6BAED6"
ColSideColors1[c(21,26)]<-"#2171B5"
colGene <- cbind(RowSideColors1)
colGene2 <- cbind(ColSideColors1)
PC_heatmap <- heatmap3(PC_matrix, Rowv = Row_dend, Colv = Col_dend, RowSideColors=colGene, ColSideColors=colGene2, col=hmcol, scale="none") #heatmap
Colv 
Rowv