免疫组库数据分析（二）：Excel 分析免疫组库数据

免疫组库数据分析（二）：Excel 分析免疫组库数据

前言

在系列文章第一篇《免疫组库数据分析（一）：windows 系统下MiXCR的安装和使用》讲解了5’RACE实验数据如何在Window系统下进行初步的比对拼接与输出文件操作。本篇将讲解如何利用Excel对CDR3克隆子输出文件进行分析。按照数据清洗、数据筛选及计算分析方面进行如下分析

1. 选择目的TCR链序列
2. 去除无效CDR3序列
3. 重新计算每一个克隆的所占比例
4. 计算克隆子中V、J基因的使用频率（Clone counts水平与Clone types水平）
5. 计算VJ 组合的使用频率（Clone counts水平与Clone types水平）
6. 计算免疫组库CDR3的氨基酸长度分布
7. 计算免疫组库CDR3的20种氨基酸使用频率
8. 计算免疫组库的CDR3的多样性

数据来源

本篇数据选用小鼠5‘RACE TCRB 免疫组库数据，数据文本格式为txt格式。直接用excel 表格打开

参数：

CloneID：每一个克隆子的ID号，按照Clone Count的数量从大往下排列

CloneCunts:每一个克隆子的couts数，可理解为Reads数。

cloneFraction: 每一个克隆子所占的比例

TargetSequnces:每一个克隆子的所检测到的DNA序列.

all V hits withScore:每一个克隆子的比对到的不同V基因及其分数

allD HitsWithScore：每一个克隆子的比对到的不同D基因及其分数

all J HitsWithScore：每一个克隆子的比对到的不同J基因及其分数

all C HitsWithScore：每一个克隆子的比对到的不同J基因及其分数

aaSeqCDR3：基于 CDR3 DNA序列翻译的氨基酸序列中。翻译从两端向中心翻译，每3个碱基翻译成一个氨基酸。*为终止密码子，‘’—""为DNA序列无法翻译成氨基酸序列.凡是有以上两个情况的均无法产生有效的CDR3氨基酸序列。

数据清洗

1. 筛选得分最高的V，D，J，C 基因序列：用Excel 表格的分列进行操作。

备注：在mixcr exportClones [options] clones.clns clones.txt可选参数中，可以选择-vAlignment、-vGene 等参数直接分选出得分最高的V基因序列。

所有数据复制到新的电子表格2中，首先all V hits withScore列数据复制到新的电子表格3中，点击Excel 窗口中点击-【数据】-【分列】-【分隔符号】-【其他】,输入星号*-点击【下一步】-【完成】即可。将所获得的新列数据直接覆盖电子表格中。

all D HitsWithScore、all J HitsWithScore、all J HitsWithScore、all C HitsWithScore 按照上述步骤进行相同操作。

数据筛选

1、选择目的TCR链序列：Excel 筛选 特定J基因序列确定目的 TCR链

在5’RACE数据中，往往存在非特异性比对到其他TCR链或者BCR链。因此需要通过 J 基因来筛选目的链所有的克隆子序列。本篇中，选取all J HitsWithScore列已经清洗过的数据，点击 【数据】-【筛选】在搜索框内输入"TRBJ",即可获得含目的TRB链的克隆子，将筛选后的数据复制到新的电子表格中。

2、去除无效CDR3序列：Excel 筛选 有效CDR3序列(In Frame sequnces)

​ aaSeqCDR3的氨基酸序列是根据TargetSequences 列中DNA序列翻译而来，从两端往中间翻译，翻译出来的氨基酸序列中*星号说明是终止密码子，“-”说明是对应的DNA序列（只剩1-2个碱基）无法翻译成氨基酸序列。此两者情况均为无效CDR3 克隆子（Out of Frame sequnces).

​ 点击 【数据】-【筛选】-【自定义筛选】-选择【不包含】选项输入“-” 与“~*”，因为 ”*"是通配符，前面加入波浪符“~”进行转义即可。 aaSeqCDR3列就能筛选出有效CDR3序列。筛选完后，复制数据至新的表格中，粘贴格式是仅为数值。

3：重新计算每一个克隆子的所占比例：SUM函数与"/"除法

将cloneFraction 列数据清空，在C2单元格，输入=B2/SUM(B:B),然后将鼠标放在单元格右下角，出现实心十字时候，双击左键，即可按照如此公式自动往下填充。

4、计算克隆子中V、J基因的使用频率（Clone counts水平）:使用SUMIF 函数

Clonecounts:指一个克隆子在整个免疫组库中的绝对量。特别是存在特异克隆扩增时候，Clonecounts尤为重要。

小鼠TCRB链V基因中选取选取有功能的V基因列进行统计，不统计假基因使用频率。文献参考[1]。使用SUMIF函数如下

SUMIF(Range,criteria,[sum_range])
参数解释：
range：需要查询的范围；
criteria：查询参数值
sum_range：求和统计范围

举例说明：M2=SUMIF($E:$E,$L2,$C:$C)

M2单元格为SUMIF条件求和值，range：查询范围为$E:$E，即E列全列，并用$符号进行绝对引用。criteria：查询参数值为 $L2 ,为相对引用，只在L列前加入$,行名不加。说明该引用无论向下自动填充公式时，只会行变，列不变，目的就是正确自动填入查询参考值。sum_range：求和统计范围，在E列查找到L2上的值（TRBV1）时，求和的范围在相同行的C列，并累加求和。

在M2单元格将鼠标放在单元格右下角，出现实心十字时候，双击左键，即可按照如此公式自动往下填充,直到遇到空白单元格为止。

如上图所示：有功能的V基因使用频率总和＜1，说明存在少量假基因未被计入统计。

克隆子中有功能的J基因的使用频率计算方式相同。

5、计算克隆子中V、J基因的使用频率（Clonetype水平）:使用CountIF 函数

Clonetype:指克隆子无论其扩增多少，有多少clone count数量，一种克隆只算一种克隆类型。使用COUNTIF函数计算目的基因出现的次数。

COUNTIF(ragne,criteria)
参数解释：
range：需要查询的范围；
criteria：查询参数值

举例说明： K2=COUNTIF($G:$G,$J2)

K2单元格为COUNTIF条件求次数值，range：查询范围为$G:$G，即G列全列，并用$符号进行绝对引用。criteria：查询参数值为 $J2,为相对引用。将鼠标放在单元格右下角，出现实心十字时候，双击左键，即可按照如此公式自动往下填充,直到遇到空白单元格为止。

如下图，TCRB的功能J基因共计12个，通过COUNTIF计算Clonetype水平中每一个J基因使用次数，并计算其所占百分比。

6、计算VJ 组合的使用频率（Clone counts水平）

将clonecounts列，V基因，J基因列复制到新的表格中，将TRB链功能V基因与J基因 构建一个空白矩阵。

• 利用SUMIFS函数计算免疫组库中，克隆子 V基因与J基因组合使用的counts总量。

SUMIFS(sum_range,ciriteris_range1,criteria1,[ciriteris_range2,criteria2],...)
参数解释：
sum_range:求和范围
ciriteris_range1：参数1的查询的范围；
criteria1：查询参数1的值
ciriteris_range2：参数2的查询的范围；
criteria2：查询参数2的值

举例说明：E2=SUMIFS($A:$A,$B:$B,$D3,$C:$C,E$2) :计算TRBV1（参数1)-TRBJ1-1(参数2）组合的克隆子的count总数。

$A:$A:求和范围-clonecounts列

$B:$B:查询范围1-V基因列

$D3：参数值1-对应到D3单元格内容：TRBV1

$C:$C：查询范围2-J基因列

E$2：参数值2-对应到E2单元格内容：TRBJ1-1

需要理解$所代表的的引用是绝对引用。根据E2单元格公 式，向下向右拉，自动填充公式，计算数值，结果如下。

• 下一步计算矩阵数值总和，求出每一个VJ组合的百分比是多少：

1.复制V-J组合矩阵到原有矩阵的下面，清空数值，

2. E28单元格内输入=100*E3/SUM($E$3:$O$24)"，计算TRBV1-TRBJ-1组合counts数量占计算所得总体Counts的百分比。
3. SUM-V列，SUM-J 行分别计算了每一种功能V基因、J基因的使用频率。

7.计算VJ 组合的使用频率（Clone tpyes水平）

利用COUNTIFS函数计算Clone tpyes水平 克隆子V-J组合的使用次数。首先复制V-J组合矩阵到表格空白，清空数值。

COUNTIFS(criteria_ragne1,criteria1,[criteria_ragne2,criteria2],....)
参数解释：
criteria_range1：参数1所查询的范围；
criteria1：查询参数1的值
ciriteris_range2：参数2的查询的范围；
criteria2：查询参数2的值

举例说明：
’S5=COUNTIFS($B:$B,$R5,$C:$C,S$2) 单元格S4计算TRBV3（参数1）-TRBJ1-1（参数2）组合的次数。

接下来计算 V-J组合使用频率（百分比），方法同5，结果如下。

8.计算免疫组库CDR3的氨基酸长度分布：LEN函数与SUMIF函数

8.1 先计算每一个CDR3氨基酸的长度：LEN函数

C2=LEN(B2) 计算B2单元格中氨基酸的长度

8.2 计算CDR3氨基酸长度分布：SUMIF函数

先在空白列构建1-30的等差数列，在相邻行单元格利用SUMIF函数，写入对应公式。

举例：F15=SUMIF($C:$C,$E15,$A:$A)*100

说明：该单元格计算E15（CDR3 aa长度=14）的使用频率。查询范围是C列（氨基酸长度），参数值是E15（14），求和范围是A列（克隆子比例），结果如下。

9.计算免疫组库CDR3的20种氨基酸的使用频率

9.1 计算每一个克隆子CDR3中20种氨基酸的个数

SUBSTITUTE(text,old_text,new_text) 
参数：Text 为需要替换其中字符的文本。 
Old_text 为需要替换的旧文本。 
New_text 用于替换 old_text 的文本。

举例说明：H2=LEN（$B2)-LEN(SUBSTITUTE($B2,H$1,""))

SUBSTITUTE($B2,H$1," "):将单元格B2中CDR3氨基酸序列的氨基酸R(H1单元格内容提供)替换成空

LEN(SUBSTITUTE($B2,H$1,""))：将氨基酸R替换为空白后，CDR3的氨基酸长度

LEN（$B2)：B2中CDR3氨基酸长度

LEN（$B2)-LEN(SUBSTITUTE($B2,H$1,""))：B2中CDR3氨基酸序列中氨基酸R的个数

注意：公式往下往右拉，自动填充时，需要理解$符号的应用。

另外AB列利用SUM函数，验证该公式的正确性。CDR3氨基酸总个数前后一致。

9.2 计算免疫组库中，每一个克隆子CDR3中，单个氨基酸的使用频率。

如上图 所示，第二个克隆子CDR3氨基酸长度是12，氨基酸R 出现3次，因此氨基酸R在这条CDR3氨基酸中使用频率是3/12=0.25，而第二个克隆子的使用频率（Clone Fraction)为0.0005108，因此氨基酸R在CDR3整体免疫组库中使用频率为 0.0005108*0.25=0.0001277。

批量计算如下：

第一个克隆子CDR3中氨基酸R的使用频率：AD2单元格输入“=$A2\ast (H2/$C2)”

9.3 计算免疫组库中20种氨基酸的使用频率：

在右侧空白单元格处复制20种氨基酸名称，利用SUM函数，分别计算AD-AW列对应氨基酸的频率总和。

例如氨基酸R的使用频率总和,在单元格AY中输入“=SUM（AD：AD)"。然后右拉自动填充公式，计算其他19种氨基酸的使用频率总和。

10.计算免疫组库的CDR3的多样性

10.1 多样性指数介绍

1）香农（Shannon）指数。计算公式为：
$$H\prime =-\sum （pi㏑pi）$$
在式中，H′为香农指数；pi为第i个种在全体物种中的比例。如以个体数量而言，ni为第i个种的个体数量，N为总个体数量，则有pi=ni/N。㏑为自然对数，底数e=2.7182838。

2）辛普森指数（Simpson’s Index ）
它反映的是在同一个样本中随机的抽取2个个体，这两个个体来自同一个类的概率。
$$D=\sum (n/N{)}^2$$
n为第i个种的个体数量，N为总个体数量。 D值在0-1之间。0表示无限多样，1表示没有多样性。也就是说D值越大，多样性越低，由于与直觉相反，通常利用Simpson’s Reciprocal Index： 1 / D。即值越大，多样性越高

3）玛格列夫（Margalef）指数。仅考虑群落的物种数量和总个体数，将一定大小的样本中的物种数量定义为多样性指数。计算公式为：
$$D=（S-1）/㏑N$$
在式中，D为玛格列夫指数；N为总个体数量；S为总物种数量。

4）贝格-派克（Berger-Parker）指数。本指数仅考虑总种数及数量最优势种两个因素，易计算，侧重群落中优势种的作用。计算公式为：
$$d=Nmax/N$$
在式中，d为贝格-派克指数；Nmax为数量最优势种个体数量；N为总个体数。

10.2 举例说明

将原数据中cloneId列,cloneCounts 列;cloneFraction 列数据复制到新的表格中，其中cloneFraction为克隆子counts比例，等同于Pi。免疫组库CDR3多样性指数计算如下。