导语
GUIDE ╲
免疫球蛋白(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 ）【了解CDR3重排请戳】，然后生成克隆型（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包使用
代码语言：javascript复制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()用于获取一个整体视图（序列的总体计数、框内和框外的数量和比例等）。它返回核苷酸和氨基酸克隆型计数，以及读计数的汇总:
代码语言：javascript复制cloneset.stats(twb)（2）
代码语言：javascript复制repseq.stats(twb)2. 高丰度的克隆型统计Most abundant clonotypes statistics
（1）clonal.proportion()
该函数是通过读取克隆类型的计数来获得最丰富的数目。例如，计算占总“Read.count”的25%(大约)的克隆型的数量。举例如下：
代码语言：javascript复制clonal.proportion(twb, 25)（2）top.proportion()
要获得最丰富的clonotypes的reads总和所占一个集合中reads总数的比例，可以使用top，即top克隆型reads/所有克隆型reads。
例：计算 top10的clonotypes的reads总和所占一个集合中reads总数的比例
代码语言：javascript复制top.proportion(twb, 10)（3）vis.top.proportions()
vis.top.proportions()用来可视化（2）中计算的比例
举例：
代码语言：javascript复制vis.top.proportions(twb)（3）tailbound.proportion()
该函数使用.col和.bound得到具有列.col的值≤.bound的特点的clonotypes的子集，并计算这种子集的 reads和占整个数据框的比例。
举例：得到“Read.count”<= 100的 reads总和占总的reads的比例
代码语言：javascript复制tailbound.proportion(twb, 100)3. 克隆空间内稳态Clonal space homeostasis
克隆空间内稳态显示了克隆型按特定比例占据了多少空间。
（1）例：计算按照 Low (0 < X <= 0.05) High (0.05 < X <= 1)的比例所占空间
代码语言：javascript复制clonal.space.homeostasis(twb, c(Low = .05, High = 1))（2）例：使用默认参数（会把比例划分成5个不同区域）
代码语言：javascript复制clonal.space.homeostasis(twb[[1]])（3）可视化：
代码语言：javascript复制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序列。
代码语言：javascript复制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序列.
举例：
代码语言：javascript复制imm.in <- get.frames(twb, 'in') # 等同于'get.inframes(twb)'，获取框内序列
count.frames(twb[[1]], 'all') # 所有序列数，即返回行数代码语言：javascript复制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序列
代码语言：javascript复制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'
代码语言：javascript复制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距离
代码语言：javascript复制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距离
代码语言：javascript复制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）例：
代码语言：javascript复制imm1.vs <- geneUsage(.data=twb[[1]], .genes=HUMAN_TRBV)
#输入数据为twb，计算基因为HUMAN_TRBV
head(imm1.vs)（2）
代码语言：javascript复制imm.vs.all <- geneUsage(twb, HUMAN_TRBV)
#输入数据为twb
imm.vs.all[1:10, 1:4]（3）计算V-J基因联合计数
代码语言：javascript复制imm1.vj <- geneUsage(twb[[1]], list(HUMAN_TRBV, HUMAN_TRBJ))
imm1.vj[1:5, 1:5]（4）函数vis.gene.usage可视化
①
代码语言：javascript复制vis.gene.usage(.data=twb, .genes=HUMAN_TRBJ,
.main = 'twb J-usage dodge',
.dodge = T)
#.dodge = T所有结果在一个条形图展示②
代码语言：javascript复制vis.gene.usage(.data=twb, .genes=HUMAN_TRBJ,
.main = 'twb J-usage dodge',
.dodge = F, .ncol = 2)
#.dodge = F是每个数据框的结果形成单独的条形图展示③展示twb第一个数据框中，基因HUMAN_TRBV的频率
代码语言：javascript复制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
END