looping not working in R

Question

I am trying to get a function to solve a small problem. I have to two list, and each list comprise n samples. Each sample has a variable amount of identifiers of bacteria (in the example letters, in my problem bacterial identifiers like OTU1-OTUn, in both cases are “character”). One list comprise samples from diet, and the another list samples from gut contents. I want to know for each sample of list gut, how many bacteria from diet are in the gut and how many bacteria in the gut do not come from diet. This was easily done when working with phyloseq object were diet and gut are both phyloseq objects with n samples each.

Bacteria_from_diet<-length(intersect(taxa_names(gut),taxa_names(diet))
Bacteria_not_diet<-length(taxa_names(diet)- Bacteria_from_diet

However, this “summarizes” the result over the n samples of gut and diet, I mean, like if I collapse data by sample, and I need some measure of variation. I have tried the following code in R:

diet<-list(DL1=c("A","B","C"),DL2=c("A","C","D"),DL3=c("B","D","E"),DL4=c("B","D","E"))
gut<-list(DL5=c("A","F","G"),DL6=c("B","F","H"),DL7=c("D","H","J"),DL8=c("A","G","F"))


gut_vs_diet <- function(a,b) ## a is diet and b is gut
{
xx<-10 
gut = numeric(xx) 
diet = numeric(xx)
all<-unlist(lapply(b,length))  ### get the number of elements of each element of list b 
for(i in seq_along(b)){ #### loop over b (gut) to get:
diet<-length(intersect(b[[i]],a[[i]])) ### the number of elements of diet are present in gut 
gut = all-diet ## the number of elements of gut that not come from diet
}
gutvsdiet = data.frame(all,gut,diet)
return(gutvsdiet)
} 

When running the funtion I obtain this result, which is not correct

gut_vs_diet(diet,gut)
    all gut diet
DL5   3   3    0
DL6   3   3    0
DL7   3   3    0
DL8   3   3    0

In some cases, I was able to get some value in diet column, but the function randomly choose the diet sample.

I do not know where could be the mistake.Anyway, I would like to do this iteratively, I mean, get the values for each sample of gut compared with all samples of diet. Alternatively, I can run a replicate(10,gut_vs_diet(sample(diet),sample(gut)) to get random comparisons and avoid somekind of bias.

Thank you very much for your help

Manuel

Answer 1

Here is my version of your code:

diet <- list(DL1=c("A","B","C"), DL2=c("A","C","D"), DL3=c("B","D","E"), DL4=c("B","D","E"))
gut  <- list(DL5=c("A","F","G"), DL6=c("B","F","H"), DL7=c("D","H","J"), DL8=c("A","G","F"))

gut_vs_diet <- function(a, b) ## a is diet and b is gut
{
  all  <- lengths(b)  ### get the number of elements of each element of list b 
  diet <- mapply(function(ai, bi) length(intersect(ai, bi)), a, b)
  # diet <- lengths(mapply(intersect, a, b))  ## a variant
  data.frame(all, gut=all-diet, diet)
} 
gut_vs_diet(diet,gut)
# > gut_vs_diet(diet,gut)
#     all gut diet
# DL5   3   2    1
# DL6   3   3    0
# DL7   3   2    1
# DL8   3   3    0

Answer 2

Just for completion, with a for loop it would look like this. Note that you need to subtract all[[i]] - diet and construct the dataframe inside the loop, otherwise you will just fill it with the last result of the loop, which is, data.frame(all = c(3,3,3,3), gut = 3, diet = 0)

diet <- list(DL1 = c("A", "B", "C"), DL2 = c("A", "C", "D"), DL3 = c("B", "D", "E"), DL4 = c("B", "D", "E"))
gut <- list(DL5 = c("A", "F", "G"), DL6 = c("B", "F", "H"), DL7 = c("D", "H", "J"), DL8 = c("A", "G", "F"))


gut_vs_diet <- function(a, b)
   {
   all <- lengths(b)
   gutvsdiet <- NULL
   for (i in seq_along(b)) {
      diet <- length(intersect(b[[i]], a[[i]]))
      gut <- all[[i]] - diet
      resultForThisListElement <- c(all[[i]], gut, diet)
      gutvsdiet <- rbind(gutvsdiet, resultForThisListElement)
   }
   colnames(gutvsdiet) <- c("all", "gut", "diet")
   return(gutvsdiet)
}

gut_vs_diet(diet, gut)

Answer 3

As @jogo suggested in a comment, you can use mapply instead of your for-loop:

FOO <- function(x, y){
  all <- lengths(y)
  diet <- mapply(function(a, b){
    length(intersect(b, a))
  }, x, y)
  gut <- all - diet
  return(data.frame(all, gut, diet))
}

> FOO(diet, gut)
    all gut diet
DL5   3   2    1
DL6   3   3    0
DL7   3   2    1
DL8   3   3    0