tcR包：T细胞受体和免疫球蛋白数据进行高级分析和可视化（二）

导语

GUIDE ╲

免疫球蛋白(IG)和T细胞受体(TR)在适应性免疫应答过程中起着关键的抗原识别作用。上一次我们介绍到tcR包：T细胞受体和免疫球蛋白数据进行高级分析和可视化（一）。今天小编继续为大家介绍分析T细胞受体库的R包：tcR包，可以对TR序列进行多样性评估、共享T细胞受体序列识别、基因usage统计计算等。

R包使用

三、基因usage

2. 基因usage比较（Gene usage comparing）

为了评估组库中V-和J基因usage情况，该程序包实现了两种分析方法的子程序：信息论方法和PCA（主成分分析）方法。

（1）香农熵（Shannon entropy）和J-S散度（Jensen-Shannon divergence）

使用entropy函数可以评估基因usage的多样性。用函数kl.div实现Kullback-Leibler非对称测度，函数js.div和js.div.seg实现用Jensen-Shannon对称测度来评估不同组库的基因usage的距离，函数js.div用于计算给定分布之间的JS差异，而js.div.seg用于计算两个克隆集或列表的基因分布之间的JS差异。vis.radarlike 用来可视化距离。

①计算分布的香农熵

entropy(0:100, .laplace = 1)
#用拉普拉斯修正方法将0:100分布转换
#.laplace是拉普拉斯修正，.laplace = 1意思是在转换前对每个值加“1”

②计算列表中每个数据框的 V区片段usage的香农熵

entropy.seg(twb, HUMAN_TRBV)

③计算两个数据框之间的V-usage的JS差异

js.div.seg(twb[1:2], HUMAN_TRBV, .verbose = F)
#.verbose是否输出数据处理进度条

④可视化

imm.js <- js.div.seg(twb, HUMAN_TRBV, .verbose = F)
#计算距离
vis.radarlike(imm.js, .ncol = 2)
#每个数据框与其他数据框的距离结果绘制距离图

（2）主成分分析Principal Component Analysis (PCA)

主成分分析(PCA)是一种将一组观测值转换为一组特殊值进行分析的统计过程。使用pca.segments函数在 V-usage 或 J-usage上的基因片段频率数据执行PCA，返回PCA对象或绘制结果。函数pca.segments.2D是在VJ-usage上执行PCA。vis.pca绘制PCA结果。

①计算V-segment usage的PCA并绘图

pca.segments(twb, .cast.freq.seg = T)  
#.cast.freq.seg = T，将codegeneUsage应用于提供的数据

②返回数值

pca.segments(twb, .cast.freq.seg = T,.do.plot = F)  
#.do.plot如果是T绘制一个图形，否则返回一个pca对象

③vis.pca绘制PCA结果

tmp = geneUsage(twb)
vis.pca(prcomp(t(tmp[,-1])))

四、组库的重叠分析（Repertoire overlap analysis）

tcR提供了许多基于clonotypes之间共享的cloneset来评估相似度的函数，处理数据框数据。

1. 重叠量化（Overlap quantification）

repOverlap 函数是计算cloneset重叠系数的所有函数的通用接口

（1）共享克隆型的数量

评估两个克隆集相似性最直接但也最有效的方法是计算共享克隆类型的数量。函数intersectClonesets (repOverlap(your_data, 'exact'))默认使用“CDR3.nucleotide”计算共享克隆型的数量，但是用户可以通过使用参数 .type 或 .col来更改目标列。在函数find.clonotypes中，用户可以选择将哪种方法应用于元素：元素的精确匹配（exact）、Hamming距离匹配或Levenshtein距离匹配。参数.norm是通过将共享克隆类型的数量除以克隆集大小的乘积来执行标准化(强烈建议不要这样做，否则结果将与克隆集的大小相关)。

①CDR3评估两个克隆集相似性

repOverlap(twb[1:2], 'exact', .seq ='nuc', .verbose = F)
#使用 'exact'匹配方法
#.seq是选择序列类型， "nuc"代表使用CDR3的核苷酸序列
#.seq="aa" 代表使用CDR3的氨基酸序列
#.verbose是否输出程序进程#比较twb的前两个数据框

②twb数据框两两评估相似性

repOverlap(twb, 'exact', 'nuc', .norm = T, .verbose = F)
#.norm是否标准化

③使用氨基酸序列+Vgene评估

repOverlap(twb, 'exact', .seq='aa', .vgene = T, .verbose = F)
#.seq="aa" 代表使用CDR3的氨基酸序列
#.vgene = T，使用V基因计算共享或相似克隆型

④可视化：绘制相似性值热图

vis.heatmap(repOverlap(twb, 'exact', .seq='aa', 
                       .vgene = T, .verbose = F),
            .title = 'twb - (ave)-intersection',
            .labs = '')
#.labs横纵坐标名称

⑤其他函数

函数intersectCount, intersectLogic和intersectIndices在选择匹配的列上相对来说更灵活。它们都具有参数.col，用于指定将在交集计算中使用的列的名称。函数intersectCount返回相似元素的数量；intersectIndices(x, y)返回两列矩阵，第一列表示给定x中一个元素的索引，第二列表示y中的与x中的相对元素相似的元素的索引；intersectLogic(x, y)返回length(x)或nrow(x)的逻辑向量，其中位置i为TRUE表示在y中找到了索引为{i}的元素。

例：

imm.1.2 <- intersectLogic(twb[[1]], twb[[2]],
           .col = c('CDR3.amino.acid.sequence', 'V.gene'))  
#得到共享元素的逻辑向量，其中元素是指CDR3核苷酸序列的元组和相应的V-segment
head(twb[[1]][imm.1.2, c('CDR3.amino.acid.sequence', 'V.gene')])
#获取twb[[1]]和twb[[2]]中同时存在的元素

（2）“Top cross”

在最丰富的克隆型中，共有克隆型的数量可能与那些具有较少计数的克隆型显著不同。top.cross函数可以发现这些，应用tcR::intersectClonesets，比如.n == seq(1000, 100000, 1000)代表，到前1000个克隆型，2000个，3000个等等直到前100000个克隆。

举例：卡定不同数的克隆型，比较共享克隆型的数量

twb.top <- top.cross(.data = twb, .n = seq(500, 10000, 500), 
                     .verbose = F, .norm = T)
top.cross.plot(twb.top)

（3）更复杂的相似性度量

Tversky 指数 (克隆集用repOverlap(your_data, 'tversky')；向量用tversky.index)是集合上的非对称相似性度量，用于比较变体和原型。如果使用默认参数，它类似于Dice’s系数。

重叠系数 (克隆集用repOverlap(your_data, 'overlap')；向量用 overlap.coef) 是度量两个集合之间重叠的相似性度量，定义为交集的大小除以两个集合大小中较小的那个。

Jaccard 指数(克隆集用repOverlap(your_data, 'jaccard')；向量用 jaccard.index)是用来比较样本集的相似性和多样性的统计量。

Morisita’s重叠指数(克隆集用repOverlap(your_data, 'morisita')；向量用morisitas.index)是对群体中个体离散度的统计度量，用于比较样本之间的重叠。这个公式是基于增加样本的大小将增加多样性的假设，因为它将包括不同的栖息地(例如不同的动物群体)。

例：对每一对repertoires应用Morisitas重叠指数，使用V gene计算。比如当且仅当它们的CDR3 aa序列相等且它们的V基因相等时，一个CDR3克隆型与另一个CDR3克隆型是相等的。

repOverlap(twb, .method='morisita', .seq='aa', 
           .quant='read.count',
           .vgene = T, .verbose = F)
#对每一对repertoires应用Morisitas重叠指数
#方法可用.method = c("exact", "hamm", "lev", "jaccard","morisita", "tversky", "overlap", "horn")
#.quant使用哪一列来表示克隆类型的数量，
#.quant可选 c("read.count", "umi.count", "read.prop", "umi.prop")

2. 重叠统计量和测试（Overlap statistics and tests）

计算给定的重叠矩阵中的值的OZ-scores(“重叠Z分数”)，即对于每个值，计算离矩阵平均值的标准偏差数。

mat <- repOverlap(twb)
ozScore(mat)

3. 共享库（Shared repertoire）

函数shared.repertoire能够计算共享克隆型组库。shared.representation能够计算每个集合的共享克隆类型数量，以确定共享的程度(比如发现有一定克隆数量的人群数量)。函数shared.summary 相当于repOverlap(, 'exact')，但适用于共享的数据框。cosine.sharing函数利用共享序列计数向量的余弦相似度衡量集合之间的距离。

①例：计算在两个或两个以上的人中发现的氨基酸CDR3序列和V基因的共享库，并从输入列表中的每个数据框中返回此类克隆型的Read.count列。

imm.shared <- shared.repertoire(.data = twb, 
                 .type = 'avrc',
                 .min.ppl = 2, .verbose = F)
#.type表示如何创建共享库，是一个四个字母的字符串
#'avrc'中的第一个字母a表示使用CDR3氨基酸序列，若换成n表示核苷酸序列
#'avrc'中的第二个字母v表示是否使用V.gene列，若换成0代表不使用
#'avrc'中的第三个字母r表示选择带有数字字符的列时使用UMIs还是reads
#'avrc'中的第四个字母c表示要选择的列的名称作为序列的数字特征。
# "c" 代表"Umi.count"，"p"代表"Umi.proportion"， "r"代表"Rank"列，"i"代表"Index"列。
#.min.ppl代表至少有多少人必须拥有一个序列才能将这个序列保留在共享库中
head(imm.shared)

②计算共享序列的数量

shared.representation(imm.shared)

五、 多样性评估

一些函数用于评估克隆型在给定集合中的分布情况。用函数diversity 和inverse.simpson评估多样性。用函数gini 和 gini.simpson评估克隆分布的偏态性。函数repDiversity是用户应该使用哪个函数来估计克隆集的多样性的接口。函数diversity (repDiversity(your_clonesets, "div"))计算生态多样性指数。函数inverse.simpson (repDiversity(your_clonesets, "inv.simp"))计算逆辛普森指数。Inverse Simpson Index（如两个相似克隆类型的逆概率）。函数gini (repDiversity(your_clonesets, "gini")) 计算克隆分布的经济基尼指数Gini index。函数gini.simpson (repDiversity(your_clonesets, "gini.simp"))计算 Gini-Simpson index。函数chao1 (repDiversity(your_clonesets, "chao1")) 计算Chao1 index。函数repDiversity接受单个clonesets以及clonesets列表。函数.quant指定要使用哪一列来计算分集。

举例: 用生态多样性指数评估克隆的多样性

repDiversity(twb, .method='div', .quant='read.count')
#'div'是生态多样性指数
#.method = c("chao1", "gini.simp", "inv.simp", "gini","div", "entropy")
#.quant使用哪一列来表示克隆类型的数量

六、可视化

1. CDR3长度和读数分布图

vis.count.len绘制CDR3核苷酸序列长度图，vis.number.count绘制counts直方图。输入数据数据框或数据列表。

例：绘制CDR3核苷酸序列长度图

vis.count.len(twb[[1]], .name = "twb[[1]] CDR3 lengths", 
              .col = "Read.count")
#.col数据框的列数

2. 顶部比例条形图

函数vis.top.proportions可以实现对最丰富的克隆型比例的可视化。

例：不同的top数的克隆所占比例

vis.top.proportions(twb, .head=c(10, 500, 3000, 10000),
                   .col = "Read.count")

3. 克隆空间稳态条形图

vis.clonal.space函数可以可视化每组克隆类型占用了多少空间，并按数据中的比例将其分成组。可以将clonal.space.homeostasis的输出作为输出。

twb.space <- clonal.space.homeostasis(twb)
vis.clonal.space(twb.space)

4. 热图

集合的配对距离或相似度可以表示为二元矩阵，其中每一行和每一列表示一个克隆集。vis.heatmap用来可视化。

twb.shared <- repOverlap(twb, "exact", .norm = F, .verbose = F)
#使用 'exact'匹配来评估两个克隆集相似性
vis.heatmap(twb.shared, .title = "Twins shared nuc clonotypes",
            .labs = c("Sample in x", "Sample in y"),
            .legend = "# clonotypes")
#.legend图例名称

5. 雷达类图

vis.radarlike可以用来可视化距离，命名为雷达类图（因为它不是精确的雷达图）。

twb.js <- js.div.seg(twb, HUMAN_TRBV, .verbose = F)
#计算Jensen-Shannon差异来评估不同repertoires的基因usage的距离
vis.radarlike(twb.js, .ncol = 2)

6. 基因usage直方图

(这部分之前已经介绍过了)

7. PCA可视化

(这部分之前已经介绍过了)

8. Logo-like图

km <- get.kmers(twb[[1]]$CDR3.amino.acid.sequence, 
                .head = 100, .k = 7, .verbose = F)
#从给定的字符向量或数据框中获取kmers的向量
#K是代表kmer的大小，kmers是指将序列分为k个碱基的字符串
d <- kmer.profile(km)
#返回给定字符向量或数据框具有相同长度序列的配置文件
vis.logo(d)

七、突变网络

突变网络(或突变图)是一个图，顶点代表核苷酸或框内氨基酸序列（框外氨基酸序列在创建突变网络的时候会被过滤掉），边代表用hamming距离连接(parameter .method =‘hamm’) 或edit距离 (parameter .method =‘lev’) ，它们之间的距离不超过mutation.network函数中.max指定的距离。

twb.shared <- shared.repertoire(twb, .head = 1000, .verbose = F)
#twb数据框的共享集合
G <- mutation.network(twb.shared)
G

小编总结

tcR包是一款非常先进功能非常强大的T细胞受体分析工具包，让我们来总结一下它都能实现哪些功能吧！

（1）可以直接接受多种分析工具（如 MiTCR、MiGEC、 VDJtools、ImmunoSEQ、IMSEQ 和MiXCR）的输出数据，作为输入进行直接分析

（2）数据操作(框内/框外序列子集设置，克隆型motif搜索)

（3）进行一些描述性统计(读数、克隆型数、基因片段usage)

（4）可统计共享克隆型(共享克隆型的数量，是否使用V基因并入计算；最丰富的克隆类型之间的连续的交集(“top-cross”)

（5）进行组库的比较(Jaccard指数、Morisita重叠指数、Horn’s指数、Tversky指数、重叠系数)

（6）V基因和J基因usage和分析(PCA, Shannon Entropy, Jensen-Shannon Divergence)。

（7）多样性评估(生态多样性指数、Gini指数、逆Simpson指数、稀疏度分析)

（8）生成人工组库(目前仅beta链)

（9）免疫扫描图谱分析

（10）分析结果进行可视化

（11）形成突变网络(网络图中顶点表示CDR3核苷酸/氨基酸序列，边代表序列之间的相似度，使用 low hamming或edit distance距离测量)

END