单细胞转录组高级分析八：整合V(D)J数据

单细胞天地 2022-06-07

The following article is from 生信会客厅 Author Kinesin

上期专题我们介绍了单细胞转录组数据的基础分析，然而那些分析只是揭开了组织异质性的面纱，还有更多的生命奥秘隐藏在数据中等待我们发掘。本专题将介绍一些单细胞转录组的高级分析内容：多样本批次校正、转录因子分析、细胞通讯分析、基因集变异分析和更全面的基因集富集分析。不足之处请大家批评指正，欢迎添加Kinesin微信交流探讨!

10X V(D)J测序简介

研究免疫的朋友一定对V(D)J基因重排表达TCR/BCR的过程了然于心，作为免疫学外行的我就不在此赘述了，我们把重点放在10X V(D)J测序的原理和优势上。

10X V(D)J测序原理

10X V(D)J测序与其转录组测序建库流程大部分相同，都要经历微流控构建单细胞油包水反应体系，mRNA逆转录为cDNA并扩增的过程。与常规10X单细胞转录组测序把barcode和UMI序列放在3'端不同，V(D)J测序要把barcode和UMI序列放在5'端。

逆转录的cDNA经过扩增后可以一分为二，一部分做5'端scRNA测序，另一部分用富集引物扩增V(D)J序列测序。

人的V(D)J序列长几百至一千多个碱基，二代测序怎么得到全长序列呢？大家先看看下面的V(D)J测序文库构建原理图：

请注意上图中的红色文字，扩增得到的全长V(D)J序列会被酶切为长短不一的片段。因为每条V(D)J序列都有很多个PCR扩增的拷贝，并且这些拷贝有相同的UMI标签，所以每条V(D)J序列都会得到一组测序reads（通过UMI确定同一来源），如下图所示：

Read1基本固定检测V(D)J序列5‘端150bp的序列，reads2则是随机覆盖V(D)J序列的各个位置（随机打断位置150bp的序列）。只要有足够的数据量，最终测序数据将覆盖全部的V(D)J序列。

10X V(D)J测序优势

scRepertoire简介

scRepertoire由华盛顿大学的Nick Borcherding博士开发，是一款针对10X V(D)J数据的免疫组库分析工具。它直接导入10X cellranger的输出结果，具有丰富的分析项目和美观的可视化结果。

scRepertoire目前发布在github，正式版只支持R4.0，R3.6版本只能使用开发版。开发版相比正式版应该有一些功能上的差异，我运行官方教程时发现有些函数开发版不支持。

library(devtools)#正式版，需要4.0版R语言支持install_github("ncborcherding/scRepertoire")#开发版，适用4.0版以下的R语言install_github("ncborcherding/scRepertoire@dev")

V(D)J数据分析

导入文件说明

scRepertoire需要导入后缀为contig_annotations.csv的CellRanger输出文件，要求barcode没有前缀且以列表的方式传递给scRepertoire。本教程使用scRepertoire包自带的VDJ数据，可通过data("contig_list")加载。

library(Seurat)library(scRepertoire)library(tidyverse)library(patchwork)rm(list=ls())dir.create("VDJ")#加载scRepertoire自带的contig_annotations文件data("contig_list")

查看contig_list[[1]]的前5行，可以发现这是标准的contig注释数据

cellranger输出的文件中，没有细胞barcode与clonotype直接对应的数据，因此不能直接整合到seurat中。scRepertoire分析的第一步是把上面的文件转换为barcode与clonotype对应的数据，它通过combineTCR函数实现。

combined <- combineTCR(contig_list, samples = c("PY", "PY", "PX", "PX", "PZ","PZ"), ID = c("P", "T", "P", "T", "P", "T"), cells ="T-AB")#参数说明：

#samples：用于指定样本名称的参数，因为contig_list包含了6个样本的VDJ数据，所以这里传递了一个6个字符元素的向量；

#ID：用于指定样本名称之外的其他分类信息，例如分组信息；#cells：指定VDJ类型，"T-AB"代表TCR的α和β链配对，"T-GD"代表γ和δ配对；

#还有removeNA, removeMulti, filterMulti都选择了默认值，有兴趣了解的可以查看帮助文档

输出结果combined依然是一个列表，我们看看combined[[1]]的前5行

可以看到barcode都加上了由samples和ID组成的前缀，并且把每个细胞配对的两条链放在一行，之后的分析都基于配对信息定义的clonotype。下面分析所用的函数，一般都有两个通用参数：

exportTable：默认赋值F，函数分析结果为图形；改为T之后函数分析结果输出为表格。
cloneCall：可选"aa", "nt", "gene", "gene+nt"，代表用CDR3氨基酸序列、CDR3核苷酸序列、VDJ基因还是VDJ基因+CDR3核苷酸序列中的哪一种来定义克隆型。以下分析主要用"aa"定义clonotype，大家可以根据自己的需要选择其他类型。

##展示每个样本的克隆型数量p1 <- quantContig(combined, cloneCall="gene+nt", scale = F)p2 <- quantContig(combined, cloneCall="gene+nt", scale = T)plotc = p1 + p2ggsave('VDJ/quantContig.png', plotc, width = 8, height = 4)#设置exportTable = T，则输出表格而非图形quantContig(combined, cloneCall="gene+nt", scale = T, exportTable = T # contigs values total scaled# 1 2692 PY_P 3208 83.91521# 2 1513 PY_T 3119 48.50914# 3 823 PX_P 1068 77.05993# 4 928 PX_T 1678 55.30393# 5 1147 PZ_P 1434 79.98605# 6 764 PZ_T 2768 27.60116

左图展示的是每个样本克隆型总数，右图展示的是克隆型总数/克隆种数

##CDR3序列的长度分布，“aa”代表统计氨基酸序列长度p1 <- lengthContig(combined, cloneCall="aa", chains = "combined") p2 <- lengthContig(combined, cloneCall="aa", chains = "single")plotc = p1 + p2ggsave('VDJ/lengthContig.png', plotc, width = 8, height = 4)

左图展示的是α和β链一起统计，右图展示的是α和β链分开统计

##对比两个样本的克隆型p1 = compareClonotypes(combined, numbers = 10, samples = c("PX_P", "PX_T"), cloneCall="aa", graph = "alluvial")ggsave('VDJ/compareClonotypes.png', p1, width = 8, height = 4)

展示排名前10的克隆型在两个样本间的对比情况

##Clonal Space Homeostasis，这个不知道怎么翻译clonalHomeostasis(combined, cloneCall = "aa")#此分析将克隆型按其相对丰度分为rare, small, medium, large, hyperexpanded#5大类，并统计各个类别的占比

##克隆型分类占比，与上一个分析相似，只是分类方法调整了clonalProportion(combined, cloneCall = "aa")

##Overlap Analysis，分析样本相似性clonalOverlap(combined, cloneCall = "aa", method = "morisita")

##克隆多样性指数clonalDiversity(combined, cloneCall = "aa", group = "samples")

V(D)J与scRNA整合分析

前面介绍的分析内容只能依据样本或分组信息（可通过ID参数指定）开展，其实我们还可以通过barcode将VDJ数据与scRNA数据联系起来，这样就可以将clonotype显示在降维图上，也可以基于cluster展示clonotype的情况。

官方提供的scRNA数据需要下载：https://drive.google.com/open?id=1np-EzG7U9W_Fz_SchBrsAhtqE3_rB_H9

下载之后会得到一个seurat2.rda的文件，这是一个降维聚类后的seurat对象，请保存在分析相关目录下。

seurat <- get(load("Data/VDJ/seurat2.rda"))#查看UMAP图DimPlot(seurat, label = T) + NoLegend()

#查看原始seurat的meta.data有哪些内容names(seurat@meta.data)# [1] "nCount_RNA" "nFeature_RNA" "integrated_snn_res.0.5" "seurat_clusters"# [5] "Patient" "Type" "RawBarcode" ##将VDJ数据添加到seurat对象的meta.data中seurat <- combineExpression(combined, seurat, cloneCall="aa", groupBy = "sample")#查看添加VDJ数据后seurat的meta.data有哪些内容names(seurat@meta.data) [1] "nCount_RNA" "nFeature_RNA" "integrated_snn_res.0.5" [4] "seurat_clusters" "Patient" "Type" [7] "RawBarcode" "barcode" "CTgene" [10] "CTnt" "CTaa" "CTstrict" [13] "Frequency" "cloneType"

#可以发现导入了双链的相关信息，增加了克隆型的频率信息（Frequency），并且将频率划分了层次（cloneType）

##UMAP图展示cloneType的分布seurat$cloneType <- factor(seurat$cloneType, levels = c("Hyperexpanded (100 < X <= 500)", "Large (20 < X <= 100)", "Medium (5 < X <= 20)", "Small (1 < X <= 5)", "Single (0 < X <= 1)", NA))colorblind_vector <- colorRampPalette(c("#FF4B20", "#FFB433", "#C6FDEC", "#7AC5FF", "#0348A6"))DimPlot(seurat, group.by = "cloneType") + scale_color_manual(values = colorblind_vector(5), na.value="grey")

##UMAP图展示特定克隆型的分布seurat <- highlightClonotypes(seurat, cloneCall= "aa", sequence = c("CAVNGGSQGNLIF_CSAEREDTDTQYF", "NA_CATSATLRVVAEKLFF"))DimPlot(seurat, group.by = "highlight")

##桑基图展示特定克隆型的来源alluvialClonotypes(seurat, cloneCall = "aa", y.axes = c("Patient", "cluster", "Type"), color = "CAVNGGSQGNLIF_CSAEREDTDTQYF") + scale_fill_manual(values = c("grey", colorblind_vector(1)))

##桑基图展示克隆型在病例-cluster-组织类型之间的关系alluvialClonotypes(seurat, cloneCall = "aa", y.axes = c("Patient", "cluster", "Type"), color = "cluster")

##按cluster分析克隆型combined2 <- expression2List(seurat, group = "cluster")p1 = clonalDiversity(combined2, cloneCall = "aa") + ggtitle("clonalDiversity")p2 = clonalHomeostasis(combined2, cloneCall = "aa") + ggtitle("clonalHomeostasis")p3 = clonalProportion(combined2, cloneCall = "aa") + ggtitle("clonalProportion")p4 = clonalOverlap(combined2, cloneCall="aa", method="overlap") + ggtitle("clonalOverlap") plotc = (p2|p3)/(p1|p4)ggsave('VDJ/clonal.png', plotc, width = 10, height = 8)