【单细胞中性粒】慢性病毒性肝炎(复现fig1)

生信菜鸟团

发布于 2024-06-11 17:34:33

10200

代码可运行

文章被收录于专栏：生信菜鸟团

运行总次数：0

代码可运行

Single-cell atlas of the liver myeloid compartment before and after cure of chronic viral hepatitis - Journal of Hepatology

今天来复现这篇文章的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)