新一代测序(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_TRBV
head(imm1.vs)
(2)
imm.vs.all <- geneUsage(twb, HUMAN_TRBV)
#输入数据为twb
imm.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