scRepertoire||单细胞免疫组库分析：R语言应用
前情回顾
10× Genomics单细胞免疫组库VDJ分析必知必会
免疫组库数据分析||immunarch教程：快速开始
免疫组库数据分析||immunarch教程：克隆型分析
免疫组库数据分析||immunarch教程：探索性数据分析
免疫组库数据分析||immunarch教程：载入10X数据
免疫组库数据分析||immunarch教程：GeneUsage分析
免疫组库数据分析||immunarch教程：Diversity 分析
免疫组库数据分析||immunarch教程：Clonotype tracking
免疫组库数据分析||immunarch教程：Clonotypes annotation
免疫组库数据分析||immunarch教程：Kmer 与 Motif 分析
在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 and preprint 2）。这些样本由成对的外周血和肿瘤浸润组成，有效地为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()
克隆空间内稳态
通过检查克隆空间，我们可以有效地观察特定比例下克隆所占的相对空间。另一种思考方法是把整个免疫受体测序看作一个量杯(总量)。在这个杯子里，我们将填充不同粘度的液体-或不同数量的克隆比例。克隆空间内稳态是问以不同比例(或者是不同粘度的液体，类推)填充的克隆杯子的百分比是多少。比例切点在函数中cloneType变量下设置，可以进行调整，以下为基准:
cloneTypes + Rare = .0001 + Small = .001 + Medium = .01 + Large = .1 + Hyperexpanded = 1
?clonalHomeostasi
clonalHomeostasis(combined, cloneCall = "gene")+ scale_fill_npg()
克隆的比例
与上述克隆空间稳态一样，克隆比例作用是将克隆放入单独的箱（bin）中。关键的区别在于，clonalProportion()函数不会查看克隆与总数的相对比例，而是根据总数对克隆进行排序，并将它们放入箱子(bin)中。
split表示克隆型按拷贝数或出现频率的排序，即1:10为每个样本中前10个克隆型。默认的bin位于函数中的split变量之下，可以进行调整，但在基线上它们如下所示。
split + 10 + 100 + 1000 + 10000 + 30000 + 100000
clonalProportion(combined, cloneCall = "gene") + scale_fill_npg()
如果您对加载到scRepertoire中的样本之间的相似性度量感兴趣，那么使用clonalOverlap()可以帮助可视化。目前有两种方法可用于clonalOverlap()重叠系数和莫里西塔指数。前者是在较小的样本中以独特的克隆型的长度来衡量克隆型的重叠。森里西塔指数更为复杂，它是对一个种群中个体分散程度的一种生态测量，包含了种群的规模。
clonalOverlap(combined, cloneCall = "gene+nt", method = "morisita")
scRepertoire的另一个新特性是使用powerTCR R包改编的clonesizeDistribution()通过克隆大小分布对示例进行聚类。如果用这个函数来分析样本克隆大小分布的相似性，请阅读并引用相应的引文。在这个函数中，method是指层次聚类的方法。
clonesizeDistribution(combined, cloneCall = "gene+nt", method="ward.D2")
多样性分析
多样性也可以通过样本或其他变量来衡量。多样性是通过四个度量来计算的:1)Shannon, 2) inverse Simpson, 3) Chao1和4)based Coverage Estimator (ACE)。前两者一般用于估计基线多样性，Chao/ACE指数用于估计样本的丰富度。
?clonalDiversity
clonalDiversity(combined, cloneCall = "gene", group = "samples")+ scale_fill_npg()
与表达谱数据整合
这里我就换一套数据来展示了，因为原教程的数据没有下载到，所以这里有用不知道从哪来的一套VDJ和RNA的数据来展示。需要说明的是，这里不是多个样本，而是一个样本了。也就是常说的，做了5‘转录组和VDJ 的项目。
csv1 <- read.csv("F:\\VDJ\\filtered_contig_annotations.csv", stringsAsFactors = F)
mycsv1 = list(ex=csv1)
?combineBCR
combined <- combineBCR(mycsv1, samples = c("PY"), ID = c("P"))
head(combined$PY_P)
barcode sample ID IGH cdr3_aa1
1 PY_P_ACGCAGCCATTCCTCG-1 PY P IGHV1-18.IGHJ4.None.IGHM CARDRKDVLDHW
2 PY_P_AGTAGTCGTAATAGCA-1 PY P IGHV4-61.IGHJ2.None.IGHM CARDKSCINGVCHRSWYFDLW
3 PY_P_TACGGATGTCAGCTAT-1 PY P IGHV5-51.IGHJ4.None.IGHM CARRGNGEMATYW
4 PY_P_AACTCAGTCAAACCAC-1 PY P IGHV2-5.IGHJ4.IGHD2-21.IGHM CAHRLNTYCRGDCYASFDYW
5 PY_P_GTCTCGTCAAAGAATC-1 PY P IGHV3-23.IGHJ3.None.IGHM CAKDLRGTRVTSYDAFDIW
6 PY_P_CGGAGTCAGCACCGCT-1 PY P IGHV4-39.IGHJ4.None.IGHM CATEYSSSPALW
cdr3_nt1 IGLC
1 TGTGCGAGAGATAGGAAAGATGTCCTGGACCACTGG IGLV2-8.IGLJ2.IGLC2
2 TGTGCGAGAGACAAATCTTGTATTAATGGTGTGTGCCACAGGTCTTGGTACTTCGATCTCTGG IGLV2-8.IGLJ1.IGLC1
3 TGTGCGCGACGTGGAAATGGAGAGATGGCTACTTATTGG IGLV2-14.IGLJ1.IGLC1
4 TGTGCACACCGCCTCAACACATATTGTCGTGGTGACTGCTACGCCTCCTTTGACTACTGG IGLV2-14.IGLJ2.IGLC2
5 TGTGCGAAAGATCTTCGGGGGACTAGGGTGACTTCTTATGATGCTTTTGATATCTGG IGLV2-14.IGLJ2.IGLC2
6 TGTGCGACAGAATATAGTAGCTCGCCGGCCCTCTGG IGLV2-14.IGLJ2.IGLC2
cdr3_aa2 cdr3_nt2
1 CNSYAGGNKVF TGCAACTCATATGCAGGCGGCAACAAGGTATTC
2 CNSYVDNNSYVF TGCAACTCCTATGTAGACAACAACAGTTATGTCTTC
3 CNSYTRSTTLYVF TGCAATTCATATACACGCAGCACCACTCTTTATGTCTTC
4 CTSYTSSSTLLF TGCACCTCATATACAAGCAGCAGCACTCTGCTGTTC
5 CTSYSSTYTLVLI TGCACCTCATATTCAAGCACCTACACTCTCGTGCTGATC
6 CSSHATSSTAYVIF TGCAGCTCACACGCAACCAGCAGCACTGCCTACGTGATATTC
CTgene
1 IGHV1-18.IGHJ4.None.IGHM_IGLV2-8.IGLJ2.IGLC2
2 IGHV4-61.IGHJ2.None.IGHM_IGLV2-8.IGLJ1.IGLC1
3 IGHV5-51.IGHJ4.None.IGHM_IGLV2-14.IGLJ1.IGLC1
4 IGHV2-5.IGHJ4.IGHD2-21.IGHM_IGLV2-14.IGLJ2.IGLC2
5 IGHV3-23.IGHJ3.None.IGHM_IGLV2-14.IGLJ2.IGLC2
6 IGHV4-39.IGHJ4.None.IGHM_IGLV2-14.IGLJ2.IGLC2
CTnt
1 TGTGCGAGAGATAGGAAAGATGTCCTGGACCACTGG_TGCAACTCATATGCAGGCGGCAACAAGGTATTC
2 TGTGCGAGAGACAAATCTTGTATTAATGGTGTGTGCCACAGGTCTTGGTACTTCGATCTCTGG_TGCAACTCCTATGTAGACAACAACAGTTATGTCTTC
3 TGTGCGCGACGTGGAAATGGAGAGATGGCTACTTATTGG_TGCAATTCATATACACGCAGCACCACTCTTTATGTCTTC
4 TGTGCACACCGCCTCAACACATATTGTCGTGGTGACTGCTACGCCTCCTTTGACTACTGG_TGCACCTCATATACAAGCAGCAGCACTCTGCTGTTC
5 TGTGCGAAAGATCTTCGGGGGACTAGGGTGACTTCTTATGATGCTTTTGATATCTGG_TGCACCTCATATTCAAGCACCTACACTCTCGTGCTGATC
6 TGTGCGACAGAATATAGTAGCTCGCCGGCCCTCTGG_TGCAGCTCACACGCAACCAGCAGCACTGCCTACGTGATATTC
CTaa CTstrict cellType
1 CARDRKDVLDHW_CNSYAGGNKVF IGH.53_IGHV1-18_IGLC51_IGLV2-8 B
2 CARDKSCINGVCHRSWYFDLW_CNSYVDNNSYVF IGH.98_IGHV4-61_IGLC99_IGLV2-8 B
3 CARRGNGEMATYW_CNSYTRSTTLYVF IGH.534_IGHV5-51_IGLC489_IGLV2-14 B
4 CAHRLNTYCRGDCYASFDYW_CTSYTSSSTLLF IGH.12_IGHV2-5_IGLC13_IGLV2-14 B
5 CAKDLRGTRVTSYDAFDIW_CTSYSSTYTLVLI IGH.504_IGHV3-23_IGLC465_IGLV2-14 B
6 CATEYSSSPALW_CSSHATSSTAYVIF IGH.287_IGHV4-39_IGLC268_IGLV2-14 B
我们观察到Barcode被改变名称了，所以我们：
combined$PY_P <- stripBarcode(combined$PY_P , column = 1, connector = "_", num_connects = 3)
读入RNA数据并执行Seurat标准流程：
?Read10X
seo<- Read10X(data.dir = 'F:\\VDJ\\filtered_feature_bc_matrix')
colnames(seo)
srto <- CreateSeuratObject(counts = seo$`Gene Expression`)
srto ->pbmc
pbmc <- NormalizeData(object = pbmc)
pbmc <- FindVariableFeatures(object = pbmc)
pbmc <- ScaleData(object = pbmc)
pbmc <- RunPCA(object = pbmc)
pbmc <- FindNeighbors(object = pbmc)
pbmc <- FindClusters(object = pbmc)
pbmc <- RunTSNE(object = pbmc)
pbmc <- RunUMAP(object = pbmc,dims = 1:15)
pbmc -> seurat
seurat
An object of class Seurat
33538 features across 7231 samples within 1 assay
Active assay: RNA (33538 features, 2000 variable features)
3 dimensional reductions calculated: pca, tsne, umap
DimPlot(seurat, label = T) + NoLegend() + scale_fill_npg()
接下来，我们可以获取clonotypic信息，并使用combineExpression()函数将其附加到Seurat对象中。重要的是，附件的主要需求是匹配contig细胞条形码和seurat或SCE对象的元数据行名中的条形码。如果这些不匹配，将失败。方便起见，我们建议您对Seurat对象行名称进行更改。
?combineExpression
seurat1 <- combineExpression(combined, seurat, cloneCall="gene")
作为核心函数，这个源代码我们是要背诵的了
combineExpression
function (df, sc, cloneCall = "gene+nt", groupBy = "none",
cloneTypes = c(None = 0, Single = 1, Small = 5, Medium = 20,
Large = 100, Hyperexpanded = 500), filterNA = FALSE)
{
df <- checkList(df)
cloneCall <- theCall(cloneCall)
Con.df <- NULL
meta <- grabMeta(sc)
cell.names <- rownames(meta)
if (groupBy == "none") {
for (i in seq_along(df)) {
data <- data.frame(df[[i]], stringsAsFactors = FALSE)
data2 <- unique(data[, c("barcode", cloneCall)])
data2 <- na.omit(data2[data2[, "barcode"] %in%
cell.names, ])
data2 <- data2 %>% group_by(data2[, cloneCall]) %>%
summarise(Frequency = n())
colnames(data2)[1] <- cloneCall
data <- merge(data, data2, by = cloneCall, all = TRUE)
Con.df <- rbind.data.frame(Con.df, data)
}
}
else if (groupBy != "none") {
data <- data.frame(bind_rows(df), stringsAsFactors = FALSE)
data2 <- na.omit(unique(data[, c("barcode", cloneCall,
groupBy)]))
data2 <- data2[data2[, "barcode"] %in% cell.names,
]
data2 <- as.data.frame(data2 %>% group_by(data2[, cloneCall],
data2[, groupBy]) %>% summarise(Frequency = n()))
colnames(data2)[c(1, 2)] <- c(cloneCall, groupBy)
x <- unique(data[, groupBy])
for (i in seq_along(x)) {
sub1 <- subset(data, data[, groupBy] == x[i])
sub2 <- subset(data2, data2[, groupBy] == x[i])
merge <- merge(sub1, sub2, by = cloneCall)
Con.df <- rbind.data.frame(Con.df, merge)
}
}
Con.df$cloneType <- NA
for (x in seq_along(cloneTypes)) {
names(cloneTypes)[x] <- paste0(names(cloneTypes[x]),
" (", cloneTypes[x - 1], " < X <= ",
cloneTypes[x], ")")
}
for (i in 2:length(cloneTypes)) {
Con.df$cloneType <- ifelse(Con.df$Frequency > cloneTypes[i -
1] & Con.df$Frequency <= cloneTypes[i], names(cloneTypes[i]),
Con.df$cloneType)
}
PreMeta <- unique(Con.df[, c("barcode", "CTgene",
"CTnt", "CTaa", "CTstrict", "Frequency",
"cloneType")])
rownames(PreMeta) <- PreMeta$barcode
if (inherits(x = sc, what = "Seurat")) {
sc <- AddMetaData(sc, PreMeta)
}
else {
rownames <- rownames(colData(sc))
colData(sc) <- cbind(colData(sc), PreMeta[rownames, ])[,
union(colnames(colData(sc)), colnames(PreMeta))]
rownames(colData(sc)) <- rownames
}
if (filterNA == TRUE) {
sc <- filteringNA(sc)
}
return(sc)
}
<bytecode: 0x00000224fac95a18>
<environment: namespace:scRepertoire>
head(seurat1@meta.data)
orig.ident nCount_RNA nFeature_RNA RNA_snn_res.0.8 seurat_clusters barcode
AAACCTGAGATCTGAA-1 SeuratProject 4386 1206 3 3 <NA>
AAACCTGAGGAACTGC-1 SeuratProject 6258 1723 4 4 <NA>
AAACCTGAGGAGTCTG-1 SeuratProject 4243 1404 1 1 <NA>
AAACCTGAGGCTCTTA-1 SeuratProject 5205 1622 4 4 <NA>
AAACCTGAGTACGTTC-1 SeuratProject 7098 2123 11 11 <NA>
AAACCTGGTATCACCA-1 SeuratProject 1095 604 10 10 <NA>
CTgene CTnt CTaa CTstrict Frequency cloneType highlight
AAACCTGAGATCTGAA-1 <NA> <NA> <NA> <NA> NA <NA> <NA>
AAACCTGAGGAACTGC-1 <NA> <NA> <NA> <NA> NA <NA> <NA>
AAACCTGAGGAGTCTG-1 <NA> <NA> <NA> <NA> NA <NA> <NA>
AAACCTGAGGCTCTTA-1 <NA> <NA> <NA> <NA> NA <NA> <NA>
AAACCTGAGTACGTTC-1 <NA> <NA> <NA> <NA> NA <NA> <NA>
AAACCTGGTATCACCA-1 <NA> <NA> <NA> <NA> NA <NA> <NA>
为什么那么多NA？ 我很是困惑，不会做错了吧？
table(rownames(seurat1@meta.data) %in% combined$PY_P$barcode)
FALSE TRUE
6469 762
head(seurat1@meta.data %>% filter(CTgene != "<NA>"))
orig.ident nCount_RNA nFeature_RNA RNA_snn_res.0.8
AAACCTGTCACTGGGC-1 SeuratProject 7539 2215 7
AAACCTGTCAGGTAAA-1 SeuratProject 1017 532 12
AAACCTGTCTCTTATG-1 SeuratProject 5974 1838 7
AAAGATGAGCGATTCT-1 SeuratProject 5352 1434 7
AAAGCAAAGAGGTTAT-1 SeuratProject 5878 1899 8
AAAGTAGCAATCCAAC-1 SeuratProject 5083 1463 7
seurat_clusters barcode
AAACCTGTCACTGGGC-1 7 AAACCTGTCACTGGGC-1
AAACCTGTCAGGTAAA-1 12 AAACCTGTCAGGTAAA-1
AAACCTGTCTCTTATG-1 7 AAACCTGTCTCTTATG-1
AAAGATGAGCGATTCT-1 7 AAAGATGAGCGATTCT-1
AAAGCAAAGAGGTTAT-1 8 AAAGCAAAGAGGTTAT-1
AAAGTAGCAATCCAAC-1 7 AAAGTAGCAATCCAAC-1
CTgene
AAACCTGTCACTGGGC-1 IGHV4-59.IGHJ5.None.IGHM_IGKV1D-33.IGKJ4.IGKC
AAACCTGTCAGGTAAA-1 NA_IGKV3-15.IGKJ2.IGKC
AAACCTGTCTCTTATG-1 IGHV7-4-1.IGHJ4.None.IGHM_IGLV3-10.IGLJ3.IGLC3
AAAGATGAGCGATTCT-1 IGHV4-59.IGHJ4.None.IGHA1_NA
AAAGCAAAGAGGTTAT-1 IGHV3-74.IGHJ3.None.IGHM_IGLV2-23.IGLJ3.IGLC3
AAAGTAGCAATCCAAC-1 IGHV3-23.IGHJ4.None.IGHG2_IGLV1-44.IGLJ3.IGLC3
CTnt
AAACCTGTCACTGGGC-1 TGTGCGAGAGGCGGGAACAGTGGCTTAGACCCCTGG_TGTCAACAGTATGATAATCTCCCGCTCACTTTC
AAACCTGTCAGGTAAA-1 NA_TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT
AAACCTGTCTCTTATG-1 TGTGCGAGAGGCGAACTAGAGAGGACGGTGACTACGGACTACTGG_TGTTACTCAACAGACAGCAGTGATAATCTGGGGGTGTTC
AAAGATGAGCGATTCT-1 TGTGCGAGGGAAGGCGCGCAGTACCTGGTACTTGAGTACTGG_TGTCAGCAGTATAATAACTGGCCTTTCACCTTC
AAAGCAAAGAGGTTAT-1 TGTGCAAGAGATCTTCCCGGTGCTTTTGATATCTGG_TGCTGCTCATATGCAGGTAGTAGCACCCGGGTGTTC
AAAGTAGCAATCCAAC-1 TGTGCGAAAGATTGGGTGGACAACAGTGAGTGGTCCTCATTCGCCGTTTTTGACTACTGG_TGTGCAACATGGGATGACAGCCTGAAAGGTTGGGTGTTT
CTaa
AAACCTGTCACTGGGC-1 CARGGNSGLDPW_CQQYDNLPLTF
AAACCTGTCAGGTAAA-1 NA_CQQYDNWPPYTF
AAACCTGTCTCTTATG-1 CARGELERTVTTDYW_CYSTDSSDNLGVF
AAAGATGAGCGATTCT-1 CAREGAQYLVLEYW_CQQYNNWPFTF
AAAGCAAAGAGGTTAT-1 CARDLPGAFDIW_CCSYAGSSTRVF
AAAGTAGCAATCCAAC-1 CAKDWVDNSEWSSFAVFDYW_CATWDDSLKGWVF
CTstrict Frequency
AAACCTGTCACTGGGC-1 IGH.1_IGHV4-59_IGLC1_IGKV1D-33 1
AAACCTGTCAGGTAAA-1 NA_NA_IGLC2_IGKV3-15 12
AAACCTGTCTCTTATG-1 IGH.2_IGHV7-4-1_IGLC3_IGLV3-10 1
AAAGATGAGCGATTCT-1 IGH.3_IGHV4-59_IGLC4_IGKV3-15 1
AAAGCAAAGAGGTTAT-1 IGH.4_IGHV3-74_IGLC5_IGLV2-23 1
AAAGTAGCAATCCAAC-1 IGH.5_IGHV3-23_IGLC6_IGLV1-44 1
cloneType highlight
AAACCTGTCACTGGGC-1 Single (0 < X <= 1) <NA>
AAACCTGTCAGGTAAA-1 Medium (5 < X <= 20) Clonotype2
AAACCTGTCTCTTATG-1 Single (0 < X <= 1) <NA>
AAAGATGAGCGATTCT-1 Single (0 < X <= 1) <NA>
AAAGCAAAGAGGTTAT-1 Single (0 < X <= 1) <NA>
AAAGTAGCAATCCAAC-1 Single (0 < X <= 1) <NA>
所以很可能是我们的数据中没有那么多B细胞，很多匹配不上啊。
slot(seurat1, "meta.data")$cloneType <- factor(slot(seurat1, "meta.data")$cloneType,
levels = c("Hyperexpanded (100 < X <= 500)", "Large (20 < X <= 100)",
"Medium (5 < X <= 20)", "Small (1 < X <= 5)",
"Single (0 < X <= 1)", NA))
DimPlot(seurat1, group.by = "cloneType") +
scale_color_manual(values = colorblind_vector(5), na.value="grey")
为了验证我们的猜想，我们做一个差异分析看看特征基因不就知道了，主要是7,9 群，如果他们果然是B细胞那就很好了。或者做一个B细胞的marker 点图。
library(Seurat)
library(SeuratData)
library(ggplot2)
data("pbmc3k.final")
modify_vlnplot<- function(obj,
feature,
pt.size = 0,
plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm"),
...) {
p<- VlnPlot(obj, features = feature, pt.size = pt.size, ... ) +
xlab("") + ylab(feature) + ggtitle("") +
theme(legend.position = "none",
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.y = element_line(),
axis.title.y = element_text(size = rel(1), angle = 0, vjust = 0.5),
plot.margin = plot.margin )
return(p)
}
## main function
StackedVlnPlot<- function(obj, features,
pt.size = 0,
plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm"),
...) {
plot_list<- purrr::map(features, function(x) modify_vlnplot(obj = obj,feature = x, ...))
plot_list[[length(plot_list)]]<- plot_list[[length(plot_list)]] +
theme(axis.text.x=element_text(), axis.ticks.x = element_line())
p<- patchwork::wrap_plots(plotlist = plot_list, ncol = 1)
return(p)
}
my36colors <-c('#E5D2DD', '#53A85F', '#F1BB72', '#F3B1A0', '#D6E7A3', '#57C3F3', '#476D87',
'#E95C59', '#E59CC4', '#AB3282', '#23452F', '#BD956A', '#8C549C', '#585658',
'#9FA3A8', '#E0D4CA', '#5F3D69', '#C5DEBA', '#58A4C3', '#E4C755', '#F7F398',
'#AA9A59', '#E63863', '#E39A35', '#C1E6F3', '#6778AE', '#91D0BE', '#B53E2B',
'#712820', '#DCC1DD', '#CCE0F5', '#CCC9E6', '#625D9E', '#68A180', '#3A6963',
'#968175'
)
classmk <-c("IL7R", "CCR7", "IL7R", "S100A4","CD14", "LYZ","MS4A1","CD8A","FCGR3A", "MS4A7", "GNLY",
"NKG7","FCER1A", "CST3","PPBP")
StackedVlnPlot(seurat1, c(classmk), pt.size=0, cols=my36colors)
我们知道，MS4A1是B细胞的一个明显的marker，果不其然，这俩群是B细胞没错了。然后我们再看一下差异分析的结果：
seurat1@misc$mk <- FindAllMarkers(seurat1,only.pos = T)
seurat1@misc$mk %>% filter(cluster == 9) %>% top_n(20,avg_logFC ) -> topgene
topgene
p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
TCL1A 0.000000e+00 2.670484 0.947 0.014 0.000000e+00 9 TCL1A
CD79A1 0.000000e+00 2.245267 0.986 0.093 0.000000e+00 9 CD79A
IGHM1 0.000000e+00 2.178957 0.961 0.069 0.000000e+00 9 IGHM
FCER21 0.000000e+00 2.031561 0.870 0.033 0.000000e+00 9 FCER2
MS4A11 0.000000e+00 1.987625 0.975 0.084 0.000000e+00 9 MS4A1
LINC009261 0.000000e+00 1.871576 0.863 0.060 0.000000e+00 9 LINC00926
PLPP51 0.000000e+00 1.708546 0.803 0.093 0.000000e+00 9 PLPP5
IGHD1 0.000000e+00 1.635674 0.722 0.027 0.000000e+00 9 IGHD
CD79B1 1.925974e-299 1.784911 0.937 0.151 6.459332e-295 9 CD79B
HLA-DQA11 5.403560e-293 1.778771 0.901 0.139 1.812246e-288 9 HLA-DQA1
HLA-DQB11 8.293830e-278 1.872821 0.975 0.189 2.781585e-273 9 HLA-DQB1
HLA-DRA3 3.962852e-230 1.945424 0.996 0.279 1.329061e-225 9 HLA-DRA
HLA-DQA23 2.739902e-219 1.611850 0.954 0.231 9.189083e-215 9 HLA-DQA2
HLA-DRB53 8.598615e-212 1.690653 0.986 0.272 2.883803e-207 9 HLA-DRB5
HLA-DRB13 4.137877e-208 1.765784 1.000 0.305 1.387761e-203 9 HLA-DRB1
HLA-DPB12 2.529860e-181 1.673158 0.993 0.369 8.484643e-177 9 HLA-DPB1
CD742 1.774900e-165 2.103404 1.000 0.847 5.952659e-161 9 CD74
IGHV3-21 1.128949e-111 1.904073 0.144 0.005 3.786270e-107 9 IGHV3-21
IGHV3-30 3.301247e-62 1.627276 0.113 0.007 1.107172e-57 9 IGHV3-30
IGKV3-201 1.622764e-42 1.774558 0.120 0.013 5.442425e-38 9 IGKV3-20
seurat1@misc$mk %>% filter(cluster == 7) %>% top_n(20,avg_logFC ) -> topgene
topgene
clu2featur(seurat1,idents = 7 ,features = topgene$gene[1:8])
差异分析的结果，也证实了我们的猜想：主要有cluster 7，9 是B细胞
现在我们可以通过首先将克隆型作为一个因素排序来查看克隆型容器的分布，这可以防止颜色按字母顺序排列。接下来，我们使用Seurat中的DimPlot()函数调用和scale_color_manual附加图层。
seurat1 <- highlightClonotypes(seurat1, cloneCall= "aa", sequence = c("CAKDRAEWLPQYYYGMDVW_CQQYDNLPPLFTF", "NA_CQQYDNWPPYTF"))
DimPlot(seurat1, group.by = "highlight")
由于我们的克隆型实在太少，所以看的不是很清楚：
table(table(seurat1$CTgene))
1 2 3 4 5 12 14
646 30 4 2 2 1 1
我们还可以使用occupiedscRepertoire()函数并选择x.axis来显示细胞对象的元数据中的集群或其他变量，从而查看分配到特定频率范围的集群的单元数。
occupiedscRepertoire(seurat1, x.axis = "seurat_clusters")
最后，在修改了所有元数据之后，我们可以使用alluvialClonotypes()函数查看跨多个类别的clonotypes。要理解这种绘图方法的基本概念，我强烈推荐阅读这篇文章this post本质上我们能够使用绘图来检查分类变量的流向。因为这个函数将生成一个图形，其中每个克隆型按所谓的分层进行排列，这将花费一些时间，具体时间取决于细胞总数的大小。为了加快速度，我们实际上将在使用alluvialClonotypes()之前对seurat对象进行子集化。
alluvialClonotypes(seurat1, cloneCall = "gene",
y.axes = c("Frequency", "cloneType"),
color = "seurat_clusters")
对于希望在Seurat对象中使用元数据来执行scRepertoire提供的分析的用户来说，还有一个选项是使用expression2List()函数，该函数将获取元数据并按亚群将数据输出为列表。也就是人们常说的，细胞类型特异的克隆型分析。
combined2 <- expression2List(seurat1, group = "seurat_clusters")
length(combined2) #now listed by cluster
clonalDiversity(combined2, cloneCall = "nt")
clonalHomeostasis(combined2, cloneCall = "nt")
clonalOverlap(combined2, cloneCall="aa", method="overlap")
这是对scRepertoire从最初处理和可视化到附加到Seurat对象中的mRNA表达值的能力的一般概述。如果您有任何问题、意见或建议，请访问github：https://ncborcherding.github.io/vignettes/vignette.html
