scRepertoire||单细胞免疫组库分析：R语言应用(一)
scRepertoire||单细胞免疫组库分析：R语言应用(一)
分享是一种态度
作者 | 周运来
男，
一个长大了才会遇到的帅哥，
稳健，潇洒，大方，靠谱。
一段生信缘，一棵技能树，
一枚大型测序工厂的螺丝钉，
一个随机森林中提灯觅食的津门旅客。
前情回顾
10× Genomics单细胞免疫组库VDJ分析必知必会(https://www.jianshu.com/p/db4831091a5c) 免疫组库数据分析||immunarch教程：快速开始(https://www.jianshu.com/p/9d7711879bf5) 免疫组库数据分析||immunarch教程：克隆型分析(https://www.jianshu.com/p/287f890d7ef4) 免疫组库数据分析||immunarch教程：探索性数据分析(https://www.jianshu.com/p/dd4fcfb63627) 免疫组库数据分析||immunarch教程：载入10X数据(https://www.jianshu.com/p/7379d0a809a8) 免疫组库数据分析||immunarch教程：GeneUsage分析(https://www.jianshu.com/p/0dbdd6733b34) 免疫组库数据分析||immunarch教程：Diversity 分析(https://www.jianshu.com/p/8b846094c092) 免疫组库数据分析||immunarch教程：Clonotype tracking(https://www.jianshu.com/p/79ee2c5871a7) 免疫组库数据分析||immunarch教程：Clonotypes annotation(https://www.jianshu.com/p/effc2ad05f47) 免疫组库数据分析||immunarch教程：Kmer 与 Motif 分析(https://www.jianshu.com/p/f2b7d0153432)
在immunarch 教程中我们比较细致地分析了免疫组库数据，但是没有把免疫组库数据和单细胞转录组数据结合在一起。今天介绍一下scRepertoire，可以弥补这个缺憾。其实在immunarch 中也可以根据barcode信息将两者整合到一起，只是需要代码处理。这里仅根据其文章的配图看看scRepertoire的主要功能。图文来自：https://f1000research.com/articles/9-47/v1
单细胞测序是免疫学和肿瘤学领域的新兴技术，它允许研究人员将RNA定量和其他模式数据结合起来，如在单个细胞水平上的免疫细胞受体分析。许多工作流和软件包已经被创建来处理和分析单细胞转录数据。这些软件包允许用户将基于单细胞的实验中产生的海量数据提取新颖的见解。而单细胞免疫组库目前还缺乏成熟的数据分析软件。通过构建scRepertoire，用户可以很容易地结合mRNA和免疫组分析，从而处理来自10x Genomics的T细胞受体库(TCR)和免疫球蛋白(Ig)工作流程的免疫分析数据，并与流行的Seurat R包兼容。scRepertoire 包和处理过的数据都是开源的，可以在GitHub上获得，那里提供了关于这个包功能的详细教程。
https://ncborcherding.github.io/vignettes/vignette.html
scRepertoire是在R v3.5.1中构建和测试的（我在4.0.2的R上也安装运行了，没事）。scRepertoire分析的灵感来自于没有代码派生的bulk 免疫分析工具 tcR(v2.2.4) R包。克隆型可以通过结合免疫位点基因(一种更敏感的方法)或互补决定区3 (CDR3)的核苷酸/氨基酸序列来命名。除了R中的基本函数之外，还使用dplyr (v0.8.3)和rempe2 (v1.4.3) R包执行数据处理。可视化是使用ggplot2 (v3.2.1)和ggalluvial(v0.11.1)R包生成的，颜色来自于使用colorRamps (v2.3)和RColorBrewer (v1.1.2)R包。多样性度量使用vegan(v2.5-6) R包计算。函数的可视化输出存储为几何或统计ggplot分层对象，允许用户方便地修改表示。
(A)根据患者和标本类型的TCRs序列的总数量，使用quantContig函数来计算独特的克隆型(外周,P; 肿瘤,T)。
(B)利用abundanceContig 函数按样品和类型划分的克隆型总丰度。
(C)利用外周血和肿瘤样本密度比较克隆型的相对丰度。
(D)利用lengthContig函数对样本进行CDR3核苷酸长度分析。该曲线的双峰性质是一个函数调用克隆型细胞的一个和两个免疫受体测序。
(A)克隆稳态空间在所有6个样本中使用基因和CDR3 AA序列进行克隆型识别。
(B)利用基因和CDR3 AA序列进行克隆型识别，在所有6个样本中特定克隆型所占的相对比例。
(C)所有6个样品中克隆型的Morisita overlap。
(D)利用Shannon、Inverse Simpson、Chao和基于丰度的覆盖率估计(abundance-based coverage estimator ，ACE)指数按样本类型基于克隆类型的多样性度量。
(A) UMAP将ccRCC T细胞(n= 12911)投射成12个不同的簇。
(B) UMAP投影，突出外周血(红色)和肿瘤(蓝色)人群，并伴随每个簇的相对比例组成，分别以外周血和肿瘤细胞总数比例。
(C)使用combineSeurat函数，将单个细胞按克隆类型的数量分组，然后可以在UMAP投影上叠加显示。
(D) combineSeurat计算克隆型的频率，可用于检查亚群组成，如箱线图所示。
(E)将克隆类型信息与Seurat对象进行梳理后，利用highlightClonotypes可以利用序列信息具体突出感兴趣的单个克隆类型
(F)多类间克隆型的相互作用可以利用alluvialGraph函数可视化。
这个最后这一个fig可以说是scRepertoire的特色了，把单细胞免疫组库信息mapping到单细胞RNA降维空间上。
scRepertoire: An R-based toolkit for single-cell immune receptor analysis
https://ncborcherding.github.io/vignettes/vignette.html
好了，到了你们家运来哥哥展示搬砖工基本功的时候了。
scRepertoire旨在获取来自10x Genomics Cell Ranger管道的过滤contig输出(filtered_contig_annotations.csv)，处理这些数据以基于两个TCR或Ig链分配克隆型，并分析克隆型动力学。后者可分为
- 1)单独克隆型分析功能，如独特克隆型或克隆空间定量;
- 2)与mRNA表达数据的相互作用，使用Seurat、singlecelperperor Monocle 3包。
加载R包和数据
一如往常我们先把环境配置好：
# 没有安装的话先安装
suppressMessages(library(scRepertoire))
suppressMessages(library(Seurat))
scRepertoire提供了一组示例数据，这些数据来自于三个肾透明细胞癌患者的T细胞，目的是为了展示R包的功能。更多关于数据集的信息可以在预印1和预印2中找到（preprint 1(https://www.biorxiv.org/content/10.1101/478628v1.abstract) and preprint 2(https://www.biorxiv.org/content/10.1101/824482v1.abstract) )。这些样本由成对的外周血和肿瘤浸润组成，有效地为T细胞受体(TCR)富集创造了6个不同的runs。我们可以使用head函数预览列表中的元素，并查看第一个contig注释。：这里注意条码贴上PX_P_ # # # # # # # # # # # # # -这是指患者X (PX)和外周血(P)。
如果你是filtered_contig_annotation.csv文件加载到R环境创建列表,您还需要调用stringsAsFactors 为 FALSE ,这将防止分类变量的转换为内置的因素和必要的一些scRepertoire的函数。代码应该看起来像这样:
# 这个数据你没有也没关系。示例用不到啊
csv1 <- read.csv("F:\\VDJ\\filtered_contig_annotations.csv", stringsAsFactors = F)
mycsv1 = list(ex=csv1)
对象“contig_list”是从10x Genomics Cell Ranger中的6个filtered_contig_annotation.csv文件中创建的。该对象是用contig_list <- list(csv1, csv2，…)创建的。
data("contig_list")
head(contig_list[[1]])
barcode is_cell contig_id high_confidence
1 AAACCTGAGAGCTGGT TRUE AAACCTGAGAGCTGGT-1_contig_1 TRUE
2 AAACCTGAGAGCTGGT TRUE AAACCTGAGAGCTGGT-1_contig_2 TRUE
3 AAACCTGAGCATCATC TRUE AAACCTGAGCATCATC-1_contig_1 TRUE
4 AAACCTGAGCATCATC TRUE AAACCTGAGCATCATC-1_contig_2 TRUE
5 AAACCTGAGCATCATC TRUE AAACCTGAGCATCATC-1_contig_5 TRUE
6 AAACCTGAGTGGTCCC TRUE AAACCTGAGTGGTCCC-1_contig_1 TRUE
length chain v_gene d_gene j_gene c_gene full_length productive
1 705 TRB TRBV20-1 TRBD1 TRBJ1-5 TRBC1 TRUE TRUE
2 502 TRB None None TRBJ1-5 TRBC1 FALSE None
3 693 TRB TRBV5-1 TRBD2 TRBJ2-2 TRBC2 TRUE TRUE
4 567 TRA TRAV12-1 None TRAJ37 TRAC TRUE TRUE
5 361 TRB None None TRBJ1-5 TRBC1 FALSE None
6 593 TRB TRBV7-9 TRBD1 TRBJ2-5 TRBC2 TRUE TRUE
cdr3 cdr3_nt reads
1 CSASMGPVVSNQPQHF TGCAGTGCTAGCATGGGACCGGTAGTGAGCAATCAGCCCCAGCATTTT 16718
2 None None 6706
3 CASSWSGAGDGELFF TGCGCCAGCAGCTGGTCAGGAGCGGGAGACGGGGAGCTGTTTTTT 26719
4 CVVNDEGSSNTGKLIF TGTGTGGTGAACGATGAAGGCTCTAGCAACACAGGCAAACTAATCTTT 18297
5 None None 882
6 CASSPSEGGRQETQYF TGTGCCAGCAGCCCCTCCGAAGGGGGGAGACAAGAGACCCAGTACTTC 11218
umis raw_clonotype_id raw_consensus_id
1 6 clonotype96 clonotype96_consensus_1
2 3 clonotype96 None
3 11 clonotype97 clonotype97_consensus_2
4 6 clonotype97 clonotype97_consensus_1
5 4 clonotype97 None
6 6 clonotype98 clonotype98_consensus_1
一些工作流中将有标准条形码标签额外的标签。在继续之前，我们将使用函数stripBarcode()来避免任何标签问题。重要的是，stripBarcode()用于删除来自其他管道的条形码上的前缀。所以你需要先看看有什么标签，以及标签的格式，同时建议对比10X的filtered_contig_annotation文件。
如果条形码为+ AAACGGGAGATGGCGT-1 + AAACGGGAGATGGCGT，则不需要stripBarcode功能。
for (i in seq_along(contig_list)) {
contig_list[[i]] <- stripBarcode(contig_list[[i]], column = 1, connector = "_", num_connects = 3)
}
合并多样本Contigs
CellRanger的输出是对TCRA和TCRB链的量化，下一步是通过细胞条形码创建一个带有TCR基因和CDR3序列的单一列表对象。这是使用combineTCR()执行的，其中的输入的是contig_list。还可以根据样本和身份信息重新标记条形码，以防止重复。如果是BCR该用什么呢？combineBCR。
cells + T-AB - T cells, alpha-beta TCR + T-GD - T cells, gamma-delta TCR
removeNA + TRUE -这是一个严格的过滤器，用于移除至少有一个NA值的细胞条码+ FALSE -包含和合并NA值为1的细胞的默认设置。
removeMulti + TRUE -这是一个严格的过滤器，可以移除任何超过2个免疫受体链的细胞条码+ FALSE -包含和合并带有> 2链的细胞的默认设置。
filterMulti + TRUE -用多个链分离细胞条形码中前2个表达的链+ FALSE -包含和合并细胞与> 2链的默认设置。
?combineTCR
combined <- combineTCR(contig_list, samples = c("PY", "PY", "PX", "PX", "PZ","PZ"), ID = c("P", "T", "P", "T", "P", "T"), cells ="T-AB")
head(combined$PY_P)
barcode sample ID TCR1 cdr3_aa1
1 PY_P_AAACCTGAGAGCTGGT PY P <NA> <NA>
2 PY_P_AAACCTGAGCATCATC PY P TRAV12-1.TRAJ37.TRAC CVVNDEGSSNTGKLIF
4 PY_P_AAACCTGAGTGGTCCC PY P <NA> <NA>
5 PY_P_AAACCTGCAAACGCGA PY P TRAV29DV5.TRAJ22.TRAC CAASGYGSARQLTF
7 PY_P_AAACCTGGTGAAAGAG PY P TRAV2.TRAJ42.TRAC CAAGGGGSQGNLIF
9 PY_P_AAACGGGCACGCTTTC PY P TRAV8-1.TRAJ13.TRAC CAVNAGTGGYQKVTF
cdr3_nt1
1 <NA>
2 TGTGTGGTGAACGATGAAGGCTCTAGCAACACAGGCAAACTAATCTTT
4 <NA>
5 TGTGCAGCAAGCGGTTACGGTTCTGCAAGGCAACTGACCTTT
7 TGTGCTGCGGGAGGGGGAGGAAGCCAAGGAAATCTCATCTTT
9 TGTGCCGTGAATGCGGGGACTGGGGGTTACCAGAAAGTTACCTTT
TCR2 cdr3_aa2
1 TRBV20-1.TRBJ1-5.TRBD1.TRBC1 CSASMGPVVSNQPQHF
2 TRBV5-1.TRBJ2-2.TRBD2.TRBC2 CASSWSGAGDGELFF
4 TRBV7-9.TRBJ2-5.TRBD1.TRBC2 CASSPSEGGRQETQYF
5 TRBV2.TRBJ1-6.TRBD1.TRBC1 CASRVQGNRGSPLHF
7 TRBV3-1.TRBJ2-5.TRBD1.TRBC2 CASSQNRVETQYF
9 <NA> <NA>
cdr3_nt2
1 TGCAGTGCTAGCATGGGACCGGTAGTGAGCAATCAGCCCCAGCATTTT
2 TGCGCCAGCAGCTGGTCAGGAGCGGGAGACGGGGAGCTGTTTTTT
4 TGTGCCAGCAGCCCCTCCGAAGGGGGGAGACAAGAGACCCAGTACTTC
5 TGTGCCAGCAGGGTACAGGGTAATAGGGGTTCACCCCTCCACTTT
7 TGTGCCAGCAGCCAAAACAGGGTAGAGACCCAGTACTTC
9 <NA>
CTgene
1 NA_TRBV20-1.TRBJ1-5.TRBD1.TRBC1
2 TRAV12-1.TRAJ37.TRAC_TRBV5-1.TRBJ2-2.TRBD2.TRBC2
4 NA_TRBV7-9.TRBJ2-5.TRBD1.TRBC2
5 TRAV29DV5.TRAJ22.TRAC_TRBV2.TRBJ1-6.TRBD1.TRBC1
7 TRAV2.TRAJ42.TRAC_TRBV3-1.TRBJ2-5.TRBD1.TRBC2
9 TRAV8-1.TRAJ13.TRAC_NA
CTnt
1 NA_TGCAGTGCTAGCATGGGACCGGTAGTGAGCAATCAGCCCCAGCATTTT
2 TGTGTGGTGAACGATGAAGGCTCTAGCAACACAGGCAAACTAATCTTT_TGCGCCAGCAGCTGGTCAGGAGCGGGAGACGGGGAGCTGTTTTTT
4 NA_TGTGCCAGCAGCCCCTCCGAAGGGGGGAGACAAGAGACCCAGTACTTC
5 TGTGCAGCAAGCGGTTACGGTTCTGCAAGGCAACTGACCTTT_TGTGCCAGCAGGGTACAGGGTAATAGGGGTTCACCCCTCCACTTT
7 TGTGCTGCGGGAGGGGGAGGAAGCCAAGGAAATCTCATCTTT_TGTGCCAGCAGCCAAAACAGGGTAGAGACCCAGTACTTC
9 TGTGCCGTGAATGCGGGGACTGGGGGTTACCAGAAAGTTACCTTT_NA
CTaa
1 NA_CSASMGPVVSNQPQHF
2 CVVNDEGSSNTGKLIF_CASSWSGAGDGELFF
4 NA_CASSPSEGGRQETQYF
5 CAASGYGSARQLTF_CASRVQGNRGSPLHF
7 CAAGGGGSQGNLIF_CASSQNRVETQYF
9 CAVNAGTGGYQKVTF_NA
CTstrict
1 NA_NA_TRBV20-1.TRBJ1-5.TRBD1.TRBC1_TGCAGTGCTAGCATGGGACCGGTAGTGAGCAATCAGCCCCAGCATTTT
2 TRAV12-1.TRAJ37.TRAC_TGTGTGGTGAACGATGAAGGCTCTAGCAACACAGGCAAACTAATCTTT_TRBV5-1.TRBJ2-2.TRBD2.TRBC2_TGCGCCAGCAGCTGGTCAGGAGCGGGAGACGGGGAGCTGTTTTTT
4 NA_NA_TRBV7-9.TRBJ2-5.TRBD1.TRBC2_TGTGCCAGCAGCCCCTCCGAAGGGGGGAGACAAGAGACCCAGTACTTC
5 TRAV29DV5.TRAJ22.TRAC_TGTGCAGCAAGCGGTTACGGTTCTGCAAGGCAACTGACCTTT_TRBV2.TRBJ1-6.TRBD1.TRBC1_TGTGCCAGCAGGGTACAGGGTAATAGGGGTTCACCCCTCCACTTT
7 TRAV2.TRAJ42.TRAC_TGTGCTGCGGGAGGGGGAGGAAGCCAAGGAAATCTCATCTTT_TRBV3-1.TRBJ2-5.TRBD1.TRBC2_TGTGCCAGCAGCCAAAACAGGGTAGAGACCCAGTACTTC
9 TRAV8-1.TRAJ13.TRAC_TGTGCCGTGAATGCGGGGACTGGGGGTTACCAGAAAGTTACCTTT_NA_NA
cellType
1 T-AB
2 T-AB
4 T-AB
5 T-AB
7 T-AB
9 T-AB
combineTCR()的输出将是一个contig数据帧列表，该列表将被简化为与单个细胞条码相关联的reads。它还会通过核苷酸序列(CTnt)、氨基酸序列(CTaa)、基因序列(CTgene)或核苷酸和基因序列的组合(CTstrict)将多个读码组合成克隆型调用。B细胞的类似函数，combineBCR()函数与2个主要注意事项类似:1)每个条码最多只能有2个序列，如果有更大的存在，选择reads值最高的2个。2)严格定义克隆型(clonotype, CTstrict)是基于v基因和>85%标准化的核苷酸序列汉明距离。汉明距离计算所有恢复的BCR序列，而不考虑运行runs。
其他功能
如果除了sample和ID之外还需要添加更多的变量呢?我们可以使用addVariable()函数来添加它们。我们需要的只是要添加的变量的名称和特定的字符或数字值(变量)。作为一个例子，这里我们添加样品被处理和测序的批次。
example <- addVariable(combined, name = "batch", variables = c("b1", "b1", "b2", "b2", "b2", "b2"))
example[[1]][1:5,ncol(example[[1]])] # This is showing the first 5 values of the new column added
同样，我们可以使用subsetContig()函数在combineTCR()之后删除特定的列表元素。为了进行子集化，我们需要确定要用于子集化的向量(名称)和要子集化的变量值(变量)。下面你可以看到我们从PX和PY中分离出4个测序结果。
subset <- subsetContig(combined, name = "sample", variables = c("PX", "PY"))
可视化
cloneCall +“gene”-使用由TCR/Ig +“nt”组成的基因-使用CDR3区域的核苷酸序列+“aa”-使用由CDR3区域的氨基酸序列+“gene+nt”-使用由TCR/Ig组成的基因+ CDR3区域的核苷酸序列。
关于克隆型的正确定义。
需要注意的是，克隆型基本上是利用两个位点的基因组合或nt/aa CDR3序列来命名的。在scRepertoire实现中，clonotype调用没有在CDR3序列中包含小的变化。因此，gene方法将是最敏感的，而使用nt或aa则是适度敏感的，并且对克隆型最具特异性的则是gene+nt。此外，克隆型calling试图合并两个位点，即TCRA和TCRB链，以及如果单个细胞条码有多个序列被识别(即一个细胞中表达2个TCRA链)。使用10x方法有一个条形码子集，只返回一个免疫受体链，未返回链被分配一个NA值。
研究克隆类型的第一个函数是quantContig()，它返回唯一克隆类型的总数或相对数量。scale + TRUE -独特克隆型的相对百分比由克隆型库的总大小scale + FALSE -报告独特克隆型的总数量。
library(ggsci)
quantContig(combined, cloneCall="gene+nt", scale = T) + scale_fill_npg()
在每个分析函数中，都可以导出用于创建可视化的数据框。要获得导出的值，使用exportTable == t。它将把数据返回到标记为functionName_output的全局环境中
quantContig_output <- quantContig(combined, cloneCall="gene+nt", scale = T, exportTable = T)
quantContig_output
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
我们也可以通过丰度来检查克隆型的相对分布。在这里，abundanceContig=()将生成一个线图，其中包含由示例或运行中的实例数量决定的总clonotypes数量。与上面一样，我们还可以使用函数中的group变量根据contig对象中的向量对其进行分组。
abundanceContig(combined, cloneCall = "gene", scale = F) + scale_fill_lancet()
abundanceContig(combined, cloneCall = "gene", exportTable = T)
# A tibble: 7,260 x 3
CTgene Abundance values
<chr> <int> <chr>
1 NA;TRAV1-1.TRAJ34.TRAC_TRBV28.TRBJ1-2.None.TRBC1 1 PY_P
2 NA;TRAV1-2.TRAJ20.TRAC_TRBV2.TRBJ2-7.TRBD2.TRBC2 1 PY_P
3 NA;TRAV1-2.TRAJ32.TRAC_TRBV2.TRBJ2-7.None.TRBC2 1 PY_P
4 NA;TRAV1-2.TRAJ37.TRAC_TRBV11-3.TRBJ2-1.TRBD2.TRBC2 1 PY_P
5 NA;TRAV1-2.TRAJ7.TRAC_TRBV20-1.TRBJ1-5.None.TRBC1 1 PY_P
6 NA;TRAV12-1.TRAJ30.TRAC_TRBV5-4.TRBJ2-2.None.TRBC2 1 PY_P
7 NA;TRAV12-1.TRAJ41.TRAC_TRBV27.TRBJ2-7.TRBD1.TRBC2 1 PY_P
8 NA;TRAV12-1.TRAJ49.TRAC_TRBV12-4.TRBJ2-5.TRBD1.TRBC2 1 PY_P
9 NA;TRAV12-1.TRAJ52.TRAC_TRBV6-5.TRBJ2-1.TRBD2.TRBC2 1 PY_P
10 NA;TRAV12-2.TRAJ13.TRAC_TRBV10-3.TRBJ2-7.TRBD1.TRBC2 1 PY_P
# ... with 7,250 more rows
我们可以通过调用lengtheContig()函数来查看CDR3序列的长度分布。重要的是，与其他基本可视化不同，克隆符号只能是“nt”或“aa”。由于调用clonotypes的方法如上所述，长度应该显示多模态曲线，这是对未返回的链序列和单个barcode中的多个链使用NA的结果。
lengthContig(combined, cloneCall="aa", chains = "combined") + scale_fill_npg()
lengthContig(combined, cloneCall="nt", chains = "single") + scale_fill_npg()
我们还可以使用compareClonotypes()函数来查看样本之间的克隆类型和动态变化。可以使用samples +根据列表元素的名称隔离特定的示例 图形+“alluvial”-下图图像+“area”-按各自克隆型面积绘制的图形 数+最上面的克隆型数来作图，这将根据个体样本的频率来计算。这也可以留空。clonotypes +可用于隔离特定的clonotype序列，确保调用匹配您想要可视化的序列。
compareClonotypes(combined, numbers = 10, samples = c("PX_P", "PX_T"),
cloneCall="aa", graph = "alluvial") + scale_fill_npg()
基本分析可视化的最后一点是TCR的vgenes的用法，使用vizVgenes()。T-GD: TCR1 == TCRG, TCR2 == TCRD
分面 ：x和方面。y +组合()对象中的任何变量
填充+颜色的bar通过变量
vizVgenes(combined, TCR="TCR1", facet.x = "sample", facet.y = "ID") + scale_fill_npg()
未完待续......