Seurat教程 || 分析Cell Hashing数据

时间:2022-07-23
本文章向大家介绍Seurat教程 || 分析Cell Hashing数据，主要内容包括其使用实例、应用技巧、基本知识点总结和需要注意事项，具有一定的参考价值，需要的朋友可以参考一下。

作者 | 周运来

Cell Hashing是与NYGC的技术创新小组合作开发的，它使用低聚标记（oligo-tagged）抗体来标记表面蛋白，在每个单细胞上放置一个“样本条形码（sample barcode）”，使得不同的样本可以被多路复用并在单个实验中运行。更多信息，可以参阅Cell Hashing with barcoded antibodies enables multiplexing and doublet detection for single cell genomics。（https://genomebiology.biomedcentral.com/articles/10.1186/s13059-018-1603-1）

其实我们在去年十月份的时候就关注过这个技术：Cell Hashing||单细胞多模态分析（https://www.jianshu.com/p/6ae3cc09d335）。案牍劳形，一直也没有再写它的文章了。今天我们依然跟着Seurat的官网来看看这个是如何分析的。请注意，如果需要看教程，请看官网。一则受个人的能力所限，一则官网才是更新快和及时的。

这里将简要演示如何在Seurat中处理由Cell Hashing 生成的数据。应用于两个数据集，我们可以成功地将细胞分离到它们的原始样本，并识别出交叉样本双峰（cross-sample doublets）。

多路复用函数HTODemux()实现了以下过程:

  • 我们对标准化的HTO值执行K -medoid聚类，它最初将细胞分离为K(# of samples )+1聚类。
  • 我们计算了HTO的一个“负”分布。对于每个HTO，我们使用平均值最低的群作为背景组。
  • 对于每个HTO，我们对负的聚类拟合一个负的二项分布。我们使用这个分布的0.99分位数作为阈值。
  • 根据这些阈值，每个细胞被划分为阳性或阴性的HTO。
  • 对于一个以上hto呈阳性的细胞被标注为双峰。

数据集描述

  • 数据代来自个不同献血者的外周血单核细胞(PBMCs)。
  • 每个供体的细胞都使用CD45作为哈希抗体进行唯一标记。
  • 随后在10X Chromium v2系统的单一lane上运行。
  • 你可以在这里（https://www.dropbox.com/sh/ntc33ium7cg1za1/AAD_8XIDmu4F7lJ-5sp-rGFYa?dl=0）下载RNA和HTO的计数矩阵，或者从GEO（https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE108313）下载FASTQ文件

原教程用的是他们处理好的rds文件，而这个我并没有成功下载，就从GEO中下载了表达谱，自己来构建Seurat的对象，所以会有所不同。这里要注意HTO数据一般和RNA的数据是对应的，对于RNA的我们很熟悉了，但是HTO的数据可能并不熟悉，这就要求我们看看 CITE-seq-Count（https://github.com/Hoohm/CITE-seq-Count） 的数据处理过程及其输出格式。

CITE-Seq 一般的序列结构：

For CITE-seq-Count, the output looks like this:

OUTFOLDER/
-- umi_count/
-- -- matrix.mtx.gz
-- -- features.tsv.gz
-- -- barcodes.tsv.gz
-- read_count/
-- -- matrix.mtx.gz
-- -- features.tsv.gz
-- -- barcodes.tsv.gz
-- unmapped.csv
-- run_report.yaml

好了我们了解了这个格式其实和10X的基本是一致的，这样我们就不用担心了，先读hto数据来看看。请注意，我们下载的是表达谱。

library(Seurat)
library(readr)

hto<- read_csv("GSM2895283_Hashtag-HTO-count.csv.gz")
Parsed with column specification:
cols(
  .default = col_double(),
  X1 = col_character()
)
See spec(...) for full column specifications.
|=================================================================| 100%    2 MB

hto[1:4,1:4]
# A tibble: 4 x 4
  X1                  GGCGACTAGAGGACGG CATCAAGGTCTTGTCC AAACCTGAGTGATCGG
  <chr>                          <dbl>            <dbl>            <dbl>
1 BatchA-AGGACCATCCAA               30                4               12
2 BatchB-ACATGTTACCGT               16               39               15
3 BatchC-AGCTTACTATCC               26                0               19
4 BatchD-TCGATAATGCGA             2698               22                2
> yhto <- as.data.frame(hto)

rownames(yhto) <- yhto[,1] ; yhto <- yhto[,-1] # 整理行名
yhto[1:4,1:4]
                    GGCGACTAGAGGACGG CATCAAGGTCTTGTCC AAACCTGAGTGATCGG
BatchA-AGGACCATCCAA               30                4               12
BatchB-ACATGTTACCGT               16               39               15
BatchC-AGCTTACTATCC               26                0               19
BatchD-TCGATAATGCGA             2698               22                2
                    TGAGGGAGTACTTAGC
BatchA-AGGACCATCCAA               26
BatchB-ACATGTTACCGT               20
BatchC-AGCTTACTATCC               13
BatchD-TCGATAATGCGA               24

rownames(yhto)
# BatchA-AGGACCATCCAA，BatchB-ACATGTTACCGT，BatchC-AGCTTACTATCC，
#  BatchD-TCGATAATGCGA，BatchE-GAGGCTGAGCTA，BatchF-GTGTGACGTATT，
#  BatchG-ACTGTCTAACGG，BatchH-TATCACATCGGT

# 还有特殊的三行：
# bad_struct;no_match;total_reads

同样的，我们读入RNA的umi数据。

umi  <-  read_tsv("GSM2895282_Hashtag-RNA.umi.txt.gz") 

Parsed with column specification:
cols(
  .default = col_double(),
  GENE = col_character()
)
See spec(...) for full column specifications.
|=================================================================| 100% 3902 MB
yumi <- as.data.frame(umi)
rownames(yumi) <- yumi[,1]
yumi<-yumi[,-1]
 dim(yumi)
[1] 40899 50000

仅仅是为了好抄代码，我们改一下变量名。

pbmc.umis <- yumi
 pbmc.htos <- yhto

rm(yumi)  #  释放我可怜的内存
rm(yhto)
gc()

取出RNA和HTO共有的barcode.

# Select cell barcodes detected by both RNA and HTO In the example datasets we have already
# filtered the cells for you, but perform this step for clarity.
joint.bcs <- intersect(colnames(pbmc.umis), colnames(pbmc.htos))

length(joint.bcs)
[1] 39842

# Subset RNA and HTO counts by joint cell barcodes
pbmc.umis <- pbmc.umis[, joint.bcs]
pbmc.htos <- as.matrix(pbmc.htos[, joint.bcs])

# Confirm that the HTO have the correct names
rownames(pbmc.htos)

[1] "BatchA-AGGACCATCCAA" "BatchB-ACATGTTACCGT" "BatchC-AGCTTACTATCC"
 [4] "BatchD-TCGATAATGCGA" "BatchE-GAGGCTGAGCTA" "BatchF-GTGTGACGTATT"
 [7] "BatchG-ACTGTCTAACGG" "BatchH-TATCACATCGGT" "bad_struct"         
[10] "no_match"            "total_reads"

创建Seurat对象，并添加HTO数据。

# Setup Seurat object
pbmc.hashtag <- CreateSeuratObject(counts = pbmc.umis)

# Normalize RNA data with log normalization
pbmc.hashtag <- NormalizeData(pbmc.hashtag)
# Find and scale variable features
pbmc.hashtag <- FindVariableFeatures(pbmc.hashtag, selection.method = "mean.var.plot")
pbmc.hashtag <- ScaleData(pbmc.hashtag, features = VariableFeatures(pbmc.hashtag))

添加HTO数据作为一个独立的assay.

# Add HTO data as a new assay independent from RNA
pbmc.hashtag[["HTO"]] <- CreateAssayObject(counts = pbmc.htos)
# Normalize HTO data, here we use centered log-ratio (CLR) transformation
pbmc.hashtag <- NormalizeData(pbmc.hashtag, assay = "HTO", normalization.method = "CLR")

head(pbmc.hashtag@meta.data)
                    orig.ident nCount_RNA nFeature_RNA nCount_HTO nFeature_HTO
TGCCAAATCTCTAAGG SeuratProject      84194         8167        350           11
ACTGCTCAGGTGTTAA SeuratProject      72747         7305       7704           11
TTCTTAGTCAACACGT SeuratProject      80714         8601        947           10
CATTCGCCACAGTCGC SeuratProject      71829         7226        759           11
TGGACGCCATCGACGC SeuratProject      68720         7500        444           11
CGAGAAGTCAACTCTT SeuratProject      63290         7395        597           10

看看数据的基本分布：

VlnPlot(pbmc.hashtag, features = c("nFeature_RNA", "nCount_RNA", "nCount_HTO","nFeature_HTO"), ncol = 2)

我们现在先不做过滤，看看是什么结果，然后再来思考我们的数据。看一下基于HTO富集的多复路细胞，这里使用Seurat函数HTODemux()函数将细胞分配回它们的样本上。

# If you have a very large dataset we suggest using k_function = 'clara'. This is a k-medoid
# clustering function for large applications You can also play with additional parameters (see
# documentation for HTODemux()) to adjust the threshold for classification Here we are using the
# default settings
pbmc.hashtag <- HTODemux(pbmc.hashtag, assay = "HTO", positive.quantile = 0.99)

# 下面时提示，不是代码，不要运行
Cutoff for BatchA-AGGACCATCCAA : 62 reads
Cutoff for BatchB-ACATGTTACCGT : 54 reads
Cutoff for BatchC-AGCTTACTATCC : 85 reads
Cutoff for BatchD-TCGATAATGCGA : 88 reads
Cutoff for BatchE-GAGGCTGAGCTA : 69 reads
Cutoff for BatchF-GTGTGACGTATT : 99 reads
Cutoff for BatchG-ACTGTCTAACGG : 157 reads
Cutoff for BatchH-TATCACATCGGT : 91 reads
Cutoff for bad-struct : 6 reads
Cutoff for no-match : 27 reads
Cutoff for total-reads : 362 reads

按照我们的惯例，我们看看这个函数的帮助文档：

?HTODemux 
HTODemux {Seurat}   R Documentation
Demultiplex samples based on data from cell 'hashing'
Description
Assign sample-of-origin for each cell, annotate doublets.

Usage
HTODemux(
  object,
  assay = "HTO",
  positive.quantile = 0.99,
  init = NULL,
  nstarts = 100,
  kfunc = "clara",
  nsamples = 100,
  seed = 42,
  verbose = TRUE
)

这时候我们看看数据metadata的变化，主要是这个函数做了什么。

head(pbmc.hashtag@meta.data)
                    orig.ident nCount_RNA nFeature_RNA nCount_HTO nFeature_HTO
TGCCAAATCTCTAAGG SeuratProject      84194         8167        350           11
ACTGCTCAGGTGTTAA SeuratProject      72747         7305       7704           11
TTCTTAGTCAACACGT SeuratProject      80714         8601        947           10
CATTCGCCACAGTCGC SeuratProject      71829         7226        759           11
TGGACGCCATCGACGC SeuratProject      68720         7500        444           11
CGAGAAGTCAACTCTT SeuratProject      63290         7395        597           10
                           HTO_maxID        HTO_secondID HTO_margin
TGCCAAATCTCTAAGG          bad-struct            no-match 0.24136778
ACTGCTCAGGTGTTAA BatchC-AGCTTACTATCC BatchB-ACATGTTACCGT 2.41473062
TTCTTAGTCAACACGT BatchF-GTGTGACGTATT          bad-struct 0.77349516
CATTCGCCACAGTCGC BatchD-TCGATAATGCGA BatchH-TATCACATCGGT 0.24264615
TGGACGCCATCGACGC          bad-struct BatchH-TATCACATCGGT 0.32077486
CGAGAAGTCAACTCTT BatchF-GTGTGACGTATT BatchE-GAGGCTGAGCTA 0.02388553
                                      HTO_classification
TGCCAAATCTCTAAGG                     bad-struct_no-match
ACTGCTCAGGTGTTAA BatchB-ACATGTTACCGT_BatchC-AGCTTACTATCC
TTCTTAGTCAACACGT          bad-struct_BatchF-GTGTGACGTATT
CATTCGCCACAGTCGC                     BatchD-TCGATAATGCGA
TGGACGCCATCGACGC                              bad-struct
CGAGAAGTCAACTCTT BatchE-GAGGCTGAGCTA_BatchF-GTGTGACGTATT
                 HTO_classification.global             hash.ID
TGCCAAATCTCTAAGG                   Doublet             Doublet
ACTGCTCAGGTGTTAA                   Doublet             Doublet
TTCTTAGTCAACACGT                   Doublet             Doublet
CATTCGCCACAGTCGC                   Singlet BatchD-TCGATAATGCGA
TGGACGCCATCGACGC                   Singlet          bad-struct
CGAGAAGTCAACTCTT                   Doublet             Doublet

运行HTODemux()的输出保存在对象元数据(metadata)中。我们可以看到有多少细胞被分为单细胞、双细胞和阴性/模糊细胞。

> # Global classification results
> table(pbmc.hashtag$HTO_classification.global)

 Doublet Negative  Singlet 
   22650    14520     2672      
# 可见，我们的数据是需要过滤的了，Singlet 居然是最少的，为什么？

用山脊图对选定的hto进行可视化富集的结果。

# Group cells based on the max HTO signal
Idents(pbmc.hashtag) <- "HTO_maxID"
RidgePlot(pbmc.hashtag, assay = "HTO", features = rownames(pbmc.hashtag[["HTO"]])[1:2], ncol = 2)

将HTO信号对象化，以确认singlets中的互斥性.

 levels(pbmc.hashtag)
 [1] "bad-struct"          "BatchA-AGGACCATCCAA" "BatchB-ACATGTTACCGT"
 [4] "BatchC-AGCTTACTATCC" "BatchD-TCGATAATGCGA" "BatchE-GAGGCTGAGCTA"
 [7] "BatchF-GTGTGACGTATT" "BatchG-ACTGTCTAACGG" "BatchH-TATCACATCGGT"
[10] "no-match"          

FeatureScatter(pbmc.hashtag, feature1 = "hto_BatchA-AGGACCATCCAA", feature2 = "hto_BatchE-GAGGCTGAGCTA")

这个hto_BatchA-AGGACCATCCAA是哪里来的，或者它在哪里？

噢噢，我们注意到：

pbmc.hashtag@assays
$RNA
Assay data with 40899 features for 39842 cells
Top 10 variable features:
 mt-Nd1, Malat1, mt-Cytb, IGLC3, IGLC2, IGHA1, mt-Nd4, IGHG2, IGJ, IGHA2 

$HTO
Assay data with 11 features for 39842 cells
First 10 features:
 BatchA-AGGACCATCCAA, BatchB-ACATGTTACCGT, BatchC-AGCTTACTATCC,
BatchD-TCGATAATGCGA, BatchE-GAGGCTGAGCTA, BatchF-GTGTGACGTATT,
BatchG-ACTGTCTAACGG, BatchH-TATCACATCGGT, bad-struct, no-match 

FeatureScatter(pbmc.hashtag, feature1 ="rna_Malat1" , feature2  = "rna_IGLC3")

从哪里取特征的标记是rna_和hto_。

比较Doublet Negative Singlet 细胞的UMIs数目。

Idents(pbmc.hashtag) <- "HTO_classification.global"
VlnPlot(pbmc.hashtag, features = "nCount_RNA", pt.size = 0.1, log = TRUE)

生成一个用于hto的二维UMAP嵌入。为了简单起见，我们在这里通过 singlets and doublets 进行分组。

# First, we will remove negative cells from the object
pbmc.hashtag.subset <- subset(pbmc.hashtag, idents = "Negative", invert = TRUE)

# Calculate a distance matrix using HTO
hto.dist.mtx <- as.matrix(dist(t(GetAssayData(object = pbmc.hashtag.subset, assay = "HTO"))))

# Calculate tSNE embeddings with a distance matrix
#pbmc.hashtag.subset <- RunTSNE(pbmc.hashtag.subset, distance.matrix = hto.dist.mtx, perplexity = 100)
# tsne 实在是跑不动，等不住，换umap吧
 pbmc.hashtag.subset <- RunPCA(pbmc.hashtag.subset)
pbmc.hashtag.subset <- RunUMAP(pbmc.hashtag.subset,dims=1:10)

DimPlot(pbmc.hashtag.subset)

计算一下各个标签下细胞的差异基因，推测它们的差异会带来什么影响。

pbmc.singlet@misc$mk <- FindAllMarkers(pbmc.,only.pos=T)

 pbmc.hashtag.subset@misc$mk %>% filter(cluster == 'bad-struct')  %>% top_n(10,avg_logFC)
               p_val avg_logFC pct.1 pct.2    p_val_adj    cluster    gene
CALD1   4.269893e-28  1.305611 0.213 0.051 1.746343e-23 bad-struct   CALD1
PTTG1   5.510165e-18  1.308636 0.230 0.075 2.253602e-13 bad-struct   PTTG1
PDHA1   5.757841e-18  1.331648 0.196 0.058 2.354900e-13 bad-struct   PDHA1
SNRNP25 9.296534e-18  1.427055 0.213 0.066 3.802189e-13 bad-struct SNRNP25
MDK     1.299665e-17  1.296846 0.213 0.067 5.315500e-13 bad-struct     MDK
HMGA1   3.069179e-14  1.396448 0.230 0.086 1.255263e-09 bad-struct   HMGA1
HSPB1   4.392379e-10  1.386510 0.248 0.110 1.796439e-05 bad-struct   HSPB1
HMGN2   4.458862e-10  1.458616 0.243 0.108 1.823630e-05 bad-struct   HMGN2
CCT3    1.635818e-09  1.396913 0.243 0.112 6.690333e-05 bad-struct    CCT3
TIMM8B  1.135882e-06  1.335207 0.230 0.117 4.645644e-02 bad-struct  TIMM8B

HTO_classification.global

DimPlot(pbmc.hashtag.subset,group.by="HTO_classification.global")

同样，看两者的差异基因如何。

pbmc.hashtag.subset@misc$mk2 <- FindMarkers(pbmc.hashtag.subset,only.pos=T,ident.1= 'Doublet',group.by="HTO_classification.global")


pbmc.hashtag.subset@misc$mk2  %>% top_n(10,avg_logFC)
                      p_val avg_logFC pct.1 pct.2     p_val_adj
CD52           0.000000e+00  1.370377 0.627 0.102  0.000000e+00
CD74          5.492774e-210  1.491273 0.379 0.079 2.246489e-205
HLA-DRA       8.923196e-155  1.686647 0.266 0.035 3.649498e-150
AIF1          2.244058e-139  1.590655 0.235 0.025 9.177971e-135
FOS           1.498210e-126  1.985102 0.235 0.037 6.127527e-122
LST1          2.400314e-119  1.501267 0.205 0.020 9.817044e-115
LYZ           7.770229e-119  2.013469 0.222 0.033 3.177946e-114
S100A9        1.729465e-101  1.765233 0.226 0.051  7.073339e-97
RP11-1143G9.4  1.297476e-90  1.634986 0.168 0.020  5.306549e-86
S100A8         3.501114e-90  1.719860 0.236 0.071  1.431920e-85


# You can also visualize the more detailed classification result by running Idents(object) <-
# 'HTO_classification' beofre plotting. Here, you can see that each of the small clouds on the
# tSNE plot corresponds to one of the 28 possible doublet combinations.

HTO热图

# To increase the efficiency of plotting, you can subsample cells using the num.cells argument
HTOHeatmap(pbmc.hashtag, assay = "HTO", ncells = 5000)

使用通常的scRNA-seq工作流和可视化细胞，并检查潜在的批次效应。

# Extract the singlets

Idents(pbmc.hashtag)  <- "HTO_classification.global"
pbmc.singlet <- subset(pbmc.hashtag, idents = "Singlet")

# Select the top 1000 most variable features
pbmc.singlet <- FindVariableFeatures(pbmc.singlet, selection.method = "mean.var.plot")

# Scaling RNA data, we only scale the variable features here for efficiency
pbmc.singlet <- ScaleData(pbmc.singlet, features = VariableFeatures(pbmc.singlet))

# Run PCA
pbmc.singlet <- RunPCA(pbmc.singlet, features = VariableFeatures(pbmc.singlet))

# We select the top 10 PCs for clustering and tSNE based on PCElbowPlot
pbmc.singlet <- FindNeighbors(pbmc.singlet, reduction = "pca", dims = 1:10)
pbmc.singlet <- FindClusters(pbmc.singlet, resolution = 0.6, verbose = FALSE)
# pbmc.singlet <- RunTSNE(pbmc.singlet, reduction = "pca", dims = 1:10)
pbmc.singlet <- RunUMAP(pbmc.singlet, reduction = "pca", dims = 1:10)
# Projecting singlet identities on UMAPvisualization
DimPlot(pbmc.singlet, group.by = "HTO_classification")



DimPlot(pbmc.singlet)

接下来该找差异基因了。

pbmc.singlet@misc$mk <- FindAllMarkers(pbmc.singlet,only.pos=T)
Calculating cluster 0

  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04s  
Calculating cluster 1
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s  
Calculating cluster 2
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  
Calculating cluster 3
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  
Calculating cluster 4
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04s  
Calculating cluster 5
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  
Calculating cluster 6
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=18s  


 pbmc.singlet@misc$mk %>% filter(cluster == "6") %>% top_n(10,avg_logFC)
               p_val avg_logFC pct.1 pct.2     p_val_adj cluster   gene
Rps2    0.000000e+00  3.828804 1.000 0.008  0.000000e+00       6   Rps2
Rps18  1.774313e-299  3.590482 1.000 0.010 7.256762e-295       6  Rps18
Rplp0  5.575846e-297  3.671938 1.000 0.010 2.280465e-292       6  Rplp0
Rps27  3.979973e-290  3.524766 1.000 0.011 1.627769e-285       6  Rps27
Rps8   3.735265e-256  3.513067 1.000 0.013 1.527686e-251       6   Rps8
Ftl1   6.754453e-253  3.719029 1.000 0.014 2.762504e-248       6   Ftl1
Rps19  8.025777e-177  3.566633 1.000 0.025 3.282463e-172       6  Rps19
Fth1   3.847843e-173  3.519539 1.000 0.025 1.573729e-168       6   Fth1
Lgals1 9.363616e-158  3.608453 1.000 0.029 3.829625e-153       6 Lgals1
Malat1 4.534014e-119  4.348673 0.929 0.036 1.854367e-114       6 Malat1

说一说，这一群的差异基因有什么特点。

https://satijalab.org/seurat/v3.2/hashing_vignette.html

随机文章
本站知识点必读
最近更新