immunarch — Fast and Seamless Exploration of Single-cell and Bulk T-cell/Antibody Immune Repertoires in R
10× Genomics单细胞免疫组库VDJ分析必知必会
免疫组库数据分析||immunarch教程：克隆型分析
免疫组库数据分析||immunarch教程：探索性数据分析
免疫组库数据分析||immunarch教程：载入10X数据
免疫组库数据分析||immunarch教程：快速开始
免疫组库数据分析||immunarch教程：GeneUsage分析
免疫组库数据分析||immunarch教程：Diversity 分析
今天，我们继续我们的免疫组库数据分析的Demos，这一次我们来谈谈Clonotype tracking分析。像我这样刚入门免疫组库的人首先会问什么是Clonotype tracking？
克隆型追踪(Clonotype tracking)是监测疫苗接种和癌症免疫感兴趣的克隆型频率变化的常用方法。例如，研究人员可以在疫苗接种前和接种后的不同时间点跟踪克隆类型，或分析肿瘤样本中恶性克隆型的生长。所谓的追踪用到的可视化方法类似我们提到过的桑基图，具体地：桑基图在单细胞数据探索中的应用。immunarch中集成了多种克隆型跟踪方法。目前有三种方法可供选择。trackClonotypes的输出可以立即用vis函数可视化。
追踪最丰富的克隆型
最简单的方法是从一个输入免疫序列中选择最丰富的克隆类型，并批量跟踪所有的免疫序列。参数.which和.col用于选择免疫序列、从中获取的克隆类型数量以及使用的列。从第一个库中选择5个最丰富的克隆型，并利用其CDR3核苷酸序列对其进行跟踪:
library(immunarch); data(immdata) # Load the package and the test dataset
?trackClonotypes
tc1 <- trackClonotypes(immdata$data, list(1, 5), .col = "nt")
tc1
CDR3.nt
1: TGCGCCAGCAGCCAAGAAGGGACAGGGTATTCCGGGGAGCTGTTTTTT
2: TGCGCCAGCAGCTACAGGGTTGGCACAGATACGCAGTATTTT
3: TGTGCCACCAGCACCAACAGGGGCGGAACCCCAGCAGATACGCAGTATTTT
4: TGTGCCACCAGCATCGGAGGCGGGAGCTACGAGCAGTACTTC
5: TGTGCCAGCAGTCCTTGGACAGGGAGTATGGCCCTCCACTTT
A2-i129 A2-i131 A2-i133 A2-i132 A4-i191 A4-i192 MS1
1: 0.020352941 0 0 0 0 0 0
2: 0.019176471 0 0 0 0 0 0
3: 0.007764706 0 0 0 0 0 0
4: 0.006352941 0 0 0 0 0 0
5: 0.005647059 0 0 0 0 0 0
MS2 MS3 MS4 MS5 MS6
1: 0 0 0 0 0
2: 0 0 0 0 0
3: 0 0 0 0 0
4: 0 0 0 0 0
5: 0 0 0 0 0
which参数(第二个参数)的值list(1,5)意味着从输入字符串immdata$data的第一个列表中选择5个clonotypes。col参数的值“nt”意味着函数应该只接受CDR3核苷酸序列。
从“MS1”库中选择10个最丰富的氨基酸克隆型序列及其V基因进行跟踪:
tc2 <- trackClonotypes(immdata$data, list("MS1", 10), .col = "aa+v")
tc2
CDR3.aa V.name A2-i129 A2-i131 A2-i133
1: CASSFEGAMDTQYF TRBV7-6 0 0 0
2: CASSLGDSTYEQYF TRBV5-6 0 0 0
3: CASSLGLREQGETQYF TRBV28 0 0 0
4: CASSLQAGGNTDTQYF TRBV7-2 0 0 0
5: CASSLYSNEQFF TRBV7-9 0 0 0
6: CASSVYSTISEQYF TRBV9 0 0 0
7: CSARDLANSYEQYF TRBV20-1 0 0 0
8: CSTEEDSYNEQFF TRBV20-1 0 0 0
9: CSVELRTESGYEQYF TRBV29-1 0 0 0
10: CSYRTGGPEQYF TRBV29-1 0 0 0
A2-i132 A4-i191 A4-i192 MS1 MS2
1: 0 0 0 0.008941176 0.0000000000
2: 0 0 0 0.011529412 0.0000000000
3: 0 0 0 0.007058824 0.0001176471
4: 0 0 0 0.063529412 0.0000000000
5: 0 0 0 0.004470588 0.0000000000
6: 0 0 0 0.037647059 0.0000000000
7: 0 0 0 0.009529412 0.0000000000
8: 0 0 0 0.024000000 0.0000000000
9: 0 0 0 0.017764706 0.0000000000
10: 0 0 0 0.007294118 0.0000000000
MS3 MS4 MS5 MS6
1: 0.0000000000 0.0000000000 0 0
2: 0.0000000000 0.0001176471 0 0
3: 0.0000000000 0.0000000000 0 0
4: 0.0000000000 0.0000000000 0 0
5: 0.0000000000 0.0000000000 0 0
6: 0.0001176471 0.0001176471 0 0
7: 0.0000000000 0.0000000000 0 0
8: 0.0001176471 0.0000000000 0 0
9: 0.0000000000 0.0000000000 0 0
10: 0.0000000000 0.0000000000 0 0
参数的值list("MS1"， "10")，它的意思是选择10个clonotypes从命名为"MS1"的指令集在指令集immdata$data的输入列表中。col的“aa+v”值意味着该功能应该同时取CDR3氨基酸序列和最丰富的克隆型的v基因片段。
同样经典而又傻瓜式地操作
p1 <- vis(tc1)
p2 <- vis(tc2)
p1 / p2
用特定的核苷酸或氨基酸序列追踪克隆型
为了跟踪特定的clonotype序列，可以提供核苷酸或氨基酸序列作为which参数，同时提供列.col，以指定在哪些列中搜索序列。例如，要跟踪下面指定的七个CDR3氨基酸序列，你需要执行以下代码:
target <- c("CASSLEETQYF", "CASSDSSGGANEQFF", "CASSDSSGSTDTQYF", "CASSLAGGYNEQFF", "CASSDSAGGTDTQYF", "CASSLDSYEQYF", "CASSSAGGYNEQFF")
tc <- trackClonotypes(immdata$data, target, .col = "aa")
vis(tc)
利用特定序列和基因片段追踪克隆型
与之前的方法相比，利用序列和基因片段的信息来追踪克隆型是可能的。你有一个特定的CDR3序列和基因片段的数据框。我们将通过从批次的第一个清单中选择10个最丰富的克隆类型来模拟这一点:
target <- immdata$data[[1]] %>%
select(CDR3.aa, V.name) %>%
head(10)
target
# A tibble: 10 x 2
CDR3.aa V.name
<chr> <chr>
1 CASSQEGTGYSGELFF TRBV4-1
2 CASSYRVGTDTQYF TRBV4-1
3 CATSTNRGGTPADTQYF TRBV15
4 CATSIGGGSYEQYF TRBV15
5 CASSPWTGSMALHF TRBV27
6 CASQGDSFNSPLHF TRBV4-1
7 CASSQDMGGRNTGELFF TRBV4-1
8 CASSEEPRLFGYTF TRBV2
9 CASSQPGQGGGDEQFF TRBV4-1
10 CASSWVARGPYEQYF TRBV6-6
tc <- trackClonotypes(immdata$data, target)
vis(tc)
请注意，您可以使用目标数据框中的任何列，例如CDR3核苷酸和氨基酸序列以及任何基因片段。所以如何确定哪一列呢？
可视化跟踪
根据你的研究和审美需求，有三种可视化克隆型追踪的方法。要选择图形的类型，需要提供“。- .plot = "smooth" -默认使用，使用光滑线条和堆叠条形图进行可视化;- .plot = "area" -使用丰度线下面的区域来可视化丰度;- .plot =“line”-只可视化线，连接同一克隆型在时间点之间的丰度水平。
target <- c("CASSLEETQYF", "CASSDSSGGANEQFF", "CASSDSSGSTDTQYF", "CASSLAGGYNEQFF", "CASSDSAGGTDTQYF", "CASSLDSYEQYF", "CASSSAGGYNEQFF")
tc <- trackClonotypes(immdata$data, target, .col = "aa")
p1 <- vis(tc, .plot = "smooth")
p2 <- vis(tc, .plot = "area")
p3<- vis(tc, .plot = "line")
library(patchwork)
p1+ p2 +p3
改变样本的顺序
vis函数的order参数控制可视化中样本的顺序。您可以通过您计划可视化的样本索引或样本名称。
# Passing indices
names(immdata$data)[c(1, 3, 5)] # check sample names
vis(tc, .order = c(1, 3, 5))
# You can change the order
vis(tc, .order = c(5, 1, 3))
当然，这样也是可以的，图我就不贴了啊。
# Passing sample names
vis(tc, .order = c("A2-i129", "A2-i133", "A4-i191"))
如果元数据（metadata）包含有关时间的信息，如接种疫苗或肿瘤样本的时间点，则可以使用它相应地对样本重新排序。在我们的示例中，immdata$meta不包含关于时间点的信息，因此我们将模拟这种情况。
immdata$meta$Timepoint <- sample(1:length(immdata$data))
immdata$meta
# A tibble: 12 x 7
Sample ID Sex Age Status Lane Timepoint
<chr> <chr> <chr> <dbl> <chr> <chr> <int>
1 A2-i129 C1 M 11 C A 2
2 A2-i131 C2 M 9 C A 11
3 A2-i133 C4 M 16 C A 7
4 A2-i132 C3 F 6 C A 3
5 A4-i191 C8 F 22 C B 1
6 A4-i192 C9 F 24 C B 4
7 MS1 MS1 M 12 MS C 10
8 MS2 MS2 M 30 MS C 6
9 MS3 MS3 M 8 MS C 12
10 MS4 MS4 F 14 MS C 8
11 MS5 MS5 F 15 MS C 5
12 MS6 MS6 F 15 MS C 9
sample_order <- order(immdata$meta$Timepoint)
immdata$meta$Timepoint[sample_order]
immdata$meta$Sample[sample_order]
vis(tc, .order = sample_order)
vis(tc, .order = order(immdata$meta$Timepoint))
改变调色板
在R控制台中运行?scale_fill_brewer，了解ColorBrewer及其配色方案的更多信息
p1<- vis(tc) + scale_fill_brewer(palette = "Spectral")
p2 <- vis(tc) + scale_fill_brewer(palette = "RdBu")
p1 + p2
我们为什么要研究Clonotype tracking啊？是为了看Clonotype 在不同样本中的分布情况进一步看Clonotype的可能的迁移和转化状态，所以关键的是样本分组的生物学意义啊。
