library(vegan) library("metagenomeSeq") library(biom) ##make function to read OTU matrix from vegan into metagenomeSeq meta_internal <- function (object) { obj <- list(counts = as.data.frame(object), taxa = as.data.frame(rownames(object))) return(obj) } ##read OTU table OTUFinal90<-load_meta(file="LP_90_final.txt") OTUFinal95<-load_meta(file="LP_95_final.txt") OTUFinal97<-load_meta(file="LP_97_final.txt") OTUFinal99<-load_meta(file="LP_99_final.txt") ##read mapping file phenotypeData_A <- load_phenoData(file = "MappingFile_LP.txt") ord90 = match(colnames(OTUFinal90$counts),rownames(phenotypeData_A)) ord95 = match(colnames(OTUFinal95$counts),rownames(phenotypeData_A)) ord97 = match(colnames(OTUFinal97$counts),rownames(phenotypeData_A)) ord99 = match(colnames(OTUFinal99$counts),rownames(phenotypeData_A)) ##pull out the matched rows from mapping file OTU_pheno_match90 = AnnotatedDataFrame(phenotypeData_A[ord90,]) OTU_pheno_match95 = AnnotatedDataFrame(phenotypeData_A[ord95,]) OTU_pheno_match97 = AnnotatedDataFrame(phenotypeData_A[ord97,]) OTU_pheno_match99 = AnnotatedDataFrame(phenotypeData_A[ord99,]) #combine OTU table with matched mapping file data, file type is MRexperiment ObjData90<-newMRexperiment(OTUFinal90$counts,phenoData = OTU_pheno_match90) ObjData95<-newMRexperiment(OTUFinal95$counts,phenoData = OTU_pheno_match95) ObjData97<-newMRexperiment(OTUFinal97$counts,phenoData = OTU_pheno_match97) ObjData99<-newMRexperiment(OTUFinal99$counts,phenoData = OTU_pheno_match99) #set sequencing depths to test, set list of OTU tables to run for loops seq.depth = c(100, 1000, 5000, 10000, 20000) sample.size = c(5,10,20,60) #Sampling_combined <- rep(NA, 3*length(seq.depth)) OTU.tables = list(ObjData90, ObjData95, ObjData97, ObjData99) sims <- 1000 #Number of bootstrap simulation to run #set up results that you want mBC.size <- rep(NA, length(seq.depth)) lowBC.size <- rep(NA, length(seq.depth)) highBC.size <- rep(NA, length(seq.depth)) mJac.size <- rep(NA, length(seq.depth)) lowJac.size <- rep(NA, length(seq.depth)) highJac.size <- rep(NA, length(seq.depth)) mBCPerm.size <- rep(NA, length(seq.depth)) lowBCPerm.size <- rep(NA, length(seq.depth)) highBCPerm.size <- rep(NA, length(seq.depth)) mJacPerm.size <- rep(NA, length(seq.depth)) lowJacPerm.size <- rep(NA, length(seq.depth)) highJacPerm.size <- rep(NA, length(seq.depth)) R_BC <- rep(NA,sims) R_Jac <- rep(NA,sims) PermR_BC <- rep(NA,sims) PermR_Jac <- rep(NA,sims) set.seed(5588) #Random seed for replication for (k in 1:4){ #Make list of tables Data <- OTU.tables[[k]] Tips = which(pData(Data)$Description =="Tip") obj_tips = Data[,Tips] Bases = which(pData(Data)$Description =="Base") obj_bases = Data[,Bases] #Pull out the OTU table, these are matrices LPData_mat <- MRcounts(Data, norm = FALSE, log = FALSE) #transform the matrices so that rrarefy can read it properly tLPData_mat <-t(LPData_mat) for (s in 1:length(sample.size)){ for (j in 1:length(seq.depth)){ for (i in 1:sims){ rOTU_table <- rrarefy(tLPData_mat, seq.depth[j]) rarefied_Data<-meta_internal(t(rOTU_table)) ord = match(colnames(rarefied_Data$counts),rownames(phenotypeData_A)) OTU_pheno_match = AnnotatedDataFrame(phenotypeData_A[ord,]) RareObjData<-newMRexperiment(rarefied_Data$counts,phenoData = OTU_pheno_match) sub_tips <- obj_tips[,sample(1:63,size=sample.size[s], replace = T)] sub_bases <- obj_bases[,sample(1:64,size=sample.size[s], replace = T)] subsample_names <- c(colnames(sub_tips),colnames(sub_bases)) subsample_names <-substring(subsample_names,1,4) raresubdata<-RareObjData[,subsample_names] mat = MRcounts(raresubdata, norm = FALSE, log = TRUE) tmat <- t(mat) BCData <- vegdist(tmat, method="bray", binary=FALSE, diag=TRUE, upper=TRUE) JacData <- vegdist(tmat, method="jaccard", binary=TRUE, diag=TRUE, upper=TRUE) Perm_BC<-adonis(BCData ~raresubdata$Description, permutations = 10, distance = "bray") Perm_Jac<-adonis(JacData~ raresubdata$Description, permutations = 10, distance = "jaccard") ANOSIMRData_Jac <- anosim(JacData, raresubdata$Description, permutations = 1, distance = "jaccard", strata = NULL, parallel = getOption("mc.cores")) ANOSIMRData_BC <- anosim(BCData, raresubdata$Description, permutations = 1, distance = "bray", strata = NULL, parallel = getOption("mc.cores")) # ANOSIMRData$statistic R_BC[i] <- ANOSIMRData_BC$statistic R_Jac[i] <- ANOSIMRData_Jac$statistic # PermANOVARData$statistic PermR_BC[i] <- Perm_BC$aov.tab$R2[1] PermR_Jac[i] <- Perm_Jac$aov.tab$R2[1] } #close i loop mBC.size[j] <- mean(R_BC) lowBC.size[j] <- sort(R_BC)[25] highBC.size[j] <- sort(R_BC)[975] mJac.size[j] <- mean(R_Jac) lowJac.size[j] <- sort(R_Jac)[25] highJac.size[j] <- sort(R_Jac)[975] mBCPerm.size[j] <- mean(PermR_BC) lowBCPerm.size[j] <- sort(PermR_BC)[25] highBCPerm.size[j] <- sort(PermR_BC)[975] mJacPerm.size[j] <- mean(PermR_Jac) lowJacPerm.size[j] <- sort(PermR_Jac)[25] highJacPerm.size[j] <- sort(PermR_Jac)[975] } #close j loop TB_ANOSIM_BC_seq <- c(mBC.size, lowBC.size, highBC.size) TB_ANOSIM_Jac_seq <- c(mJac.size, lowJac.size, highJac.size) TB_Perm_BC_seq <- c(mBCPerm.size, lowBCPerm.size, highBCPerm.size) TB_Perm_Jac_seq <- c(mJacPerm.size, lowJacPerm.size, highJacPerm.size) colnameBC <- paste("ANOSIM_R_BC_",sample.size[s],"res",k,sep="") colnameJac <- paste("ANOSIM_R_Jac_",sample.size[s],"res",k,sep="") colnamePermBC <- paste("Perm_R_BC_",sample.size[s],"res",k,sep="") colnamePermJac <- paste("Perm_R_Jac_",sample.size[s],"res",k,sep="") if ((s+k) < 3){ Sampling_combined<-cbind(TB_ANOSIM_BC_seq,TB_ANOSIM_Jac_seq,TB_Perm_BC_seq,TB_Perm_Jac_seq) colnames(Sampling_combined)<-c(colnameBC,colnameJac,colnamePermBC,colnamePermJac) } else { Sampling_combined_new<-cbind(TB_ANOSIM_BC_seq,TB_ANOSIM_Jac_seq,TB_Perm_BC_seq,TB_Perm_Jac_seq) colnames(Sampling_combined_new)<-c(colnameBC,colnameJac,colnamePermBC,colnamePermJac) Sampling_combined <-cbind(Sampling_combined,Sampling_combined_new) } } #close s loop Sampling_combined_new<-cbind(TB_ANOSIM_BC_seq,TB_ANOSIM_Jac_seq,TB_Perm_BC_seq,TB_Perm_Jac_seq) colnames(Sampling_combined_new)<-c(colnameBC,colnameJac,colnamePermBC,colnamePermJac) Sampling_combined <-cbind(Sampling_combined,Sampling_combined_new) } #closing for k write.csv(Sampling_combined, file = "LP_Sequencing_TaxRes.csv")