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

上期专题我们介绍了单细胞转录组数据的基础分析，然而那些分析只是揭开了组织异质性的面纱，还有更多的生命奥秘隐藏在数据中等待我们发掘。本专题将介绍一些单细胞转录组的高级分析内容：多样本批次校正、转录因子分析、细胞通讯分析、基因集变异分析和更全面的基因集富集分析。不足之处请大家批评指正，欢迎添加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 + p2
ggsave('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 + p2
ggsave('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")

VDJ与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([e-mail])
# [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([e-mail])
 [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)

后记

我是这个月才知道scRepertoire包的，感觉它的可视化效果非常好，因此第一时间与大家分享。如果您也觉得scRepertoire不错，请积极转发此文共同推广，欢迎大家与我交流scRepertoire的使用经验！

本篇是《单细胞转录组高级分析》专题的最后一篇，这两天我会把代码和数据整理好上传百度云，需要的朋友添加Kinesin微信索取。

获取帮助

本教程的目的在于把单细胞分析流程串起来，适合有一定R语言基础的朋友参考。如果您学习本教程有一定困难，可以分享此篇文章到朋友圈，截图微信发给Kinesin（文末二维码），我会抽时间组群直播讲解单细胞数据分析的全过程。本专题所用的软件、数据、代码脚本和中间结果，我打包放在了百度云上，需要的朋友添加Kinesin微信索取。