【单细胞中性粒】慢性病毒性肝炎(复现fig1)
今天来复现这篇文章的Figure 1 ，为后续深入中性粒分群打好基础——
先看看文章的图：
BCDE
是基于t-SNE的降维map
sce.all.int = readRDS('2-harmony/sce.all_int.rds')
library(ggplot2)
# 运行t-SNE降维
seurat_tsne <- RunTSNE(sce.all.int, dims = 1:20, do.fast = TRUE)
p_cell=DimPlot(seurat_tsne, reduction = "tsne",raster = F,group.by = "celltype",
label = T,repel = T)
p_Cell<- p_cell+theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank())+labs(title = NULL)
p_Cell
p_id=DimPlot(seurat_tsne, reduction = "tsne",raster = F,group.by = "donor_id",
label = F,repel = T) +labs(title = NULL)
p_ID <- p_id +theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank())+labs(title = NULL)
p_ID
p_group=DimPlot(seurat_tsne, reduction = "tsne",raster = F,group.by = "treatment",
label = F,repel = T)
p_Group <- p_group +theme(legend.position = c(0.02,0.95),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank())+labs(title = NULL)
p_Group
p_tissue=DimPlot(seurat_tsne, reduction = "tsne",raster = F,group.by = "tissue",
label = F,repel = T) +labs(title = NULL)
p_Tissue <- p_tissue+theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank())
p_Tissue
G
这里我画的确实不咋美观呢【此外，似乎原文的marker不是按照top基因来选的？】
all_markers <- FindAllMarkers(object = sce.all.int)
View(all_markers[1:10,])
save(all_markers,file = "all_markers.RData")
top5 <- all_markers %>% group_by(cluster) %>% top_n(5, avg_log2FC)
width <- 8+0.2*length(unique(Idents(sce.all.int)));width
height <- 8+0.1*length(top5);height
pG <- DotPlot(sce.all.int, features = unique(top5$gene) ,
assay='RNA' ) +
coord_flip() + #翻转
theme(panel.grid = element_blank(),
axis.text.x=element_text(angle = 45, hjust = 1,size = 10),
axis.text.y = element_text(size = 6))+ #轴标签
labs(x=NULL,y=NULL) +
guides(size = guide_legend("Percent Expressed") )+ #legend
scale_color_gradientn(colours = c("#F8E7E2", "#8B0000")) #颜色
pG
F
cg=c('PTPRC','ITGAM','CD68','CD14','FCGR3A','FCGR3B','CD1C','CLEC9A','IL5RA','IL3RA',"CD163","CD5L","ISG15","IFITM1","HLA-DRB1")
cg
# 创建 FeaturePlot
p1 <- FeaturePlot(seurat_tsne,
reduction = "tsne",
features = unique(str_to_upper(cg)),
cols = c("#F5E0B1", "#1A690D"),
ncol = 5)
似乎不如人家的美观，改进一下呢
# 将p1转换为列表以便逐个应用主题
# 循环遍历p1中的每个子图，并将经过主题设置的子图添加到plots_list中
for (i in 1:length(p1)) {
p <- p1[[i]]
plots_list[[i]] <- p + theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank())
}
# 使用wrap_plots将所有子图组合在一起
combined_plot <- wrap_plots(plots_list)
# 显示图形
print(combined_plot)
后期再AI一下应该就差不多了。
H
作者接下来的大部分篇章都是着眼于ISG来的，遗憾的是并未提到ISG基因集的来源（也有可能是我粗心的问题），总之我去2019年一篇nature 子刊上找到了一个ISG基因集合➡A protein-interaction network of interferon-stimulated genes extends the innate immune system landscape。
library(xlsx)
dat <- read.xlsx("./StudyFiles/ISG_from30833792.xlsx",sheetIndex = 3)
ISG <- list(dat$NA.)
# AddModuleScore 打分
sce2 <- AddModuleScore(sce.all.int,
features = ISG,
name = "ISG")
head(sce2@meta.data)
#这里就得到了基因集评分结果，注意列的位置
colnames(sce2@meta.data)[25] <- 'ISG score'
接下来画分组小提琴图并添加p值，参考了ggplot2绘制分组小提琴图并添加统计学显著性标识_cibersort小提琴图-CSDN博客
library(ggrastr)
library(ggpubr)
library(tidyverse)
library(tidyr)
data <- sce2@meta.data %>% .[,c(7,10,11,25)]
rownames(data) <- NULL
colnames(data)
#这一部分筛选出每个细胞类型中最大的值，为添加P值定位而准备的
location <- data %>% group_by(celltype) %>% slice_max(`ISG score`)
location$x <- seq(1,10,by=1)
p <- ggplot(data,aes(celltype,`ISG score`,fill=treatment))+
geom_violin(scale = "width",alpha=0.8,width=0.5,size=0.8)+ #画小提琴图
scale_fill_manual(values = c("#8AB7E9","#FF3362"))+ #分组添加颜色
stat_compare_means(aes(group=treatment), #按分组进行统计检验
method = "wilcox.test",
paired = F, #配对t检验
symnum.args = list(cutpoint=c(0,0.001,0.01,0.05,1),
symbols=c("***","**","*","ns")),
label = "p.signif",
label.y = location$`ISG score`+0.02, #添加显著性符号的位置
size=4.5)+ #显著性符号的大小
geom_segment(data=location, #在显著性符号下面添加一条短线
aes(x=x,y=`ISG score`,
xend=x+0.3,yend=`ISG score`),
size=1)+
xlab("Cell type")+ #X标签
ylab("ISG score")+ #y标签
theme_bw()+
theme(panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.border=element_rect(size=0.5), #给边框加粗
axis.text.x = element_text(angle=60,size=10,vjust = 1,hjust =1,color = "black"),
axis.text.y = element_text(size =10),
legend.position = c(0.90,0.25) )+ ##图例位置
coord_flip() ##坐标轴翻转
p
感觉还缺点啥！
细看的话，会发现原文中是箱线图叠加了小提琴图，
那就再叠加一个箱线图——
p <- ggplot(data,aes(celltype,`ISG score`,fill=treatment)) +
geom_violin(position = position_dodge(width = 0.9), alpha = 0.5) +
geom_boxplot(position = position_dodge(width = 0.9), width = 0.2, outlier.shape = NA) +
theme_minimal()+
scale_fill_manual(values = c("#8AB7E9","#FF3362"))+ #分组添加颜色
stat_compare_means(aes(group=treatment), #按分组进行统计检验
method = "wilcox.test",
paired = F, #配对t检验
symnum.args = list(cutpoint=c(0,0.001,0.01,0.05,1),
symbols=c("***","**","*","ns")),
label = "p.signif",
label.y = location$`ISG score`+0.02, #添加显著性符号的位置
size=4.5)+ #显著性符号的大小
geom_segment(data = location, # 在显著性符号下面添加一条短线
aes(x = x, y = `ISG score`,
xend = x + 0.3, yend = `ISG score`, # 增加短线的长度
group = interaction(celltype, treatment)), # 确保短线对应分组
size = 1)+
xlab("Cell type")+ #X标签
ylab("ISG score")+ #y标签
theme_bw()+
theme(panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.border=element_rect(size=0.5), #给边框加粗
axis.text.x = element_text(angle=60,size=10,vjust = 1,hjust =1,color = "black"),
axis.text.y = element_text(size =10),
legend.position = c(0.90,0.25) )+ ##图例位置
coord_flip() ##坐标轴翻转
print(p)
合理怀疑作者的图是分开画完以后再组装的，R原生图似乎做不到？就看大家能不能忍，不能的话就再加工一下~
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2024-06-22，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录