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

导语

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

背景介绍

免疫组库是指T细胞受体和B细胞受体（也称为免疫球蛋白）的总和，它们构成了机体的适应性免疫系统。这些高度多样化的抗原受体可以识别“异己”并产生免疫反应。

新一代测序(NGS)技术开辟了IG和TR序列研究的新时代，应运而生了一些方法：

标准的IMGT/HighV-QUEST网页分析（IMGT是一个用新一代测序和深度测序分析IG和TR的web）
基于原始IG / TR测序数据的处理：从reads中提取互补决定区（CDR ），然后生成克隆型（clonotype是一组测序reads相同的CDR3氨基酸或核苷酸序列或V / J基因）集，并用先进的算法的校正PCR和测序错误。

进一步对IG/TR的量化分析可以反应出机体的免疫状态，目前已有一些方法来进行评价，包括MiTCR Viewer和ViDJiL，今天我们为大家分享一种新的方法：tcR包，它整合了多种计算方法，可以对个体的免疫组库进行量化及比较分析，包括：基因usage的比较，共享clonotypes的检索，频谱分析，生成随机TR，多样性的评估以及其它常用的免疫组库分析方法。

营养补充

T细胞（抗原）受体（T cell receptor ,TCR）为所有T细胞表面的特征性标志，以非共价键与CD3结合，形成TCR—CD3复合物。TCR的作用是识别抗原。TCR分子属于免疫球蛋白超家族，其抗原特异性存在于V区；V区（Vα、Vβ）又各有三个高变区CDR1、CDR2、CDR3，其中以CDR3变异最大，直接决定了TCR的抗原结合特异性。TCR是由两条不同肽链构成的异二聚体，由α、β两条肽链组成，每条肽链又可分为可变区（V区），恒定区（C区），跨膜区和胞质区等几部分；其特点是胞质区很短。

在T细胞发育过程中CDR1、2和FR区域相对保守，CDR3区由V、D和J 进行重排而形成具有功能的TCR编码基因（T细胞克隆），由于V（65~100种）、D（2种）、J（13种）基因片段本身具有多样性，此外，由于在重排的过程中，在VD及D-J的连接区经常有非模板的核苷酸的随机插入或删除，进一步增加了CDR3区的多样性。这种基因片段连接的不准确性使TCR的表达呈多样性，以识别各种不同的抗原。

R包使用

install.packages("tcR") #安装R包library(tcR) #加载

一、R包中示例数据

1. “twinsdata”数据集

包含twa.rda和twb.rda这两个列表数据，twa.rda和twb.rda分别包含4 个数据框，每个数据框10000行。每个数据框都是双胞中的一个样本降采样(downsampled，目的是生成缩略图)到10000最丰富的克隆型（alpha和beta链）的数据。

twa.rda的第一个数据框数据及解释：

1) Umi.count：barcodes的数量

2) Umi.proportion：barcodes的比例。

3) Read.count：reads的数量

4) Read.proportion：reads的比例

5) CDR3.nucleotide.sequence：CDR3核苷酸序列

6) CDR3.amino.acid.sequence：CDR3氨基酸序列

7) V.gene：比对的可变基因片段的名称

8) J.gene：比对的加入基因片段的名称

9) D.gene：比对的多种基因片段的名称

10) V.end：比对的V基因片段的最后位置(1-based，碱基的起始读数为1)

11) J.start：比对的J基因片段的开始位置(1-based)

12) D5.end：比对的D基因片段的D’5结尾位置(1-based)

13) D3.end：比对的D基因片段的D’3结尾位置(1-based)

14) VD.insertions：V-D接合处插入核苷酸的数量 (N-nucleotides)

15) DJ.insertions：D-J接合处插入核苷酸的数量 (N-nucleotides)

16) Total.insertions：总的插入核苷酸的数量

2.“genesegments”数据

genesegments是由个数据框组成的列表，每个数据框是人类alpha-beta链片段数据，

genesegments的第一个数据框数据及解释：

1) V.allelles / J.allelles / D.allelles ：V/D/J片段名字

2) CDR3.position：CDR3起始序列所在位置

3) Full.nucleotide.sequence：CDR1-2-3序列片段

4) Nucleotide.sequence ：CDR3序列片段

5) Nucleotide.sequence.P：有 P-insertions的CDR3序列

注：tcR所有字符串都属于“character”类，而不是“factor”类。因此，应该使用带有参数stringsAsFactors=FALSE的R解析函数。

二、描述性统计量

1. 克隆集汇总Cloneset summary

（1）cloneset.stats()

cloneset.stats()用于获取一个整体视图（序列的总体计数、框内和框外的数量和比例等）。它返回核苷酸和氨基酸克隆型计数，以及读计数的汇总:

cloneset.stats(twb)

（2）

repseq.stats(twb)

2. 高丰度的克隆型统计Most abundant clonotypes statistics

（1）clonal.proportion()

该函数是通过读取克隆类型的计数来获得最丰富的数目。例如，计算占总“Read.count”的25%(大约)的克隆型的数量。举例如下：

clonal.proportion(twb, 25)

（2）top.proportion()

要获得最丰富的clonotypes的reads总和所占一个集合中reads总数的比例，可以使用top，即top克隆型reads/所有克隆型reads。

例：计算 top10的clonotypes的reads总和所占一个集合中reads总数的比例

top.proportion(twb, 10)

（3）vis.top.proportions()

vis.top.proportions()用来可视化（2）中计算的比例

举例：

vis.top.proportions(twb)

（3）tailbound.proportion()

该函数使用.col和.bound得到具有列.col的值≤.bound的特点的clonotypes的子集，并计算这种子集的 reads和占整个数据框的比例。

举例：得到“Read.count”<= 100的 reads总和占总的reads的比例

tailbound.proportion(twb, 100)

3. 克隆空间内稳态Clonal space homeostasis

克隆空间内稳态显示了克隆型按特定比例占据了多少空间。

（1）例：计算按照 Low (0 < X <= 0.05) High (0.05 < X <= 1)的比例所占空间

clonal.space.homeostasis(twb, c(Low = .05, High = 1))

（2）例：使用默认参数（会把比例划分成5个不同区域）

clonal.space.homeostasis(twb[[1]])

（3）可视化：

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

4. 框内框外序列In-frame and out-of-frame sequences

执行对in-frame和out-of-frame的clonotypes计数的函数是count.inframes、count.outframes、 get.inframes和get.outframes。该函数的参数.head用于输入数据框或子设置之前的数据框的输入列表。该函数接受数据框和数据列表作为参数。

（1）举例：获取只有in-frame序列的数据框，并在该数据框的前5000行中计算out-of-frame序列。

imm.in <- get.inframes(twb)  #获取twb列表的in-frame序列count.outframes(twb, 5000)  #计算前5000行中out-of-frame序列

（2）get.frames和count.frames函数中，参数“all”代表所有序列，参数“in”代表只看in-frame序列，参数“out”代表只看out-of-frame序列.

举例：

imm.in <- get.frames(twb, 'in') # 等同于'get.inframes(twb)'，获取框内序列count.frames(twb[[1]], 'all')   # 所有序列数，即返回行数

flag <- 'out'count.frames(twb, flag, 5000)   # 等同于 'count.outframes(twb, 5000)'，计算前5000行中out-of-frame序列

5. 检索靶向CDR3序列Search for a target CDR3 sequences

使用find.clonotypes函数对序列进行的精确或模糊搜索。该函数输入参数是数据框或数据列表，目标（是有一列是序列和其他附加列的向量或数据框），一列或多列的返回值，比较两个序列(精确匹配用“exact”；用Hamming距离匹配序列用“hamm”(即当H≤1时2个序列匹配；用Levenshtein距离匹配序列用“lev”(即当L≤1时2个序列匹配)的方法，匹配的序列的列名。

首先检索 CDR3序列

cmv<-data.frame(CDR3.amino.acid.sequence=c('CASSSANYGYTF',                'CSVGRAQNEQFF', 'CASSLTGNTEAFF',                'CASSALGGAGTGELFF', 'CASSLIGVSSYNEQFF'),                V.genes = c('TRBV4-1', 'TRBV4-1',                            'TRBV4-1', 'TRBV4-1', 'TRBV4-1'),                stringsAsFactors = F)cmv

下面使用所有匹配方法（exact, hamming or levenshtein）来进行搜索匹配或未匹配V-segment（V基因体片段是免疫球蛋白或T细胞受体基因中的一种DNA序列，因胚系基因组中有多个不同的V基因体片段而呈现变异性。一个V基因体片段和一个J或DJ基因片段发生重排后，形成一个编码V结构域的外显子）。

（1）例一’exact’

twb <- set.rank(twb)#对twb创建Rank"和"Index"列cmv.imm.ex <-  find.clonotypes(.data = twb[1:2], .targets = cmv[,1],                  #选取twb数据的前两个数据框                  #目标基因为上述cmv                  .method = 'exact',                  #匹配方法为“exact”                  .col.name = c('Read.count', 'Total.insertions'),                  #要输出的列                  .verbose = F)                  #是否输出错误信息head(cmv.imm.ex)

（2）例二hamming距离

cmv.imm.hamm.v <-   find.clonotypes(twb[1:3], cmv,                  'hamm', 'Rank',                  #使用hamming距离进行匹配                  #输出“Rank”，数据框中的克隆或克隆类型的Rank                  .target.col = c('CDR3.amino.acid.sequence',                                  'V.gene'),                  #用于计算的列                  .verbose = F)head(cmv.imm.hamm.v)

（3）例三levenshtein距离

cmv.imm.lev.v <-   find.clonotypes(twb[1:3], cmv, 'lev',               #使用levenshtein距离进行匹配                  .target.col = c('CDR3.amino.acid.sequence', 'V.gene'),                  .verbose = F)head(cmv.imm.lev.v)

三、基因usage

V区和J区基因的usage（V-usage and J-usage）是免疫组库的重要特征。有一些函数可以比较 tcR之间的基因取用情况。

1. 基因usage计算Gene usage computing

使用geneUsage函数评估 tcR的基因usage情况，输入数据框或列表，计算其给定元素（如V genes）的频率或计数。人类TCR和Ig的V和J基因名存储在.rda文件genesegments.rda中。函数的输出是数据框，第一列表示一个基因，另一列表示频率。

（1）例：

imm1.vs <- geneUsage(.data=twb[[1]], .genes=HUMAN_TRBV)#输入数据为twb，计算基因为HUMAN_TRBVhead(imm1.vs)

（2）

imm.vs.all <- geneUsage(twb, HUMAN_TRBV)#输入数据为twbimm.vs.all[1:10, 1:4]

（3）计算V-J基因联合计数

imm1.vj <- geneUsage(twb[[1]], list(HUMAN_TRBV, HUMAN_TRBJ))imm1.vj[1:5, 1:5]

（4）函数vis.gene.usage可视化

①

vis.gene.usage(.data=twb, .genes=HUMAN_TRBJ,                .main = 'twb J-usage dodge',               .dodge = T)#.dodge = T所有结果在一个条形图展示

②

vis.gene.usage(.data=twb, .genes=HUMAN_TRBJ,                .main = 'twb J-usage dodge',               .dodge = F, .ncol = 2)#.dodge = F是每个数据框的结果形成单独的条形图展示

③展示twb第一个数据框中，基因HUMAN_TRBV的频率

vis.gene.usage(imm1.vs, NA,                .main = 'twb[[1]] V-usage',               .coord.flip = F)#.coord.flip是否翻转坐标

注：这个R包还没有讲解完哦还会陆续更新~

引用文献

Nazarov VI, Pogorelyy MV, Komech EA, et al. tcR: an R package for T cell receptor repertoire advanced data analysis. BMC Bioinformatics. 2015;16(1):175. Published 2015 May 28. doi:10.1186/s12859-015-0613-1