I am trying to perform differential gene expression using a t-test on spatial RNA-sequencing data. There are a couple of different annotations/groups indicating different structures (ANN2 in code): AML area, Taggregate, immatureTLS, matureTLS, and microcluster. ANN1 relates to one of the 3 different patients.

The error I get: Error in h(simpleError(msg, call)) : error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': not enough 'y' observations I don’t understand how my data has not enough y-observations, and how I could overcome this error. I have searched google and other blogs, but I wasn’t able to resolve it. The code I use (all the code I used before is shown on this website: https://bioconductor.org/packages/devel/workflows/vignettes/GeoMxWorkflows/inst/doc/GeomxTools_RNA-NGS_Analysis.html )

7.1 Differential Expression

plots<-list()
tables<-list()
labels<-list()
test<-"ttest"
mtc<-"BY"
#options: "holm"       "hochberg"   "hommel"     "bonferroni" "BH"         "BY"         "fdr" 
counter=1

comps_df<-data.frame(comp='',val='')

for (active_group1 in unique(ann$segment)) {
    for (active_group2 in unique(ann$segment)) {
     
      #supress reduncant compares
      if(active_group1==active_group2) {next}
      comp<-paste(sort(c(active_group1,active_group2)),collapse = "_")
      if(comp %in% comps_df$comp) {next}
      temp_df<-data.frame(comp=comp ,val=1)
      comps_df<-rbind(comps_df,temp_df)
      
      labels[[counter]]<-paste(active_group1," vs ", active_group2)
      group1<-log_q[,rownames(ann)[ann$segment==active_group1]]
      group2<-log_q[,rownames(ann)[ann$segment==active_group2]]
      
      #run t_tests  
      results<-as.data.frame ( apply(log_q, 1, function(x) t.test(x[colnames(group1)],x[colnames(group2)])$p.value) )
      colnames(results)<-"raw_p_value"
      
      #multiple_testing_correction
      adj_p_value<- p.adjust(results$raw_p_value,method=mtc)
      results<-cbind(results,adj_p_value)
      
      #calc_fdr
      FDR<- p.adjust(results$raw_p_value,method="fdr")
      results<-cbind(results,FDR)
     
      #fold_changes
      #as base data is already log transformed, means need to be subtracted to get FC in log space
      fchanges<-as.data.frame( apply(log_q, 1, function(x) (mean(x[colnames(group1)]) - mean(x[colnames(group2)]) ) ) )
      colnames(fchanges)<-"FC"
      #paste("FC",active_group1," / ",active_group2)
      results<-cbind(results,fchanges)
      
      #add genenames
      results$Gene<-rownames(results)
      
      #set categories based on P-value & FDR for plotting
      results$Color <- "NS or FC < 0.5"
      results$Color[results$adj_p_value < 0.05] <- "P < 0.05"
      results$Color[results$FDR < 0.05] <- "FDR < 0.05"
      results$Color[results$FDR < 0.001] <- "FDR < 0.001"
      results$Color[abs(results$FC) < 1] <- "NS or FC < 1"
      results$Color <- factor(results$Color,
                              levels = c("NS or FC < 1", "P < 0.05", "FDR < 0.05", "FDR < 0.001"))
      
      #vulcanoplot
      
      # pick top genes for either side of volcano to label
      # order genes for convenience:
      
      results$invert_P <- (-log10(results$adj_p_value)) * sign(results$FC)
      top_g <- c()
      top_g <- c(top_g,
                 results[ind, 'Gene'][
                   order(results[ind, 'invert_P'], decreasing = TRUE)[1:15]],
                 results[ind, 'Gene'][order(results[ind, 'invert_P'], decreasing = FALSE)[1:15]])
      top_g<- unique(top_g)
      results <- results[, -1*ncol(results)] # remove invert_P from matrix
      
      # Graph results
      plots[[counter]]<- ggplot(results,
                                      aes(x = FC, y = -log10(adj_p_value),
                                          color = Color, label = Gene)) +
        geom_vline(xintercept = c(1, -1), lty = "dashed") +
        geom_hline(yintercept = -log10(0.05), lty = "dashed") +
        geom_point() +
        labs(x = paste("Enriched in", active_group2," <- log2(FC) -> Enriched in", active_group1),
             y = "Significance, -log10(P)",
             color = "Significance") +
        scale_color_manual(values = c(`FDR < 0.001` = "dodgerblue",
                                      `FDR < 0.05` = "lightblue",
                                      `P < 0.05` = "orange2",
                                      `NS or FC < 0.5` = "gray"),
                           guide = guide_legend(override.aes = list(size = 4))) +
        scale_y_continuous(expand = expansion(mult = c(0,0.05))) +
        geom_text_repel(data = subset(results, FDR<0.001 & (-1>FC| FC>1)),
                        point.padding = 0.15, color = "black", size=3.5,
                        min.segment.length = .1, box.padding = .2, lwd = 2,
                        max.overlaps = 50) +
        theme_bw(base_size = 20) +
        theme(legend.position = "bottom") +
        ggtitle(paste(test, mtc,"multitest corr"))
      
      #store tables for display later
      tables[[counter]]<-results
      
       counter = counter+1
      #datatable(subset(results, Gene %in% GOI), rownames=FALSE,caption = paste("DE results ", active_group1," vs ", active_group2))
    }
  }

grid.arrange(grobs=plots,ncol=2)

#strangly does not appear in html output??

for (c in (2:counter-1)) {
  #Gene %in% GOI
  
  print(datatable( subset(tables[[c]], Color == "FDR < 0.001" ), 
           rownames=FALSE,
           extensions = 'Buttons', options = list (
              dom = 'Bftrip',
              buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
            ),
           caption = paste("DE results ", labels[[1]]),filter='top') %>% formatRound(columns=c("raw_p_value","adj_p_value","FDR","FC"), digits=3))
  cat('\n\n<!-- -->\n\n')
}            


[normalised data example][1]Data type which is used as input: 

 


  [1]: https://i.stack.imgur.com/Yt0DJ.png


Any help would be greatly appreciated! Thanks