跟着nature学习单细胞分析之CITE-seq：同时测序mRNA和蛋白的技术是CITE-seq技术

首页标题 ꄲ 跟着nature学习单细胞分析之CITE-seq：同时测序mRNA和蛋白的技术是CITE-seq技术

跟着nature学习单细胞分析之CITE-seq分析教程

通常单细胞数据集在ECA和EGA数据库上，需要和作者联系才可以下载，并且在GEO上只找到了一个符合心意的，但是居然是CITE-seq技术，并不是常规的单细胞测序技术。所以希望有专门的教程来介绍一下这个分析这方面的知识。

一、理论基础

在进行分析之前，我们有必要了解CITE-seq技术究竟是什么？

我们都知道，如果我们想要了解单个细胞的功能或者是活性，那么最直接的方法可能就是通过流式细胞技术来检测细胞蛋白表达，但是很显然这种方法是比较费时费力的，目前来说单细胞可以做到一次性检测上万个细胞的mRNA表达情况，但是却无法对每一个细胞的蛋白表达予以计算我们知道生物体功能的执行者是蛋白，而且从mRNA到蛋白会受到转录后调控，也就是说单纯从mRNA表达量推测功能并不是十分的准确，所以我们有必要同时去考虑蛋白和mRNA的表达量（在这里蛋白被称呼为ADT，即antibody-derived tags，可以近似的认为ADT就是抗体所对应的蛋白在一个细胞上的表达量）

问题：同时测序mRNA和蛋白的技术是CITE-seq技术，那么对不同样本同时进行测序的技术是什么呢？

按照ADT进行细胞划分究竟准不准确呢，下面这张图或许可以给我们答案：

可以很清楚的看到，运用scRNA-seq进行细胞划分与运用CITE-seq中ADT进行划分可以达到相似的效果（这里我们需要注意一下CITE-seq技术是由两部分组成，一部分是scRNA-seq，一部分是检测了ADT，所以这里就有一个小技巧，如果你感兴趣的数据集是CITE-seq技术，但是你又不会，那么你单纯去分析scRNA-seq的部分也是完全没有问题的）

具体的技术细节其实不需要我们了解，我们只需要知道以下内容：

1.CITE-seq技术可以一次性获得单个细胞的mRNA和蛋白的表达量（目前来说对于蛋白的数量倒是没有明确的限制，但是一次性越多数量那么价格自然越高，所以目前来说常见的数量是100-200左右）

2.单细胞测序最大的优点就是可以实现计算单个细胞的表达量，避免了bulk-seq计算的平均表达量而导致的忽略小细胞亚群的结局，但是单细胞的难点也在于细胞分类，通过mRNA表达量进行分类可能并不是十分准确，因为按照功能来说，蛋白的表达量可能会更加直观，所以如果添加了蛋白的表达信息，就可以提高细胞亚群分类的精度，同时也可以帮助我们识别罕见的细胞亚群

那么既然我们已经对CITE-seq技术有一个大概的了解了，那么下面我们就进行代码实战

二、数据来源

接下来我们要使用的数据来自下面这篇文章

这篇文章研究癫痫疾病，作者选择了11个药物难治性癫痫（DRE）样本进行CITE-seq，目的是为了揭示癫痫的免疫特征，当然具体的生物学背景我们不是这个领域的就不要过多了解了，结论就是癫痫的发病机制与神经炎症以及促炎微环境有关，然后作者也是贴心的提供了数据下载链接，我们接下来就来尝试进行一下分析

1 #准备工作

2 library(Seurat)

3 library(tidyverse)

4 library(dplyr)

5library(patchwork)

6 ##系统报错改为英文

7 Sys.setenv(LANGUAGE = "en")

8 ##禁止转化为因子

9 options(stringsAsFactors = FALSE)

10 ##清空环境

11 rm(list=ls())

工欲善其事必先利其器~

这里关于样本的信息，我们需要在做数据分析前重点明确，因为流程大部分都是标准化的，但是找寻数据的本领确是每个人所独有的，所以如果各位小伙伴并不知道下面这张可视化所代表的样本信息是什么，建议各位小伙伴可以花费一段时间去看一下GEO相关数据集或者这篇文献（OL：嗅叶；FL：额叶；TL：颞叶）

三、数据读取

1 #读取数据

2 #本篇文献一共有11个样本，这里我们为了提高效率（内存不够）以我们这里使用4个样本

3 dir = c('课题设计/项目2/GSE201048/Sample1/',

4'课题设计/项目2/GSE201048/Sample2/','课题设计/项目2/GSE201048/Sample3/', '课题设计/项目2/GSE201048/Sample4/')

5 names(dir) = c('Sample1', 'Sample2',"Sample3","Sample4") #多样本整合的时候命名barcode-ID

6 counts <- Read10X(data.dir =dir)

1 #如果我们分别有两个文件，分别提供scRNA-seq的数据以及ADT的数据，那么我们可以分别读取然后对S4对象进行添加

2 #counts <- CreateSeuratObject(counts = counts)

3#adt_assay <- CreateAssayObject(counts = adt)

4 #counts[["ADT"]] <- adt_assay

这里我们可以很清楚的看到，我们这里有两个文件，分别是表达数据以及ADT数据，也就是说这两个数据结合起来其实就是我们的CITE-seq技术

1 #但是尽管读取是一起读取并创建了counts对象，但是我们后续创建Seurat对象的时候是需要分别读取的

2 scRNA <- CreateSeuratObject(counts[["Gene Expression"]])

3 adt <- CreateAssayObject(counts[["Antibody Capture"]])

4 all.equal(colnames(counts[["Gene Expression"]]), colnames(counts[["Antibody Capture"]]))

5 scRNA[["ADT"]] <- adt

6 Assays(scRNA)

7 #[1] "RNA" "ADT"

all.equal函数的目的就是我们每一个细胞都是分别检测了mRNA的表达量以及蛋白的表达量，所以不管是顺序还是名称都应该是一样的，然后我们把读取的蛋白表达信息添加到Seurat对象之中，这样子，我们就制备好了输入数据

1 #查看数据

2 rownames(scRNA[["ADT"]])#可以获得作者都测了哪些蛋白

3 #[1] "CD3-TotalSeqB" "CD4-TotalSeqB" "CD8a-TotalSeqB"

4 #[4] "CD14-TotalSeqB" "CD16-TotalSeqB" "CD19-TotalSeqB"

5 #[7] "CD25-TotalSeqB" "CD56-TotalSeqB" "CD20-TotalSeqB"

6 #[10] "CD335-TotalSeqB" "CD69-TotalSeqB" "CD197-TotalSeqB"

7 #[13] "CD27-TotalSeqB" "HLA-DR-TotalSeqB" "CD11b-TotalSeqB"

8 #[16] "CD45-TotalSeqB" "prip-TotalSeqB" "n212-TotalSeqB"

1 #注意：我们可以通过函数快速地在两个数据之间来回切换，以便后续的分析以及可视化

2 DefaultAssay(scRNA)

2 #[1] "RNA"

可以很清楚的看到默认是RNA槽的数据，可以进行后续的分析

四、scRNA-seq数据

这里怕各位小伙伴对于单细胞这些维度数据有点遗忘，推荐大家可以看下面的推文：

毕业救星！这个R包在高分文章常见，实用！好学！

1 #scRNA@assays[["RNA"]]@counts#原始数据，就是我们读取那三个文件后获得的原始数据

2 #scRNA@assays[["RNA"]]@data#经过标准化的数据

3 #scRNA@assays[["RNA"]]@scale.data#经过标准化后执行高可变基因筛选并且归一化后的数据

然后关于后续分析处理细节的参数调整，我们则是按照文献中作者的要求来进行调节：

1 #然后接下来执行

2 scRNA-seq的标准流程scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")

3 VlnPlot(scRNA, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

4 scRNA <- subset(scRNA, subset = nFeature_RNA > 300 & nFeature_RNA < 5000 & percent.mt < 20)
5

6 scRNA <- NormalizeData(scRNA)

7 scRNA <- FindVariableFeatures(scRNA,nfeatures = 2000)

8 scRNA <- ScaleData(scRNA)

9 scRNA <- RunPCA(scRNA,dims = 1:20)
10

11 scRNA <- FindNeighbors(scRNA, dims = 1:20)

12 scRNA <- FindClusters(scRNA, resolution = 0.9, verbose = FALSE)

13 scRNA <- RunUMAP(scRNA, dims = 1:20)

14 DimPlot(scRNA, label = TRUE)
15

16 #寻找差异基因

17 scRNA.markers <- FindAllMarkers(scRNA, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

18 scRNA.markers %>%

19 group_by(cluster) %>%

20 slice_max(n = 10, order_by = avg_log2FC)

然后我们可以得到下面的可视化以及一个差异基因的scRNA.markers对象：

五、ADT数据

然后scRNA-seq的处理就先暂时结束，我们接下来继续处理ADT数据的部分：

1 # Normalize ADT data

2 DefaultAssay(scRNA) <- "ADT"

3 scRNA <- NormalizeData(scRNA, normalization.method = "CLR", margin = 2)
4

5 #visualize

6 rownames(scRNA[["ADT"]])

7 p1 <- FeaturePlot(scRNA, "CD45-TotalSeqB", cols = c("lightgrey", 6 "darkgreen")) + ggtitle("ADT")

8 DefaultAssay(scRNA) <- "RNA"

9 p2 <- DimPlot(scRNA, label = TRUE)+ggtitle("RNA")

10 p1|p2

这里我们需要知道，作者通过CD45蛋白的表达来区分免疫细胞和非免疫细胞，其实我们需要知道CITE-seq之所以要引入蛋白的表达量（ADT）其实就是为了帮助我们更好的对细胞亚群进行注释同时也可以更加精准的帮助我们划分细胞亚群，找到罕见的细胞亚群

六、数据整合

然后上述流程只不过依旧是单打单的行为，我们接下来要真正尝试去把两种数据结合起来，为此我们需要使用到WNN算法，这个算法就可以把多模态的数据整合起来

1 #设置并行加快运行速度

2 library(future)

3 plan("multiprocess", workers =4)#设置线程

4 options(future.globals.maxSize = 2000 * 1024^2)
5

6 scRNA_WNN <- CreateSeuratObject(counts[["Gene Expression"]])

7 adt_WNN <- CreateAssayObject(counts[["Antibody Capture"]])

8 # add this assay to the previously created Seurat object

9 scRNA_WNN[["ADT"]] <- adt_WNN
10

11 DefaultAssay(scRNA_WNN) <- 'RNA'

12 scRNA_WNN[["percent.mt"]] <- PercentageFeatureSet(scRNA_WNN, pattern = "^MT-")

13 scRNA_WNN <- subset(scRNA_WNN, subset = nFeature_RNA > 300 & nFeature_RNA < 5000 & percent.mt < 20)

14 scRNA_WNN <- NormalizeData(scRNA_WNN) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA()
15

16 DefaultAssay(scRNA_WNN) <- 'ADT'

17 VariableFeatures(scRNA_WNN) <- rownames(scRNA_WNN[["ADT"]])

18 scRNA_WNN <- NormalizeData(scRNA_WNN, normalization.method = 'CLR', margin = 2) %>%

19 ScaleData() %>% RunPCA(reduction.name = 'apca')

这里需要给各位小伙伴解释一下reduction.name参数为什么需要添加，RunPCA函数执行PCA结果，然后通过reduction.name参数我们可以设置结果的名称，默认是叫做PCA，但是因为我们这里是多模态数据，而且我们是针对同一个S4对象，所以我们需要添加不同的名称加以区分