#' Calculate transfer statistics including permutation testing for two phyloseq objects
#'
#'@param phy1 phyloseq object for the giving environment
#'@param phy2 phyloseq object for the recieving environment
#'@param identifier in sample_data for mathcing dyads(character)
#'@import dplyr tidyverse broom
#'@export
randpermutationTransferStats <- function(phy1, phy2, identifier, nperm = 1000, verbose = F){
# dig out the count_tables
TAXtb <- data.frame(tax_table(phy1))
o1 <- phy1 %>% otu_table()
o2 <- phy2 %>% otu_table()
if (taxa_are_rows(phy1)){
o1 <- o1 %>% t()
o2 <- o2 %>% t()
}
# find intersect of the two datasets
ID1 <- sample_data(phy1)[[identifier]]
ID2 <- sample_data(phy2)[[identifier]]
# order according to ID1 and ID2
id1 <- order(ID1)
id2 <- order(ID2)
ID1 <- ID1[id1]
ID2 <- ID2[id2]
o1 <- o1[id1,]
o2 <- o2[id2,]
# find intersect
IDmatch <- intersect(ID1,ID2)
# cut out samples
o1 <- o1[ID1 %in% IDmatch,]
o2 <- o2[ID2 %in% IDmatch,]
ID1 <- ID1[ID1 %in% IDmatch]
ID2 <- ID2[ID2 %in% IDmatch]
deptho1 <- apply(o1,1,sum)
deptho2 <- apply(o2,1,sum)
totAbuM <- sum(deptho1)
totAbuC <- sum(deptho2)
# remove non-testabel OTU'
s01 <- apply(o1==0,2,sum)
s02 <- apply(o2==0,2,sum)
nn <- dim(o2)[1]
icvar <- s01>0 & s02>0 & s01<nn & s02<nn
o1 <- o1[,icvar]
o2 <- o2[,icvar]
p <- sum(icvar)
print(dim(o1))
TAXtb <- TAXtb[icvar,]
abuM <- apply(o1,2,sum)
abuC <- apply(o2,2,sum)
abuMrel <- o1 %>% sweep(1,deptho1,'/') %>% apply(2,mean)
abuCrel <- o2 %>% sweep(1,deptho2,'/') %>% apply(2,mean)
STAT <- get2by2table(o1,o2)
STAT <- data.frame(STAT,abuM,abuMrel,abuC,abuCrel)
STAT$otu <- rownames(STAT)
# make permutation testing
if (nperm>0){
permSTAT <- array(dim = c(dim(STAT)[1],4,nperm))
for (i in 1:nperm){
stp <- get2by2table(o1[sample(dim(o1)[1]),],o2)
permSTAT[,,i] <- as.matrix(stp)
}
dimnames(permSTAT)[[1]] <- rownames(stp)
dimnames(permSTAT)[[2]] <- colnames(stp)
permSTATfisher <- array(dim = c(dim(STAT)[1],4,nperm))
for (i in 1:nperm){
if (verbose) {print(i)}
aa <- permSTAT[,,i] %>% data.frame()
for (j in 1:dim(permSTAT)[1]){
stp2 <- getFisher(aa[j,],doglm = F)
# permSTATfisher[j,,i] <- c(stp2$Fisher_estimate, stp2$Fisher_p.value)
permSTATfisher[j,,i] <- c(stp2$Fisher_estimate,stp2$or_biascorr, stp2$Gtest_p.value, stp2$Fisher_p.value)
}
}
dimnames(permSTATfisher)[[1]] <- rownames(stp)
# dimnames(permSTATfisher)[[2]] <- c('Fisher_estimate', 'Fisher_p.value')
dimnames(permSTATfisher)[[2]] <- c('Fisher_estimate','or_biascorr', 'Gtest_p.value', 'Fisher_p.value')
}
else {
permSTAT <- NULL
permSTATfisher = NULL
}
STAT <- STAT %>%
group_by(otu) %>%
do(getFisher(x = .)) %>%
ungroup %>%
left_join(STAT,by = 'otu')
STAT <- merge(STAT,TAXtb,by.x = 'otu',by.y = 'row.names')
STAT$totAbuM <- totAbuM
STAT$totAbuC <- totAbuC
STAT <- STAT %>%
mutate(Fisher_estimatetr = truncateZerosInf(Fisher_estimate),
Glm_ortr = truncateZerosInf(Glm_or))
return(list(STAT,permSTAT, permSTATfisher))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.