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