【原】SARS-CoV-2感染的雪貂支气管肺泡灌洗液单细胞转录组数据挖掘（1）降维聚类分群

健明 2022-02-19

展开全文

给学徒们收集整理了几套带GitHub源代码的文献图表合辑，让优秀者一点一滴拆解开来分享给大家。（全部的代码复制粘贴即可运行，欢迎尝试以及批评指正）

现在是雪貂支气管肺泡灌洗液单细胞转录组显示SARS-CoV-2感染期间巨噬细胞的顺序变化专辑第一讲：主要是降维聚类分群和每个单细胞亚群的生物学命名！

下面是前年实习生（日行一膳）的分享

本次复现的是于2021年7月27日发表在Nature Communications上的”Single-cell transcriptome of bronchoalveolar lavage fluid reveals sequential change of macrophages during SARS-CoV-2 infection in ferrets“

本次的数据集可在GEO中进行下载

1642516361162

背景

在当前的冠状病毒(COVID-19)大流行期间，横向研究迅速拓宽了对严重急性呼吸综合征冠状病毒(SARS-CoV-2)免疫应答的理解。目前，正在探讨治疗方式的机制，包括抗病毒和抗炎药物，并伴随临床试验。然而，由于对人类受试者的观察性研究的固有局限性，很少能获得从初始阶段到SARS-CoV-2感染消除的免疫反应的纵向描述。并且目前大多数可用的免疫细胞转录组分析都来自横断面研究，更重要的是，由于缺乏SARS-CoV-2感染前获得的数据，无法将感染状态与未感染状态进行比较。此外，BAL的侵袭性阻碍了SARS-CoV-2感染期间危重患者序贯标本的采集。这些局限性可以通过分析SARS-CoV-2感染的动物模型来克服。

雪貂(Mustela putorius furo)被广泛用作呼吸道病毒病原学研究的动物模型。自1933年发现雪貂对流感病毒的天然易感性以来，这些动物已被用于重现几种人类呼吸道病毒疾病的病程，包括副流感病毒、呼吸道合胞病毒和SARS-CoV。此外，它们的组织解剖学特征：包括上、下呼吸道长度的比例、气道腺体密度和末端细支气管结构等，为模拟人类呼吸道感染提供了最佳条件。

这里,我们执行的scRNA-seq BAL流体样品,这是用于调查的免疫学变化肺,SARS-CoV-2感染雪貂与阴性控制相比,在感染后2天(dpi)(早期SARS-CoV-2感染高峰病毒效价)和5 dpi与组织病理学(决议阶段)。对雪貂肺免疫微环境的景观分析揭示了这段时间免疫细胞比例和特征的动态变化。具体来说，根据独特的基因表达模式将巨噬细胞群体划分为十个不同的亚群，并描述了它们的转录组的时间变化。

Fig. 1 Single-cell transcriptomes of BAL fluid cells from SARS-CoV-2-infected ferrets.

1642600784950

#安装所需要的R包
install.packages('Seurat')
install.packages("tidyverse")
install.packages("ggpubr")
install.packages("cowplot")
install.packages("dplyr")
install.packages("ggplot2")
install.packages("tidyr")

#载入R包
library(Seurat)
library(tidyverse)
library(ggpubr)
library(cowplot)
library(dplyr)
library(ggplot2)
library(tidyr)

##在GEO：GSE171828中直接下载即可
MPcov <- readRDS("Seurat_object_total_cells.Rds")
MPcov
#An object of class Seurat 
#38397 features across 59138 samples within 2 assays 
#Active assay: SCT (18694 features, 3000 variable features)
# 1 other assay present: RNA
# 6 dimensional reductions calculated: pca, umap, umap_15, umap_20, umap_25, umap_30

## Fig1c复现
goi <- c("TIFAB","CD14", "CSF3R", "NRP1", 
         "CEBPE","CPA3","KLRK1","CD3E","TRDC","CD8A","CD4",
         "CD79B","IGHG4","MKI67","WFDC2", "HBA2","EPCAM")
#细胞中的markers汇总
#参考”Supplement.pdf“文件中的Supplementary Fig. 3
Idents(MPcov) <- "Annotation"
p <- DotPlot(
  MPcov,
  features = rev(goi),
  cols = c("lightgrey", "black"),
  col.min = -2.5,
  col.max = 2.5,
  dot.min = 0,
  dot.scale = 6,
  group.by = NULL,
  split.by = NULL,
  scale.by = "size",
  scale.min = NA,
  scale.max = NA
)
p$data$id = factor(p$data$id, levels(MPcov$Annotation) %>% rev)

pdf("Fig1c.dotplot.pdf",7.5*1.2,5*1.2)
p + theme(axis.text.x = element_text(angle = 90, hjust = 1))
dev.off()

1642518381084

## Fig1d复现

pdf("Fig1d.whole_umap_seurat_clusters.pdf",900/72*0.8, 688/72*0.8) 
DimPlot(MPcov,group.by = "Annotation", label=T)
dev.off()


## Fig1e whole proportion

colorder = c("Ctrl-1","Ctrl-2","Ctrl-3",
             "C2-1","C2-2","C2-3",
             "C5-1","C5-2", "C5-3","C5-4")

x <- table(MPcov$Annotation,MPcov$Sample_name)
x <- x[, colorder]
x3= t(t(x)/rowSums(t(x)))

x4 = as.data.frame(as.table(t(x3)))
colnames(x4) = c("sample","celltype","Freq")
x4$group = x4$sample %>% str_replace("-.*","")
x4$group = factor(x4$group, levels = c("Ctrl","C2","C5"))


write.csv(x4, "Fig1e.proportion_each_cluster.csv")

1642518734151

## Fig1d复现
dose<-read.csv( "Fig1e.proportion_each_cluster.csv")
which(dose$celltype=="RBC")
which(dose$celltype=="Doublet")
dose<-dose[-c(121:140),]
top<-function(x){
  return(mean(x)+sd(x)/sqrt(length(x)))
}
bottom<-function(x){
  return(mean(x)-sd(x)/sqrt(length(x)))
}
dose_Ctrl<-dose[which(dose$group=="Ctrl"),]
dose_C2<-dose[which(dose$group=="C2"),]
dose_C5<-dose[which(dose$group=="C5"),]

ggplot(data=dose_Ctrl,aes(x=celltype,y=Freq,fill=celltype))+
  stat_summary(geom = "bar",fun = "mean",
               position = position_dodge(0.9))+
  stat_summary(geom = "errorbar",
               fun.min = bottom,
               fun.max = top,
               position = position_dodge(0.9),
               width=0.2)+
  scale_y_continuous(expand = expansion(mult = c(0,0.1)))+
  theme_bw()+
  theme(panel.grid = element_blank())+
  labs(x="Celltype",y="Proportion")+
  geom_point(data=dose_Ctrl,aes(celltype,Freq),size=3,pch=19) 


ggplot(data=dose_C2,aes(x=celltype,y=Freq,fill=celltype))+
  stat_summary(geom = "bar",fun = "mean",
               position = position_dodge(0.9))+
  stat_summary(geom = "errorbar",
               fun.min = bottom,
               fun.max = top,
               position = position_dodge(0.9),
               width=0.2)+
  scale_y_continuous(expand = expansion(mult = c(0,0.1)))+
  theme_bw()+
  theme(panel.grid = element_blank())+
  labs(x="Celltype",y="Proportion")+
  geom_point(data=dose_C2,aes(celltype,Freq),size=3,pch=19) 


ggplot(data=dose_C5,aes(x=celltype,y=Freq,fill=celltype))+
  stat_summary(geom = "bar",fun = "mean",
               position = position_dodge(0.9))+
  stat_summary(geom = "errorbar",
               fun.min = bottom,
               fun.max = top,
               position = position_dodge(0.9),
               width=0.2)+
  scale_y_continuous(expand = expansion(mult = c(0,0.1)))+
  theme_bw()+
  theme(panel.grid = element_blank())+
  labs(x="Celltype",y="Proportion")+
geom_point(data=dose_C5,aes(celltype,Freq),size=3,pch=19)