#Type 3 ANOVA conducted on each individual taxa in taxa x sample table to determine whether correlated with treatment # in this code looked at effect of N+host+N*host. library(nlme) contrasts=c("contr.helmert","contr.poly");options("contrasts") map=read.csv("Desktop/Niwot/F_manuscript/edits/analyses/NiwotDGF_map3.csv",header=TRUE) map$Plot<-factor(map$Plot) map$Run<-factor(map$Run) Adata=read.csv("Desktop/Niwot/F_manuscript/edits/analyses/beta/reduced_nosingles_rarefied_557_99_L2.csv",header=TRUE) dim(Adata) #just take abundance data A<-Adata[,1:37] A2<-t(A) colnames(A2)<-Adata$Taxon #need dataset without OTUs occoring in fewer than 2 plots ind<-which(colSums(ifelse(A2>0,1,0))>2) A3<-A2[,ind] # recalculate relative abundance now that single and double occurances removed A4<-A3/rowSums(A3)*100 A5<-as.data.frame(A4) # include factor info in OTU abundance data.frame A5$host<-map$Host A5$Run<-factor(map$Run) A5$N<-map$Nitrogen A5$rem<-map$Deschampsia A5$treatment<-map$Treatment A5$treatment2<-map$Treatment2 A5$plot<-factor(map$Plot) ############################################################################################################# #exclude deschampsia removal plots ############################################################################################################# ind<-which(A5$rem=="Control") A6<-A5[ind,] ind<-which(map$Deschampsia=="Control") map1<-map[ind,] #need dataset without OTUs occoring in fewer than 2 plots n<-dim(A6)[2]-7 ind<-which(colSums(ifelse(A6[,1:n]>0,1,0))>2) A7<-A6[,ind] # recalculate relative abundance now that single and double occurances removed A8<-A7/rowSums(A7)*100 A9<-as.data.frame(A8) A9$Run<-factor(map1$Run) A9$host<-map1$Host A9$N<-map1$Nitrogen A9$treatment2<-map1$Treatment2 A9$plot<-factor(map1$Plot) ############################################################################################################# Taxa<-as.vector(colnames(A8)) n<-length(Taxa) N.p<-vector("numeric",n) N.F<-vector("numeric",n) host.p<-vector("numeric",n) host.F<-vector("numeric",n) Nxhost.p<-vector("numeric",n) Nxhost.F<-vector("numeric",n) for(i in 1:n) { A9$temptaxa<-A9[,i] a0<-lme(temptaxa~N+host+N*host,random=~1|plot,data=A9,method="REML");a0 anova0<-anova(a0,type="marginal") N.F[i]<-anova0[2,3] N.p[i]<-anova0[2,4] host.F[i]<-anova0[3,3] host.p[i]<-anova0[3,4] Nxhost.F[i]<-anova0[4,3] Nxhost.p[i]<-anova0[4,4] } N.FDRp<-p.adjust(N.p,method="fdr") host.FDRp<-p.adjust(host.p,method="fdr") Nxhost.FDRp<-p.adjust(Nxhost.p,method="fdr") #Average values for fert X,N and rem X,D to compare l=length(A9[1,])-6 treat2<-aggregate.data.frame(A9[,1:l],by=list(A9$treatment2),mean) treat2a<-treat2[,2:(l+1)] treat2b<-data.frame(t(treat2a)) colnames(treat2b)<-treat2[,1] treat2_DC.av<-treat2b$DesC treat2_DN.av<-treat2b$DesN treat2_GC.av<-treat2b$GeumC treat2_GN.av<-treat2b$GeumN #prepare stats output file p_OTU<-data.frame(Taxa,N.F,N.p,N.FDRp,host.F,host.p,host.FDRp,Nxhost.F,Nxhost.p,Nxhost.FDRp,treat2_DC.av,treat2_DN.av,treat2_GC.av,treat2_GN.av) #provide list of significant correlations sig<-which(N.p<=0.05|host.p<=0.05|Nxhost.p<=0.05) SigTaxa<-p_OTU[sig,] SigTaxa #put all these in an output file output<-as.matrix(SigTaxa) write(colnames(output),file="Desktop/Niwot/F_manuscript/edits/analyses/otu_anova/ANOVA_HxN.csv",ncolumns=14,sep = ",", append=TRUE) output<-t(output) write(output,file="Desktop/Niwot/F_manuscript/edits/analyses/otu_anova/ANOVA_HxN.csv",ncolumns=14,sep = ",",append=TRUE)