Seurat单细胞基因显著性检验函数及批量添加显著性

❝

3、付费文章集合有打包价哦！

详情请联系作者：

❞

这篇帖子的前身可以追溯到：玩转单细胞(2):Seurat批量做图修饰。当时有人问了一个问题，可以添加显著性吗？我们的回答是你需要提取相关组的表达量，进行检验后再用ggplot函数添加即可；或者直接提取数据用ggplot作图那么显著添加也就不成问题了。时隔3月，我们这里提供 了一种函数，可以进行基因在两组之间的显著性分析。同时可进行批量的基因分析。并输出dataframe结果。同时直接在Vlnplot下循环添加显著性。但缺点是只能进行两组比较分析。完整代码已上传群文件！

一般的seurat小提琴图绘制：

library(Seurat)library(ggplot2)library(ggpubr)library(dplyr)
VlnPlot(mouse_data, features = 'S100a8', group.by = 'orig.ident')+  theme_classic() +   theme(axis.text.x = element_text(size = 10,color="black"),        axis.text.y = element_text(size = 10,color="black"),        axis.title.y= element_text(size=12,color="black"),        axis.title.x = element_blank(),        legend.position='none')

显著性检验函数，有点长，可自行保存成R文件，然后每次使用的时候source一下就可以了。

singlecell_gene_test <- function(SerautObj,                            genes.use,                            group.by=NULL,                            assay = "RNA",                            comp = NULL,                            alpha_start = .05,                            Bonferroni = T,                           only_postive =F) {  p_val.out <- c()  stat.out <- c()  condition.out <- c()  gene.out <- c()  if (only_postive == F){    for (gene in genes.use){      group1_cellname = rownames([e-mail][[e-mail][[group.by]] == comp[1],])      group1_exp = SerautObj@assays[[assay]]@data[gene, group1_cellname] 
      group2_cellname = rownames([e-mail][[e-mail][[group.by]] == comp[2],])      group2_exp = SerautObj@assays[[assay]]@data[gene, group2_cellname]      t_out = t.test(group1_exp, group2_exp)      cond = paste(comp[1], comp[2], sep = "_")      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)    }  }    else{      for (gene in genes.use){        group1_cellname = rownames([e-mail][[e-mail][[group.by]] == comp[1],])        group1_exp = SerautObj@assays[[assay]]@data[gene, group1_cellname]        group1_exp <- group1_exp[which(group1_exp>0)] 

        group2_cellname = rownames([e-mail][[e-mail][[group.by]] == comp[2],])        group2_exp = SerautObj@assays[[assay]]@data[gene, group2_cellname]        group2_exp <- group2_exp[which(group2_exp>0)] 
        t_out = t.test(group1_exp, group2_exp)        cond = paste(comp[1], comp[2], sep = "_")        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)      }
  }
  if (Bonferroni == T){    new_alpha = alpha_start/(2*length(genes.use))    cat(paste("\n", "P-value for significance: p <", new_alpha, "\n"))    sig_out = p_val.out < new_alpha    dfOUT<- data.frame(gene=gene.out, condition = condition.out, p_val = p_val.out, statistic = stat.out, significant = sig_out)
    dfOUT$sig = ifelse(dfOUT$p_val > 0.05, "ns",                       ifelse(dfOUT$p_val > 0.01, '*',                              ifelse(dfOUT$p_val > 0.001, "**", "****")))
    }
  else{    dfOUT<- data.frame(gene=gene.out, condition = condition.out, p_val = p_val.out, statistic = stat.out)    dfOUT$sig = ifelse(dfOUT$p_val > 0.05, "ns",                       ifelse(dfOUT$p_val > 0.01, '*',                              ifelse(dfOUT$p_val > 0.001, "**", "****")))  }
  return(dfOUT)}

显著性检验：

A <- singlecell_gene_test(mouse_data,                     genes.use = c('S100a8','Ltf','Ncf1','Ly6g','Anxa1','Il1b'),                    group.by = 'orig.ident',                     comp = c("10X_ntph_F", "10X_ntph_M"))
A1 <- singlecell_gene_test(mouse_data,                          genes.use = c('S100a8','Ltf','Ncf1','Ly6g','Anxa1','Il1b'),                          group.by = 'orig.ident',                           comp = c("10X_ntph_F", "10X_ntph_M"),                          only_postive = T)

作图即可：

anno_pvalue <- format(A$p_val, scientific = T,digits = 3) anno_sig <- A$sig
plots_violins <- VlnPlot(mouse_data,                          cols = c("limegreen", "navy"),                         pt.size = 0,                         group.by = "orig.ident",                         features = c('S100a8','Ltf','Ncf1','Ly6g','Anxa1','Il1b'),                          ncol = 3,                          log = FALSE,                         combine = FALSE)
for(i in 1:length(plots_violins)) {  data <- plots_violins[[i]]$data  colnames(data)[1] <- 'gene'  plots_violins[[i]] <- plots_violins[[i]] +     theme_classic() +     theme(axis.text.x = element_text(size = 10,color="black"),          axis.text.y = element_text(size = 10,color="black"),          axis.title.y= element_text(size=12,color="black"),          axis.title.x = element_blank(),          legend.position='none')+    scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))+    scale_x_discrete(labels = c("Female","Male"))+    geom_signif(annotations = anno_sig[i],                y_position = max(data$gene)+0.5,                xmin = 1,                xmax = 2,                tip_length = 0)}
CombinePlots(plots_violins)

或者添加p值：

plots_violins <- VlnPlot(mouse_data,                          cols = c("limegreen", "navy"),                         pt.size = 0,                         group.by = "orig.ident",                         features = c('S100a8','Ltf','Ncf1','Ly6g','Anxa1','Il1b'),                          ncol = 3,                          log = FALSE,                         combine = FALSE)for(i in 1:length(plots_violins)) {  data <- plots_violins[[i]]$data  colnames(data)[1] <- 'gene'  plots_violins[[i]] <- plots_violins[[i]] +     theme_classic() +     theme(axis.text.x = element_text(size = 10,color="black"),          axis.text.y = element_text(size = 10,color="black"),          axis.title.y= element_text(size=12,color="black"),          axis.title.x = element_blank(),          legend.position='none')+    scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))+    scale_x_discrete(labels = c("Female","Male"))+    geom_signif(annotations = anno_sig[i],                y_position = max(data$gene)+0.5,                xmin = 1,                xmax = 2,                tip_length = 0)}
CombinePlots(plots_violins)

好了。这就是所有内容了，其实这样检验你用不用得到倒是其次，主要是这里面包含一些小的细节知识点，学会了就能和其他内容融汇贯通了，自己感悟吧！