动脉粥样硬化斑块的单细胞CITE-seq数据分析
考虑到咱们生信技能树粉丝对单细胞数据挖掘的需求，我开通了一个专栏《100个单细胞转录组数据降维聚类分群图表复现》，也亲自示范了几个，不过自己带娃，读博，时间精力有限，所以把剩余的90多个任务安排了学徒，实习生，学员。真的是太棒了，群策群力！
今天介绍的文章是2019年10月发表在Nature Medicine杂志《Single-cell immune landscape of human atherosclerotic plaques》文章链接是：
该研究主要发现AS患者的斑块中存在CD4+ T细胞的独特亚群以及被激活和分化的T细胞群，且与斑块易损性相关亚群的巨噬细胞表现出选择性激活的表型。另外，硬化斑块会高表达PD-1（程序性细胞死亡蛋白1）而使T细胞发生耗竭后的重编程。最后，作者提出对于患有潜在心血管疾病的癌症患者，PD-1抑制剂的使用可能增加心血管疾病风险。
这篇文章主要是进行免疫分析，采用了CITE-seq和scRNA-seq。数据存储在：Single-Cell-Immune-Profiling-of-Atherosclerotic-Plaques (figshare.com)，网页链接是：https://figshare.com/s/c00d88b1b25ef0c5c788
CITE-seq中的adt为蛋白序列；gex为基因表达矩阵。只提供了T细胞和巨噬细胞
下面是scRNA-seq数据；
本次要复现的文献图表如下所示：
本次数据分析是血液和斑块组织CITE-seq分析中免疫细胞分布，并不涉及scRNA-seq数据。
下面我们开始复现：
在运行环境下创建一个data_used文件夹，里面包含如下文件：
[1] "citeseq_pbmc_adt_macrophages.txt.gz" "citeseq_pbmc_adt_t_cells.txt.gz"
[3] "citeseq_plaque_adt_macrophages.txt.gz" "citeseq_plaque_adt_t_cells.txt.gz"
分别是外周血和斑块中的T细胞和巨噬细胞，也就是说是4个txt文本文件。
step1：读取抗体的序列文件，然后构建Seurat对象
把上面的4个文件批量读取，然后构建Seurat对象后保存为Rdata 文件，供后续使用
rm(list=ls())
suppressPackageStartupMessages(library(Seurat))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(clustree))
suppressPackageStartupMessages(library(cowplot))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(patchwork))
# 1-----
path = "data_used/"
fs=list.files(path = 'data_used/',
pattern = 'adt')
sceList = lapply(fs, function(f){
adt <- read.table( file.path('data_used',f))
adt[1:4,1:4]
sce = CreateSeuratObject(counts = adt)
sce$group=strsplit(f,'_')[[1]][2]
sce$celltype=gsub('.txt.gz','',strsplit(f,'_')[[1]][4])
sce
})
# pbmc和plaque来源的两种细胞
sce_all <- merge(sceList[[1]], y= sceList[-1] )
#备份数据
save(sce_all,file="sce_all.Rdata")
step2：降维聚类分群
仍然是可以参考前面的例子：人人都能学会的单细胞聚类分群注释 ，很简单的几个函数连起来使用即可！
# 2.预处理数据 ------
sce=sce_all
sce <- NormalizeData(sce)
sce = FindVariableFeatures(sce)
sce = ScaleData(sce, features = rownames(sce))
# 3.降维分析-----
# 3.1 PCA降维 ------
sce = RunPCA(sce, npcs = 50)
#3.2 umap和tsne降维 ------
sce = RunTSNE(sce, dims = 1:20)
sce = RunUMAP(sce, dims = 1:20)
# 4. 细胞聚类分析 -----
# 4.1 聚类
sce=FindNeighbors(sce, dims = 1:20, prune.SNN = 1/15)
#clustering 分析（设置不同的分辨率）
for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1)) {
sce=FindClusters(sce, graph.name = "RNA_snn", resolution = res, algorithm = 1)
}
# 4.2 可视化 ------
# 4.2.1 树 ------
(p1_tree=clustree(sce@meta.data, prefix = "RNA_snn_res."))
可以看到，最开始是一个大群，然后首先划分为两个亚群，其中左边的亚群可以继续细分，但是右边的亚群就停滞了。
# 4.2.2 res = 1 umap可视化 ------
sce <- SetIdent(sce, value = "RNA_snn_res.1")
#UMAP可视化
(p2 <- DimPlot(sce , reduction = "umap",label = T, label.box = T) + NoAxes())
(p3 <- DimPlot(sce, reduction = "umap", group.by = "celltype",
label = T, label.box = T) + NoAxes())
简单的给这些细胞亚群一个生物学命名，可以看到，比较小的那群是巨噬细胞，也就是没办法继续细分的。
step3：细胞亚群的生物学注释
文章对T细胞进行了一定程度的细分，代码如下所示：
#5 注释-----
genes_to_check <- rownames(sce)
ct = as.data.frame(sce@assays$RNA@counts)
(p4 <- DotPlot(sce, col.min=0, dot.min = 0.9,
features = unique(genes_to_check),assay='RNA') + coord_flip())
#命名
celltype=data.frame(ClusterID=0:12,celltype='unkown')
celltype[celltype$ClusterID == 3,2]='macrophages'
celltype[celltype$ClusterID %in% c(2,4,5,6,8,9,10,11),2]='CD4+'
celltype[celltype$ClusterID %in% c(0,7,12),2]='CD8+'
celltype[celltype$ClusterID == 1,2]='CD4+CD8+'
#将celltype信息添加值meta.data
sce@meta.data$CELLTYPE = "NA"
for(i in 1:nrow(celltype)){
sce@meta.data[which(sce@meta.data$RNA_snn_res.1 == celltype$ClusterID[i]),'CELLTYPE'] <- celltype$celltype[i]
}
# 细胞注释结果umap可视化 ------
(p5 <- DimPlot(sce, reduction = "umap", group.by = "CELLTYPE",
label = T, label.box = T) + NoAxes())
写在文末
咱们现在这个专栏《100个单细胞转录组数据降维聚类分群图表复现》分享的代码是到此为止，但是一般来说单细胞文章数据分析还有很多进阶图表制作，比如inferCNV看肿瘤拷贝数变异，monocle看拟时序等等。如果你也需要，可以加入我们这个专栏《100个单细胞转录组数据降维聚类分群图表复现》创作团队，获取进阶指引哦！见：急！计划招募100个单细胞爱好者，免费学全套单细胞降维聚类分群和生物学亚群注释