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

生物信息学

WGCNA:构建共表达网络的神器

2020-8-28 0:47:09

生物信息学

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

2020-8-28 0:49:53

加入Q群
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
今日签到
有新私信 私信列表
搜索