pbmc方面中文杂志

pbmc方面中文杂志

随着单细胞测序技术的发展，多种组学的单细胞数据也越来越多，要如何更好的使用多组学的数据去解析样本的细胞组成和特征情况呢？今年4月份发表在Cell杂志上单细胞多模态数据的整合分析这篇文章中介绍了WNN（Weighted-nearest neighbor,加权最近邻）的算法，seurat团队使用不同的数据集对算法模型的构建、验证及应用进行了深入浅出的说明。在对文章进行说明前，首先来了解一下10XGenomics推出的一份样本获得两种组学的产品10XATAC_GEM双组学的原理。

单细胞转录组的优势在于可以发现新的细胞类群，但难以发现分子相似，功能不同的类群，例如T细胞中，RNA量少，RNA酶多，就难以区分亚类群，而此时多组学就有了更多的优势。10X单细胞双组学的原理如下图所示：

获取细胞核后，先利用转座酶试剂对其进行转座反应，对染色质开放区的DNA进行打断和片段化，单细胞分选，油包水液滴（GEM）中并被回收，随后细胞核开始裂解并释放出DNA片段及mRNA，在GEM液滴中完成逆转录反应，同时为DNA片段及cDNA标记上该液滴中Gelbead所带有的特定Barcode标签序列。最终构建出如下图所示的两种不同的文库包括单细胞核转录组文库和ATAC的文库：

这两种文库与单独的转录组和ATAC文库无异，仅在ATAC的index序列中增加了8bp的spacer序列。

那么问题来了如何整合这两个组学的数据呢？seurat团队给我们提供了一个有效的算法和思路。

多模态分析（multimodal analysis）就是同时测量单细胞的多模态数据，它代表了单细胞基因组学的一个发展方向，同时也需要基于多种数据类型的新的计算方法来描述细胞状态。文章介绍了“加权最近邻（weighted-nearest neighbor，WNN）”分析：使用一个无监督的框架来学习每个细胞中每种数据类型的相对效用，使多种模式的整合分析成为可能。将文章的算法应用于包含几十万个人类白细胞的CITE-seq数据集以及228个抗体的panel上，以构建一个循环免疫系统的多模态参考图谱。文章表明整合分析大大提高了描述细胞状态的能力，并验证了新的淋巴亚群的存在。此外，文章还演示了如何利用这一参考快速绘制新数据集，并解释免疫接种和COVID-19的免疫反应。

文章概览如下：

对人类免疫系统中丰富多样的细胞类型进行分类鉴定，对单细胞基因组学来说是一个有力的证明，但也展现出了他的局限性。虽然单细胞转录组 (scRNA-seq)能够发现异质组织中的新细胞类型和状态，但单靠转录组学常常无法分离分子上相似但功能上不同的免疫细胞类型。尽管T细胞具有功能多样性，但不同的T细胞群，如效应细胞、调节细胞、细胞内固定细胞和黏膜相关不变T细胞(MAIT)，即使使用最敏感和最尖端的技术，通常也不能仅用scRNA-seq有效地分离它们。

多模态单细胞技术，在同一细胞中同时描述多种数据类型，代表了细胞状态发现和鉴定的新前沿。例如，最近引入了CITE-seq，它利用寡核苷酸偶联抗体，通过测序抗体衍生标签(antibody-derived tags ，ADTs)，同时量化单细胞内RNA和表面蛋白的丰度。此外，随着技术进步，现在可以在染色质可及性（ATAC）、DNA甲基化、核小体占位（nucleosome occupancy ）或空间定位的同时对转录组进行分析。这些方法都提供了一个令人兴奋的解决方案，以克服scRNA-seq固有的局限性，并探索多种细胞模式如何影响细胞状态和功能。

在这里，文章引入了“加权最近邻”（weighted-nearest neighbor，WNN）方法，这是一个分析框架，用于集成细胞内测量的多种数据类型，并获得细胞状态的联合定义。该方法是基于非监督策略来学习细胞特定模态的“权重”，它反映每个模态的信息内容，并确定其在下游分析中的相对重要性。我们证明，WNN分析大大提高了我们定义多种生物数据类型中的细胞状态的能力。我们利用这种方法，基于包含211,000人外周血单核细胞（PBMC）的CITE-seq数据集生成多模式“图谱”，具有可扩展228个抗体的大细胞表面蛋白标记panel。利用这个数据集来识别和验证人类淋巴细胞中的异质细胞状态，并探索人类免疫系统对疫苗接种和SARS-CoV-2感染的反应。WNN在开源R工具包Seurat的更新版本中实现，代表了对单细胞数据进行综合多模态分析的广泛适用的策略。

文章使用脐带血单核细胞的CITE-seq数据和10个免疫标记共检测8617个细胞来进行算法的构建。要整合分析这两种状态的数据，要求分析方法满足以下条件：第一，robust，适应不同模态的数据；第二，能够进行多模态下游分析；第三，多模态比单模态下，性能能强。基于这个数据和要求构建了WNN的算法。如下图所示，在分析转录组时，CD8+和CD4+ T细胞部分混合在一起，但在蛋白数据中清晰分离。相比之下，传统的树突状细胞（cDCs），以及罕见的红系祖细胞和小鼠类3T3对照，在分析RNA时形成不同的簇，但根据表面蛋白丰度显示存在与其他类型的细胞混合。对每个细胞，首先计算每个模态k=20个最近邻的集合，接下来分别对蛋白近邻的分子和RNA近邻的分子表达量求平均值，并将平均值与原始值进行比较。结果显示基于蛋白knn的预测比基于RNA knn的预测更准确。然后利用预测的相对准确性来计算RNA和蛋白质的权重，从而衡量每个细胞中的相对信息。

WNN工作流中，关键的步骤如下：1.获得各模态预测和跨模态预测；2.基于细胞特定带宽核（cell-specific bandwidth kernel）将这些预测转化为预测亲和力；3.使用softmax变换计算模态权重。RNA和蛋白质模态权重是非负的，对每个细胞都是唯一的，总和为1。

最后一步整合并创建一个加权最近邻图（WNN图），基于标准化后的RNA和蛋白质的加权平均值，计算一组新的knn的细胞。计算公式如下图：

验证数据集1：CITE-seq和25中抗体，共检测30672个细胞

该数据集的结果表明WNN的整合大大提高了对细胞状态的注释，相较于单一模态更加精细化，更加完善，例如T细胞组，在scRNA分析中基本被掩盖，但是却有较高的蛋白模态的权重。验证WNN的稳定性时，高斯噪音比重增加时会降低蛋白模态在数据分析中的比重。

验证数据集2：10xGenomic PBMC细胞的ATAC和转录组数据，共检测11351个细胞

该数据集结果表明，模态组合展现了更优秀的免疫亚群的分类，其中ATAC-seq数据更能分离初始CD8 +及CD4 + T细胞状态由于可靠的检测细胞特定类型开放的染色质区域。该算法能够更敏感和强劲捕获异质性，可灵活地应用于多种数据类型，进行综合多模态分析。

验证数据集3：ASAP-seq HumanPBMC细胞的ATAC数据和227个蛋白，共检测4725个细胞

验证数据集4：SHARE-seq 小鼠的皮肤细胞的ATAC数据和转录组数据，共检测34774个细胞

以上两个数据集同样证明了WNN优秀的整合分析的能力，更加的精细化。

文章应用这个分析方法研究了多个主题方向，其中之一就是人类外周血单核细胞的多模态图谱。利用CITE-seq技术以及优化的抗体panel和整合的WNN分析策略，生成人类PBMC的多模态图谱。从8名参与艾滋病毒疫苗试验的志愿者中获得了PBMC样本，年龄跨度20-49岁（中位年龄36.5岁）。每个受试者在三个时间点采集PBMCs：注射HIV疫苗前（第0天）、第3天和第7天。整个数据集由24个样本组成，并使用“Cell hash”来最小化技术批次效应。对于每个样本，我们使用10X Chromium 3 '（使用228 TotalSeq A抗体）对细胞进行分析，总共代表了161,764个细胞（平均8,003个RNA分子/细胞，5,251个ADT/细胞）。并且还使用ECCITE-seq对所有样本中共49,147个细胞进行了分析，该技术可使用10X 5 '技术对表面蛋白进行。虽然后一组实验包含了54种抗体，其中包括实验室偶联抗体和TotalSeq-C试剂，反映了在实验时商业偶联的可用性，但我们也能够对这些细胞进行免疫库图谱分析。经过NovaSeq测序、严格的质量控制和双重过滤（补充方法），我们最终的数据集包含210,911个细胞，并允许我们分析静息（未接种）和激活（接种后）免疫系统的细胞异质性。

该WNN分析中鉴定了57个类群，包括所有主要和次要的免疫细胞类型，并揭示了细胞的多样性，特别是在淋巴细胞中。除了罕见的细胞类型外，每一类群的细胞都来自全部24个样本。我们的聚类可以分为几个大类别，包括CD4 + T细胞（12类），CD8 + T细胞（12类），非传统的T细胞（7类），NK细胞（6类），B细胞，浆细胞和plasmablasts（8类），树突细胞和单核细胞（8类），和罕见的集群造血祖细胞、血小板、红细胞和循环先天淋巴细胞（ILC）。为了更好的解释聚类结果，文章为将细胞进行三个粒度越来越大的注释（级别1，8个类别；第2级，30个类别；3级，57个类别）。虽然在T细胞亚群有较大程度的异质性，我们的分析明确确定异构子集的髓细胞与最近的高分辨率scRNA-seq完全整合分析排序的数量，包括极其罕见的人群（0.02%）定义的树突状细胞表达 AXL 和SIGLEC6。

总之，WNN算法的分析有助于揭示细胞的亚种群差异。虽然我们目前对WNN分析的实现侧重于对两种模式的分析，但随着这些技术的成熟，该框架可以很容易地扩展到处理任意数量的多模态数据。因此，其为综合多模态分析提供了一种途径，可以超越细胞的局部和转录聚焦的观点，并对细胞行为、身份和功能进行统一定义。

[1] Hao Y , Hao S , Andersen-Nissen E , et al. Integrated analysis of multimodal single-cell data[J]. 2021.
[2]
[3]

2019-10-22 R语言Seurat包下游分析-1

下游分析

cellranger count 计算的结果只能作为错略观测的结果，如果需要进一步分析聚类细胞，还需要进行下游分析，这里使用官方推荐 R 包（Seurat 3.0）

流程参考官方外周血分析标准流程（ ）

Rstudio操作过程：

## 安装seurat

es('Seurat')

## 载入seurat包

library(dplyr)

library(Seurat)

## 读入pbmc数据（文件夹路径不能包含中文，注意“/“的方向不能错误，这里读取的是10x处理的文件，也可以处理其它矩阵文件，具体怎样操作现在还不知道，文件夹中的3个文件分别是：，，，文件的名字不能错，否则读取不到）

<- Read10X( = "D:/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/")

## 查看稀疏矩阵的维度,即基因数和细胞数

dim()

[1:10,1:6]

## 创建Seurat对象与数据过滤，除去一些质量差的细胞（这里读取的是单细胞 count 结果中的矩阵目录；在对象生成的过程中，做了初步的过滤；留下所有在>=3 个细胞中表达的基因 = 3；留下所有检测到>=200 个基因的细胞 = 200。)

pbmc <- CreateSeuratObject(counts = , project = "pbmc3k", = 3, es = 200)

pbmc

##计算每个细胞的线粒体基因转录本数的百分比（%）,使用[[ ]] 操作符存放到metadata中，mit-开头的为线粒体基因

pbmc[[""]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

##展示基因及线粒体百分比（这里将其进行标记并统计其分布频率，"nFeature_RNA"为基因数，"nCount_RNA"为细胞数，""为线粒体占比）

VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", ""), ncol = 3)

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "")

plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

CombinePlots(plots = list(plot1, plot2))

## 过滤细胞：根据上面小提琴图中基因数"nFeature_RNA"和线粒体数""，分别设置过滤参数，这里基因数 200-2500，线粒体百分比为小于 5%，保留gene数大于200小于2500的细胞；目的是去掉空GEMs和1个GEMs包含2个以上细胞的数据；而保留线粒体基因的转录本数低于5%的细胞,为了过滤掉死细胞等低质量的细胞数据。

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & < 5)

## 表达量数据标准化,LogNormalize的算法：A = log( 1 + ( UMIA ÷ UMITotal ) × 10000

pbmc <- NormalizeData(pbmc, = "LogNormalize", = 10000)

#pbmc <- NormalizeData(pbmc) 或者用默认的

## 鉴定表达高变基因(2000个）,用于下游分析,如PCA；

pbmc <- FindVariableFeatures(pbmc, = "vst", nfeatures = 2000)

## 提取表达量变化最高的10个基因；

top10 <- head(VariableFeatures(pbmc), 10)

top10

plot1 <- VariableFeaturePlot(pbmc)

plot2 <- LabelPoints(plot = plot1, points = top10)

CombinePlots(plots = list(plot1, plot2))

plot1<-VariableFeaturePlot(object=pbmc)

plot2<-LabelPoints(plot=plot1,points=top10,repel=TRUE)

CombinePlots(plots=list(plot1,plot2))

## PCA分析：

# PCA分析数据准备,使用ScaleData()进行数据归一化；默认只是标准化高变基因（2000个）,速度更快,不影响PCA和分群,但影响热图的绘制。

#pbmc <- ScaleData(pbmc,s ="")

## 而对所有基因进行标准化的方法如下：

<- rownames(pbmc)

pbmc <- ScaleData(pbmc, features = )

pbmc <- ScaleData(pbmc, s = "")

## 线性降维（PCA）,默认用高变基因集,但也可通过features参数自己指定；

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

## 展示 pca 结果（最简单的方法）

DimPlot(object=pbmc,reduction="pca")

## 检查PCA分群结果, 这里只展示前5个PC,每个PC只显示5个基因；

print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)

##PC_ 1 

##Positive:  RPS27, MALAT1, RPS6, RPS12, RPL13 

##Negative:  CSTA, FCN1, CST3, LYZ, LGALS2 

##PC_ 2 

##Positive:  NKG7, GZMA, CST7, KLRD1, CCL5 

##Negative:  RPL34, RPL32, RPL13, RPL39, LTB 

##PC_ 3 

##Positive:  MS4A1, CD79A, BANK1, IGHD, CD79B 

##Negative:  IL7R, RPL34, S100A12, VCAN, AIF1 

##PC_ 4 

##Positive:  RPS18, RPL39, RPS27, MALAT1, RPS8 

##Negative:  PPBP, PF4, GNG11, SDPR, TUBB1 

##PC_ 5 

##Positive:  PLD4, FCER1A, LILRA4, SERPINF1, LRRC26 

##Negative:  MS4A1, CD79A, LINC00926, IGHD, FCER2 

## 展示主成分基因分值

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

## 绘制pca散点图

DimPlot(pbmc, reduction = "pca")

## 画第1个或15个主成分的热图；

DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)

DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)

## 确定数据集的分群个数

# 鉴定数据集的可用维度，方法1：Jackstraw置换检验算法；重复取样（原数据的1%）,重跑PCA,鉴定p-value较小的PC；计算‘null distribution’(即零假设成立时)时的基因scores。虚线以上的为可用维度，也可以调整 dims 参数，画出所有 pca 查看。

#pbmc <- JackStraw(pbmc, ate = 100)

#pbmc <- ScoreJackStraw(pbmc, dims = 1:20)

#JackStrawPlot(pbmc, dims = 1:15)

# 方法2：肘部图（碎石图）,基于每个主成分对方差解释率的排名。

ElbowPlot(pbmc)

## 细胞聚类：分群个数这里选择10,建议尝试选择多个主成分个数做下游分析,对整体影响不大；在选择此参数时,建议选择偏高的数字（为了获取更多的稀有分群,“宁滥勿缺”）；有些亚群很罕见,如果没有先验知识,很难将这种大小的数据集与背景噪声区分开来。

## 非线性降维（UMAP/tSNE)基于PCA空间中的欧氏距离计算nearest neighbor graph,优化任意两个细胞间的距离权重（输入上一步得到的PC维数） 。

pbmc <- FindNeighbors(pbmc, dims = 1:10)

## 接着优化模型,resolution参数决定下游聚类分析得到的分群数,对于3K左右的细胞,设为0.4-1.2 能得到较好的结果(官方说明)；如果数据量增大,该参数也应该适当增大。

pbmc <- FindClusters(pbmc, resolution = 0.5)

## 使用Idents（）函数可查看不同细胞的分群；

head(Idents(pbmc), 5)

## 结果：AAACCTGAGGTGCTAG    AAACCTGCAGGTCCAC    AAACCTGCATGGAATA AAACCTGCATGGTAGG      AAACCTGCATTGGCGC 

               1                3                0               10                2 

Levels: 0 1 2 3 4 5 6 7 8 9 10 11

## Seurat提供了几种非线性降维的方法进行数据可视化（在低维空间把相似的细胞聚在一起）,比如UMAP和t-SNE,运行UMAP需要先安装'umap-learn'包，这里不做介绍，两种方法都可以使用，但不要混用，如果混用，后面的结算结果会将先前的聚类覆盖掉，只能保留一个。

## 这里采用基于TSNE的聚类方法。

pbmc <- RunTSNE(pbmc, dims = 1:10)

## 用DimPlot()函数绘制散点图,reduction = "tsne",指定绘制类型；如果不指定,默认先从搜索 umap,然后 tsne, 再然后 pca；也可以直接使用这3个函数PCAPlot()、TSNEPlot()、UMAPPlot()； cols,分别调整分组颜色和点的大小；

DimPlot(pbmc,reduction = "tsne",label = TRUE, = 1.5)

## 这里采用基于图论的聚类方法

pbmc<-RunUMAP(object=pbmc,dims=1:10)

DimPlot(object=pbmc,reduction="umap")

## 细胞周期归类

pbmc<- CellCycleScoring(object = pbmc, es = $, es = $)

head(x = )

DimPlot(pbmc,reduction = "tsne",label = TRUE,="Phase", = 1.5)

## 存储结果

saveRDS(pbmc, file = "D:/")

save(pbmc,file="D:/")

## 寻找cluster 1的marker

s <- FindMarkers(pbmc, ident.1 = 1, = 0.25)

head(s, n = 5)

## 结果:      p_val             avg_logFC        pct.1       pct.2        p_val_adj

MT-CO1  0.000000e+00    -0.6977083      0.985       0.996       0.000000e+00

RPS27  2.182766e-282     0.3076454       1.000       0.999        3.480202e-278

MT-CO3 2.146399e-274    -0.4866429      0.995       0.997       3.422218e-270

DUSP1  2.080878e-247    -1.7621662       0.376       0.745       3.317752e-243

RPL34  8.647733e-244     0.3367755        1.000       0.997       1.378795e-239

##寻找每一cluster的marker

s <- FindAllMarkers(pbmc, = TRUE, = 0.25, old = 0.25)

s %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)

# A tibble: 24 x 7

# Groups:   cluster [12]

       p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene 

       <dbl>     <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>

 1 2.29e-123     0.636 0.344 0.097 3.65e-119 0       CD8B 

 2 7.62e-113     0.487 0.632 0.305 1.22e-108 0       LEF1 

 3 2.04e- 74     0.483 0.562 0.328 3.25e- 70 1       LEF1 

 4 1.39e- 61     0.462 0.598 0.39  2.22e- 57 1       ITM2A

 5 0.            2.69  0.972 0.483 0.        2       GNLY 

 6 0.            2.40  0.964 0.164 0.        2       GZMB 

 7 1.31e-121     0.768 0.913 0.671 2.09e-117 3       JUNB 

 8 2.06e- 94     0.946 0.426 0.155 3.28e- 90 3       RGS1 

 9 2.05e-255     1.57  0.586 0.09  3.27e-251 4       GZMK 

10 2.94e-140     1.57  0.69  0.253 4.68e-136 4       KLRB1

# ... with 14 more rows

## 存储marker

(s,file="D:/")

## 各种绘图

## 绘制Marker 基因的tsne图

FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"),cols = c("gray", "red"))

## 绘制Marker 基因的小提琴图

VlnPlot(pbmc, features = c("MS4A1", "CD79A"))

VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)

## 绘制分cluster的热图

top10 <- s %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)

DoHeatmap(pbmc, features = top10$gene) + NoLegend()

剩下的便是寻找基因 marker 并对细胞类型进行注释(见下回分解)