Seurat-多组学数据整合

R版本与运行环境信息

  1. Date:2021-5-6
  2. R version 4.0.3 (2020-10-10)
  3. Platform: x86_64-conda-linux-gnu (64-bit)
  4. Running under: CentOS Linux 7 (Core)

多组学数据

是指对同一个细胞同时测量多种指标， 如CITE-seq ：是一种可以同时检测细胞表面蛋白和RNA的单细胞测序。一方面用带DNA标签的抗体来检测细胞表面蛋白，同时用10x的RNA单细胞测序来检测全部的mRNA。scRNA-seq+scATAC-seq，本例子主要是演示了如何在单细胞转录组cluster后的基础上去分析每个cell上的表面蛋白的含量

数据下载

Here, we analyze a dataset of 8,617 cord blood mononuclear cells (CBMCs), where transcriptomic measurements are paired with abundance estimates for 11 surface proteins,

  1. #转录组数据
  2. wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE100nnn/GSE100866/suppl/GSE100866_CBMC_8K_13AB_10X-RNA_umi.csv.gz
  3. #表面蛋白数据
  4. wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE100nnn/GSE100866/suppl/GSE100866_CBMC_8K_13AB_10X-ADT_umi.csv.gz

载入相关包并数据读入

由于数据量太大，所以上服务器

  1. library(Seurat)
  2. library(tidyverse)
  3. library(cowplot)
  4. setwd("~/multi_data/")
  5. #查看文件列表，两个文件
  6. #'GSE100866_CBMC_8K_13AB_10X-ADT_umi.csv.gz'：表面蛋白的数据
  7. #'GSE100866_CBMC_8K_13AB_10X-RNA_umi.csv.gz'：转录组的数据
  8. #读入转录组的数据
  9. cbmc.rna <- as.sparse(read.csv("GSE100866_CBMC_8K_13AB_10X-RNA_umi.csv.gz",header = T,row.names = 1))
  10. #仅仅提取人类的数据，CollapseSpeciesExpressionMatrix，当数据是由多物种组成的时候，可以提取某物种的数据，本数据分为人和老鼠，人的数据默认为HUMAN_老鼠尾MOUSE_，提取人类的前100的表达量基因
  11. cbmc.rna <- CollapseSpeciesExpressionMatrix(cbmc.rna,prefix = "HUMAN_",controls = "MOUSE_",ncontrols = 100)
  12. #读入表面蛋白的数据
  13. cmbc.adt <- as.sparse(read.csv("GSE100866_CBMC_8K_13AB_10X-ADT_umi.csv.gz",header=T,row.names = 1))
  14. #检查二者的cell是否一致
  15. >all.equal(colnames(cbmc.rna),colnames(cmbc.adt))
  16. TRUE

创建Seurat对象并将表面蛋白数据加入

  1. #创建seurat对象
  2. cbmc <- CreateSeuratObject(counts = cbmc.rna,project = "multi_test")
  3. #将表面蛋白数据创建为seurat对象CreateAssayObject,并将其写入cbmc
  4. adt_assay <- CreateAssayObject(counts = cmbc.adt)
  5. cbmc[["ADT"]] <- adt_assay
  7. ###查看cbmc内的数据
  8. >Assays(cbmc)
  9. 'RNA''ADT'
  10. >cbmc[["ADT"]]
  11. Assay data with 13 features for 8617 cells
  12. First 10 features:
  13.  CD3, CD4, CD8, CD45RA, CD56, CD16, CD10, CD11c, CD14, CD19 
  14. >cbmc[["RNA"]]
  15. Assay data with 36240 features for 8617 cells
  16. First 10 features:
  17.  ERCC-ERCC-00104, HUMAN-A1BG, HUMAN-A1BG-AS1, HUMAN-A1CF, HUMAN-A2M,
  18. HUMAN-A2M-AS1, HUMAN-A2ML1, HUMAN-A4GALT, HUMAN-A4GNT, HUMAN-AAAS

对转录组数据分析

  1. #设定分析数据
  2. DefaultAssay(cbmc) <- "RNA"
  3. #对数据标准分析，质控 归一化 高变基因筛选 标准化 降维 聚类.....
  4. cbmc <- NormalizeData(cbmc,normalization.method = "LogNormalize",scale.factor = 10000,verbose = T)
  5. cbmc <- FindVariableFeatures(cbmc,selection.method = "vst",nfeatures = 2000,verbose = T)
  6. cbmc <- ScaleData(cbmc,features = VariableFeatures(cbmc),verbose = T)
  7. cbmc <- RunPCA(cbmc,verbose = T)
  8. cbmc <- FindNeighbors(cbmc,dims = 1:30,verbose = T)
  9. cbmc <- FindClusters(cbmc,resolution = 0.8,verbose = T)
  10. cbmc <- RunUMAP(cbmc,dims = 1:30,verbose = T)
  11. p1 <- DimPlot(cbmc,label = T)
  12. p1

scRNA-Seurat-多组学数据整合 - 图1

表面蛋白数据的可视化

现在已经得到了转录组的结果，因此下一步可以通过对表面蛋白的数据进行分析，首先需要对表面蛋白的数据进行归一化处理

  1. #"CD19" %in% rownames(cbmc) 
  2. #查看CD19 RNA在cell中的表达
  3. DefaultAssay(cbmc) <- "RNA"
  4. cd19.rna <- FeaturePlot(cbmc,reduction = "umap",features = "CD19",cols = c("lightgrey", "darkgreen")) + ggtitle("CD19 RNA")
  5. #查看CD19 蛋白在cell中的表达
  6. DefaultAssay(cbmc) <- "ADT"
  7. cd19.pro <- FeaturePlot(cbmc,reduction = "umap",features = "CD19",cols = c("gray", "red")) + ggtitle("CD19 pro")
  8. cd19.rna + cd19.pror

scRNA-Seurat-多组学数据整合 - 图2

  1. #可以使用Key为ADT和RNA数据添加前缀
  2. Key(cbmc[["ADT"]])
  3. Key(cbmc[["RNA"]])
  4. ##返回结果为前缀
  5. 'adt_'
  6. 'rna_'
  7. ##无需切换数据即可查看CD19的表达情况
  8. p3 <- FeaturePlot(cbmc,features = "adt_CD19",cols = c("lightgrey", "darkgreen"))+ ggtitle("CD19 ADT")
  9. p4 <- FeaturePlot(cbmc,features = "adt_CD19")+ ggtitle("CD19 RNA")
  10. p3 + p4

scRNA-Seurat-多组学数据整合 - 图3

查看表面蛋白在不同的cluster上的表达情况

  1. #CD19 是B cell的常见marker
  2. p5 <- VlnPlot(cbmc,features = "adt_CD19")
  3. p5

scRNA-Seurat-多组学数据整合 - 图4

基因/蛋白在各个cluster的 表达情况

  1. #查看各个表面蛋白在不同cell中的表达情况
  2. adt.markers <- FindMarkers(cbmc,ident.1 = 6,assay = "ADT")
  4. >head(adt.markers)
  5. A data.frame: 7 × 5
  6. p_val    avg_log2FC    pct.1    pct.2    p_val_adj
  7. <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
  8. CD19    2.067533e-215    1.2787751    1    1    2.687793e-214
  9. CD45RA    8.106076e-109    0.4117172    1    1    1.053790e-107
  10. CD4    1.123162e-107    -0.7255977    1    1    1.460110e-106
  11. CD14    7.212876e-106    -0.5060496    1    1    9.376739e-105
  12. CD3    1.639633e-87    -0.6565471    1    1    2.131523e-86
  13. CD8    1.042859e-17    -0.3001131    1    1    1.355716e-16
  14. CD11c    8.957964e-11    -0.4382277    1    1    1.164535e-0
  15. #查看各个基因在不同cell中的表达情况
  16. rna.marker <- FindMarkers(cbmc,ident.1 = 5,assay = "RNA")
  18. >head(rna.marker)
  19. A data.frame: 6 × 5
  20. p_val    avg_log2FC    pct.1    pct.2    p_val_adj
  21. <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
  22. AC109351.1    0    0.3203893    0.265    0.005    0
  23. CTD-2090I13.1    0    2.0024376    0.972    0.062    0
  24. DCAF5    0    0.6637418    0.619    0.055    0
  25. DYNLL2    0    2.0387603    0.984    0.094    0
  26. FAM186B    0    0.3000479    0.244    0.002    0
  27. HIST2H2AB    0    1.3104432    0.812    0.013    0

查看各个基因/蛋白表达的相关性

  1. p5 <- FeatureScatter(cbmc,feature1 = "adt_CD19", feature2 = "adt_CD3")
  2. p6 <- FeatureScatter(cbmc, feature1 = "adt_CD3", feature2 = "rna_CD3E")
  3. p7 <- FeatureScatter(cbmc, feature1 = "adt_CD4", feature2 = "adt_CD8")
  4. #选择原始数据进行可视化"adt_CD4"和"adt_CD8"表达量之间的关系
  5. p8 <- FeatureScatter(cbmc, feature1 = "adt_CD4", feature2 = "adt_CD8", slot = "counts")
  6. plot_grid(p5,p6,p7,p8,labels = c("adt_CD19 vs adt_CD3","adt_CD3 vs rna_CD3E","adt_CD4 vs adt_CD8","no-Normalize"))

scRNA-Seurat-多组学数据整合 - 图5