Search code examples
rbioinformaticsseurat

Tri-nested loop in R

Currently, I am writing a function that uses a tri nested loop in R; however, it appears to have some strange behavior. I am noticing the following issues:

1) gene does not appear to be appended gene to gene.out at the end of the loop, the list I get from gene.out is gene.use.

2) cond ONLY makes duplicates (i.e 22_22, 35_35, ect)

As far as I can tell neither of these things should be happening in the third loop. Is this some weird R loop behavior or is this a coding mistake?

Here is the code in question:

for (gene in genes.use){

    for (i in groups){
      cat(paste("i: ",i, "\n"))

      i_cells = rownames([email protected][[email protected][[group.by]] == i,])
      i_vector = SerautObj@assays[[assay]]@data[gene, i_cells]

      for(j in groups){
        cat(paste("j: ",j, "\n"))

        j_cells = rownames([email protected][[email protected][[group.by]] == j,])
        j_vector = SerautObj@assays[[assay]]@data[gene, j_cells]


        cond = paste(i, j, sep = "_")
        cat(paste(gene, cond, sep = "\n"))

        #preform t-test
        t_out = t.test(i_vector, j_vector)

        #constuct outs

        condition.out <- c(condition.out, cond)



        stat.out <- c(stat.out, t_out[["statistic"]])
        p_val.out <- c(p_val.out, t_out[["p.value"]])
        gene.out <- c(gene.out, gene)
        }
    }
}

edit:

Forgot to include, When I do print(paste("i: ", i) in the i loop and print(paste("j: ", j)) I get:

i: group1
i: group2
i: group3
j: group1
j: group2
j: group3

Toy set data from https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html:

[email protected]
structure(list(orig.ident = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L), .Label = "pbmc3k", class = "factor"), nCount_RNA = c(0, 
0, 0, 0, 0, 0, 0, 0, 0, 0), nFeature_RNA = c(0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L), group = c(1, 2, 3, 4, 5, 1, 2, 3, 4, 
5)), row.names = c("AAACATACAACCAC", "AAACATTGAGCTAC", "AAACATTGATCAGC", 
"AAACCGTGCTTCCG", "AAACCGTGTATGCG", "AAACGCACTGGTAC", "AAACGCTGACCAGT", 
"AAACGCTGGTTCTT", "AAACGCTGTAGCCA", "AAACGCTGTTTCTG"), class = "data.frame")

SerautObj@assays[[assay]]@data

    new("dgCMatrix", i = integer(0), p = c(0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L), Dim = c(10L, 10L), Dimnames = list(c("AL627309.1", 
"AP006222.2", "RP11-206L10.2", "RP11-206L10.9", "LINC00115", 
"NOC2L", "KLHL17", "PLEKHN1", "RP11-54O7.17", "HES4"), c("AAACATACAACCAC", 
"AAACATTGAGCTAC", "AAACATTGATCAGC", "AAACCGTGCTTCCG", "AAACCGTGTATGCG", 
"AAACGCACTGGTAC", "AAACGCTGACCAGT", "AAACGCTGGTTCTT", "AAACGCTGTAGCCA", 
"AAACGCTGTTTCTG")), x = numeric(0), factors = list())

genes.use = c("PLEKHN1", "HES4", "NOC2L")
groups = Map(c, unique([email protected]$groups))

Thanks for reading!

Solution

  • It turned out that the issue was with:

    Map(c, unique([email protected]$groups))

    I then tried:

    as.list(unique(SerautObj[["time"]]))

    This had the same problem, but this was fixed with:

    unlist(as.list(unique(SerautObj[["time"]])))

    It appears that there is weird for loop behavior over lists and you need to go to an atomic type vector or you can end up with duplicates. I'm guessing there is some strange referencing happening or something with loops over lists.