這篇帖子的前身可以追溯到：玩轉(zhuǎn)單細(xì)胞(2):Seurat批量做圖修飾。當(dāng)時(shí)有人問(wèn)了一個(gè)問(wèn)題，可以添加顯著性嗎眷细？我們的回答是你需要提取相關(guān)組的表達(dá)量鸠按，進(jìn)行檢驗(yàn)后再用ggplot函數(shù)添加即可醒陆；或者直接提取數(shù)據(jù)用ggplot作圖那么顯著添加也就不成問(wèn)題了。時(shí)隔3月葫掉，我們這里提供 了一種函數(shù)些举，可以進(jìn)行基因在兩組之間的顯著性分析。同時(shí)可進(jìn)行批量的基因分析俭厚。并輸出dataframe結(jié)果户魏。同時(shí)直接在Vlnplot下循環(huán)添加顯著性。但缺點(diǎn)是只能進(jìn)行兩組比較分析挪挤。完整代碼已上傳群文件叼丑！ 一般的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')
顯著性檢驗(yàn)函數(shù)，有點(diǎn)長(zhǎng)扛门，可自行保存成R文件鸠信，然后每次使用的時(shí)候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(SerautObj@meta.data[SerautObj@meta.data[[group.by]] == comp[1],])
group1_exp = SerautObj@assays[[assay]]@data[gene, group1_cellname]
group2_cellname = rownames(SerautObj@meta.data[SerautObj@meta.data[[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(SerautObj@meta.data[SerautObj@meta.data[[group.by]] == comp[1],])
group1_exp = SerautObj@assays[[assay]]@data[gene, group1_cellname]
group1_exp <- group1_exp[which(group1_exp>0)]
group2_cellname = rownames(SerautObj@meta.data[SerautObj@meta.data[[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)
}
顯著性檢驗(yàn)：
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)
好了论寨。這就是所有內(nèi)容了星立，其實(shí)這樣檢驗(yàn)?zāi)阌貌挥玫玫降故瞧浯危饕沁@里面包含一些小的細(xì)節(jié)知識(shí)點(diǎn)葬凳，學(xué)會(huì)了就能和其他內(nèi)容融匯貫通了贞铣，自己感悟吧！更多精彩內(nèi)容請(qǐng)至KS科研分享與服務(wù)公眾號(hào)
