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