######################################### #### Code for real data # ######################################### # 1.Enter data with first p columns representing OTUs/ Taxa abundance counts # and last column representing the grouping variable. Example: # OTU1 OTU2 ... OTUp grp # c11 c12 ... c1p Grp1 # c21 c22 ... c2p Grp2 # .... # 2.Enter taxonomy matrix for the OTUs with each row for a single OTU # in same order as columns in the real data. Example below: # OTU_ID Taxonomy # OTU1 String1 # OTU2 String2 # ... ########################################### real.data = read.table("real.data.txt", header=T, quote="\"") colnames(real.data)[dim(real.data)[2]]="grp" tax.id = read.table("taxonomy.txt", header=T, quote="\"") par1_new=dim(real.data)[2]-1 ancom.detect=function(otu_data,n_otu,alpha){ logratio.mat=matrix(NA,nr=n_otu,nc=n_otu) for(i in 1:(n_otu-1)){ for(j in (i+1):n_otu){ data.pair=otu_data[,c(i,j,n_otu+1)] lr=log((0.001+as.numeric(data.pair[,1]))/(0.001+as.numeric(data.pair[,2]))) #lr=log((data.pair[,1])/(data.pair[,2])) logratio.mat[i,j]=wilcox.test(lr[data.pair\$grp==unique(data.pair\$grp)[1]],lr[data.pair\$grp==unique(data.pair\$grp)[2]])\$p.value } } ind <- lower.tri(logratio.mat) logratio.mat[ind] <- t(logratio.mat)[ind] logratio.mat[which(is.finite(logratio.mat)==FALSE)]=NA a=logratio.mat[upper.tri(logratio.mat,diag=F)==T] b=matrix(0,nc=n_otu,nr=n_otu) b[upper.tri(b)==T]=p.adjust(a, method = "BH") diag(b)=NA ind.1 <- lower.tri(b) b[ind.1] <- t(b)[ind.1] # W.original=apply(logratio.mat,1,function(x){ # subp=length(which(x<0.05)) # }) W=apply(b,1,function(x){ subp=length(which(x=0.10){ c.start=max(W.detected)/par1_new cutoff=c.start-c(0.05,0.10,0.15,0.20,0.25) prop_cut=rep(0,length(cutoff)) for(cut in 1:length(cutoff)){ prop_cut[cut]=length(which(W.detected>=par1_new*cutoff[cut]))/length(W.detected) } del=rep(0,length(cutoff)-1) for(i in 1:(length(cutoff)-1)){ del[i]=abs(prop_cut[i]-prop_cut[i+1]) } if(del[1]<0.02&del[2]<0.02&del[3]<0.02){nu=cutoff[1] }else if(del[1]>=0.02&del[2]<0.02&del[3]<0.02){nu=cutoff[2] }else if(del[2]>=0.02&del[3]<0.02&del[4]<0.02){nu=cutoff[3] }else{nu=cutoff[4]} up_point=min(W.detected[which(W.detected>=nu*par1_new)]) W.detected[W.detected>=up_point]=99999 W.detected[W.detected=par1*cutoff[cut]))/length(detected.diff[,4]) } del=abs(diff(prop_cut, differences = 1)) if(del[1]<0.02&del[2]<0.02&del[3]<0.02){nu=cutoff[1] }else if(del[1]>=0.02&del[2]<0.02&del[3]<0.02){nu=cutoff[2] }else if(del[2]>=0.02&del[3]<0.02&del[4]<0.02){nu=cutoff[3] }else if(del[3]>=0.02&del[4]<0.02&del[5]<0.02){nu=cutoff[4] }else {nu=cutoff[5]} up_point=min(detected.diff[which(detected.diff[,4]>=nu*par1),4]) # nu=quantile(detected.diff[,4],0.95,na.rm = TRUE) detected.diff[detected.diff[,4]>=up_point,4]=99999 detected.diff[detected.diff[,4]