10X Genomics单细胞免疫组库VDJ分析必知必会

分享是一种态度

作者 | 周运来

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树，

一枚大型测序工厂的螺丝钉，

一个随机森林中提灯觅食的津门旅客。

我们生活着的世界并非只有我们自己，而是有很多小于或大于我们的生物不断与我们交互着，有的让我们开心，有的使我们伤心。这就关系到一个本质的问题：

我是谁？

我是一个脑袋，一个躯体，一个想法。。。？

自我是指：机体胚系基因（germ-line-gene）的编码产物，在免疫学上还要加上一条：免疫系统发育早期遭遇过的物质。

免疫学的功能就是识别自我与非我，并保持自我的稳态的一门学问。那么自我的复杂度要远小于非我，有限的自我如何识别无限的非我呢？我们需要不同的屏障和内在的库存，这里的库有一部分就是我们今天说的免疫组库。

什么是免疫组库

或谈淋巴细胞抗原受体多样性的产生。

参与适应性免疫的淋巴细胞与参与固有免疫的非淋巴细胞间的一个重要区别就是，淋巴细胞具有结构多样的抗原受体储备。我们关心的是这个库是如何产生的。这就不得不看看淋巴细胞的发育过程了:

引自免疫学原理

其实B细胞和T细胞的结构还是比较相似的。在B细胞和T细胞中发现的受体分别被称为B细胞受体（B Cell Receptor ，BCR）和T细胞受体（T Cell Receptor，TCR）。抗原的检测过程因白细胞（淋巴细胞（lymphocyte）是白细胞的一种，是体积最小的白细胞）类型的不同而不同，如B细胞或T细胞。B细胞受体与自由存在的可溶性抗原结合，而T细胞受体仅在主要组织相容性复合体(MHC)上识别抗原。这是B细胞受体和T细胞受体的关键区别。

当我们把BCR/TCR的结构放大之后是这样的：

以B细胞抗原受体（BCR）为例，它包含了两条重链（H）和两条轻链（L）。重链包含了一个可变区（VH）和3个恒定区（CH1/CH2/CH3），轻链则包含了一个可变区（VL）和一个恒定区（CL）。而免疫球蛋白有IgA、IgG、IgM、IgD、IgE五种，各自可以搭配κ或λ两种轻链。免疫学的书上说，IG的恒定区决定了它的免疫原性，即重链的恒定区决定了它是5种球蛋白中的哪一种，轻链的恒定区则决定了搭配重链的是两种轻链中的哪一种；而IG的可变区则决定了它的特异性（克隆性），即决定了它和什么样的抗原结合。

引自免疫学原理

在可变区里面，还有变化相对小的4个骨架区（FR1/FR2/FR3/FR4），变化相对大的3个互补决定区（CDR1/CDR2/CDR3），其中变化最大的是CDR3区。

互补决定区是VL与VH均有3个HVR，它们共同组成Ig的抗原结合部位（antigen-binding site）。该部位因在空间结构上可与抗原决定簇形成精密的互补，故超变区又称互补性决定区(complementarity determining region，CDR)。整个抗体分子可分为恒定区和可变区两部分。在可变区内有一小部分氨基酸残基变化特别强烈，这些氨基酸的残基组成和排列顺序更易发生变异区域称高变区。在L链、H链的V区中有三个高变区(hypervariable regions，HVR)，该部位因在空间结构上可与抗原决定簇形成精密的互补，故高变区又称互补性决定区。互补决定区||百科(https://baike.baidu.com/item/%E4%BA%92%E8%A1%A5%E5%86%B3%E5%AE%9A%E5%8C%BA/8733130)

image

淋巴细胞抗原受体多样性的产生归根到底还是胚系基因的重组和重排。人体TCR和BCR胚系基因的组成：

引自免疫学原理

胚系基因的特点和基因转录：

• 四种人体受体链编码基因分属不同染色体上五个不同的基因座位，比如编码免疫球蛋白IG三种链（IgH，Igκ，Igλ）的基因，并不在同一条染色体上，编码重链的基因位于14号染色体长臂，编码轻链κ的基因位于2号染色体短臂，编码轻链λ的基因位于22号染色体长臂。
• 除了C区，多数VDJ区由多个片段组成。
• 在DNA水平，不同基因区段会相互组合。
• VDJ基因片段重组构成淋巴细胞受体结构的多样性。

那么这种组合能产生怎样规模的多样性呢？

重排过程选择的片段是随机的， 随机重排的结果可以产生108-1010种免疫球蛋白分子。其次, 淋巴细胞在成熟后还会发生体细胞突变, 又极大的增加了抗体的多样性. 综合来看, H链和L链可能的组合数预计可达1011以上。小鼠每天生成的淋巴细胞为108个，可见动物一生中还不可能使全部可能产生的基因组合得到表达。

基因重排发生在B细胞成熟前，体细胞高频突变发生在B细胞成熟后。B细胞从干细胞发育成成熟的B细胞之前，要经历重链的VDJ基因重排和轻链的VJ基因重排；成熟的B细胞，迁移到外周淋巴器官之后，在外界抗原的刺激下，可能发生体细胞高频突变（SHM）。

也许你会问(尽管很可能不会)，为什么B细胞发育成熟前会产生这种基因重排而产生结构异相的各种受体呢？这是因为只有发育中的淋巴细胞可选择表达重组激活基因（recombinase activating gene，RAG）。

人重组激活基因是免疫球蛋白（immunoglobulin，Ig）和T淋巴细胞受体（T lymphocyte receptor，TCR）基因片段重排所必需，RAG1基因突变使其编码的重组酶活性完全或部分丧失，V(D)J重组失衡，T淋巴细胞和B淋巴细胞的发育在早期被阻断，导致原发性免疫缺陷病（primary immunodeficiency, PID）,为常染色体隐性遗传。

至此，我们知道了免疫组库的来源以及基本组成。

那么，

为什么要做单细胞免疫组库

当然是为了进一步探究人体免疫机制，挖掘免疫组库与疾病的关系，促进人类健康了。免疫细胞是如何产生作用的呢？

鉴于抗原受体的多样性，不在单细胞水平上的研究只能得到模糊的宏观视角。其实单细胞水平的研究，开始的很早，只是通量很低。

• 揭示克隆性、多样性、抗原特异性和细胞环境
• 组装并注释全长V(D)J基因序列
• 从单个T细胞识别α和β链序列
• 将来自单个B细胞的重链和轻链免疫球蛋白(Ig)序列以全同型分辨率配对
• 同时测定同一细胞中TCR、B细胞Ig、细胞表面蛋白表达及5’基因表达
• 配对TCRα和β链与TCR-pMHC特异性序列
• 同时测定细胞表面蛋白和基因表达

主要应用：

• 免疫组库可以捕捉肿瘤发生时免疫微环境的变化，寻找免疫治疗的靶点，从而辅助免疫治疗更好的抗击肿瘤。
• 器官或者骨髓移植时，经常会诱发宿主排斥反应的发生，从而发生慢性移植抗宿主病。
• 自身免疫免疫性疾病是由于机体对自身抗原发生免疫反应而导致自身组织损害所引起的疾病。
• 免疫组库在感染性疾病、抗体开发、用药及疫苗评估等多个方面均有应用价值。例如通过免疫组库研究，可以检测感染类疾病过程中的免疫动态变化；在抗体开发方面，可以获得特征性的BCR（Ig）序列，缩短抗体开发的流程；也可以针对某种疾病用药后的外周血样本进行评估，确认药物是否激发免疫反应及其功效。

单细胞免疫组库如何做

一般的免疫组库测序（Immune Repertoire sequencing(IR-SEQ)）以多重PCR或5’RACE技术目的扩增决定B细胞受体（BCR）或T细胞受体（TCR）多样性的互补决定区（CDR区），再结合高通量测序技术，全面评估免疫系统的多样性，深入挖掘免疫组库与疾病的关系。

用淋巴细胞分离液分离外周血T/B淋巴细胞，提取DNA（或RNA），采用多重PCR/5'RACE对CDR3进行捕获（5'RACE还可以测CDR1、CDR2），通过 Hiseq2000（Hiseq2500、Miseq）平台进行高通量测序。

DNA水平选择多重PCR方法，侧重研究基因重组信息；RNA 水平可选择多重PCR或5'RACE方法，侧重于研究基因的表达状态。

但是这种方法并不能获得每个细胞BCR和TCR的具体状态。2015年，10× Genomics发布了基于微流控和油滴包裹技术的Chromium单细胞系统平台，可实现高通量的单细胞转录组和单细胞V(D)J测序。不但可以将TCR/BCR双链完美匹配，而且可以细化到单细胞水平，同时获得表达谱信息。

image

我们这里主要介绍这款仪器：

image

10× Genomic单细胞免疫组库测序是建立在GemCode技术上的微流体平台，将带有条形码和引物的凝胶珠与单个细胞包裹在油滴中；接下来在每个油滴内，凝胶珠溶解，细胞裂解释放mRNA，通过逆转录产生用于测序的带条形码的cDNA。液体油层破坏后，cDNA一分为二，后续同时进行基因表达和免疫组库文库构建；其中TCR或者BCR的V(D)J序列通过设计在TCR或者BCR、lg的C区域的巢式PCR引物进行富集。然后使用Illumina测序平台对文库进行测序检测，即可一次性获得大量单细胞的基因表达和免疫组库数据，实现在单细胞水平同时对基因表达和免疫组库进行研究。

熟悉10X单细胞转录组的朋友对这一套流程绝对不会陌生：捕获 ，建库测序，拆库定量，数据分析：

image

区别在于在单细胞水平上识别并扩增VDJ区域的基因。

10× Genomics提供了完整的实验流程和数据分析方案。实验方面更多的是经验的积累，这里我们还是关心一下拿到reads之后的生物信息学分析吧。其实10X单细胞的生信入门的门槛是很低的，大部分工作都被它的cellranger做完了。干这一行的，以至于有的生信工程师都没见过fastq序列长什么样，拿到序列直接朝cellranger里面一丢就可以得到几乎全套的结果了。这也反反映出作为下游的生信工程师应该注意修炼的数据挖掘的功底。

因为只是配参数和投任务的不叫工程师：

$ cd /home/jdoe/runs
$ cellranger vdj --id=sample345 \
                 --reference=/opt/refdata-cellranger-vdj-GRCh38-alts-ensembl-3.1.0 \
                 --fastqs=/home/jdoe/runs/HAWT7ADXX/outs/fastq_path \
                 --sample=mysample \
                 --localcores=8 \
                 --localmem=64 

我们肯定想知道cellranger vdj执行的一般过程以及关键节点：

image

其实我们完全是有必要cd到cellranger路径下读一读源码的，不就是python代码吗？可惜似乎永远有80%的生信工程师处于入门阶段，自己的python还没整明白呢。于是很长一段时间，我们都在纳闷这些结果是怎么产生的：

Outputs:
- Run summary HTML:                                  /home/jdoe/runs/sample345/outs/web_summary.html
- Run summary CSV:                                   /home/jdoe/runs/sample345/outs/metrics_summary.csv
- All-contig FASTA:                                  /home/jdoe/runs/sample345/outs/all_contig.fasta
- All-contig FASTA index:                            /home/jdoe/runs/sample345/outs/all_contig.fasta.fai
- All-contig FASTQ:                                  /home/jdoe/runs/sample345/outs/all_contig.fastq
- Read-contig alignments:                            /home/jdoe/runs/sample345/outs/all_contig.bam
- Read-contig alignment index:                       /home/jdoe/runs/sample345/outs/all_contig.bam.bai
- All contig annotations (JSON):                     /home/jdoe/runs/sample345/outs/all_contig_annotations.json
- All contig annotations (BED):                      /home/jdoe/runs/sample345/outs/all_contig_annotations.bed
- All contig annotations (CSV):                      /home/jdoe/runs/sample345/outs/all_contig_annotations.csv
- Filtered contig sequences FASTA:                   /home/jdoe/runs/sample345/outs/filtered_contig.fasta
- Filtered contig sequences FASTQ:                   /home/jdoe/runs/sample345/outs/filtered_contig.fastq
- Filtered contigs (CSV):                            /home/jdoe/runs/sample345/outs/filtered_contig_annotations.csv
- Clonotype consensus FASTA:                         /home/jdoe/runs/sample345/outs/consensus.fasta
- Clonotype consensus FASTA index:                   /home/jdoe/runs/sample345/outs/consensus.fasta.fai
- Clonotype consensus FASTQ:                         /home/jdoe/runs/sample345/outs/consensus.fastq
- Concatenated reference sequences:                  /home/jdoe/runs/sample345/outs/concat_ref.fasta
- Concatenated reference index:                      /home/jdoe/runs/sample345/outs/concat_ref.fasta.fai
- Contig-consensus alignments:                       /home/jdoe/runs/sample345/outs/consensus.bam
- Contig-consensus alignment index:                  /home/jdoe/runs/sample345/outs/consensus.bam.bai
- Contig-reference alignments:                       /home/jdoe/runs/sample345/outs/concat_ref.bam
- Contig-reference alignment index:                  /home/jdoe/runs/sample345/outs/concat_ref.bam.bai
- Clonotype consensus annotations (JSON):            /home/jdoe/runs/sample345/outs/consensus_annotations.json
- Clonotype consensus annotations (CSV):             /home/jdoe/runs/sample345/outs/consensus_annotations.csv
- Clonotype info:                                    /home/jdoe/runs/sample345/outs/clonotypes.csv
- Barcodes that are declared to be targeted cells:   /home/jdoe/runs/sample345/out/cell_barcodes.json
- Loupe V(D)J Browser file:                          /home/jdoe/runs/sample345/outs/vloupe.vloupe

Pipestance completed successfully!

首先我们可以看到的就是web_summary.html，这里面包含了我们样本的基本信息，也就是后期数据挖掘的前提。下面是TCR的结果，能看出TCRα chain or “TRA,” and TCRβ chain or “TRB”的cell占比：

image

image

然后是BCR的结果，可以看出IgH，Igκ，Igλ的注释结果：

image

image

算法概述

image

上图显示了10x V(D)Jread-pairs aligned到一个组装的contig，说明了read的结构。每个V(D)J链捕获1到几个umi。一轮以C-region5 '端为靶点的富集PCR反应，接着是酶的裂解反应，结果产生了来自同一转录本的分子池（pool of molecules ）。分子携带相同的10x条形码和UMI序列，但插入长度不同，导致R2起始点不同。R2起始点的多样性使每个转录本的目标部分得到完全覆盖，一般约为650bp。

组装过程将单个条形码的reads作为输入，并将这些reads“粘接”(glues)在一起，产生一组contig作为输出，这些contig代表对当前转录序列的最佳估计。此外，每一群中的每个基地都被赋予了一个质量值。我们还跟踪的UMIs的数量和reads支持的contig。

由于数据中存在多种形式的“噪声”，使得生成contig的问题十分复杂。这些原因包括背景(细胞外)mRNA、细胞双重态、细胞内转录错误、cDNA反转录错误、测序过程中的随机错误、测序过程中的指标跳变等。

组装过程在某些地方使用参考序列，如下所述，除非管道在denovo模式下运行。参见注释算法，其中部分内容在组装算法中使用。

组装算法的步骤如下:

可以看出使用的方法是De Bruijn graph结构来组装，当年学的宏基因组组装：从what 到how(https://www.jianshu.com/p/034ccfe097b5)派上用场了。

Targeted Cell Calling Algorithm

在10x系统中，液滴(GEMs)的数量很多，其中一些液滴含有一个细胞，另一些液滴含有一个靶细胞(T或B)。

目标细胞的检测依赖于其V(D)J转录本的鉴定和计数。一些T和B细胞对这些转录本的表达水平很低，因此可能无法检测到。相反，细胞外足够高水平的mRNA可能导致一些条形码被错误地识别为目标细胞。因此，目标细胞调用算法的目标是近似包含目标细胞的一组条形码。

该算法作为汇编算法的一部分执行。要被识别为目标细胞，条形码必须满足以下三个要求:

• 必须有一个多样的contig，如果只有一个这样的contig，必须有一个以上的UMI支持它的连接区域。(在denovo的情况下，我们只要求有一群。)虽然其他类型的细胞也能在TCR和BCR位点上转录，但只有T和B细胞能产生包含V和C片段的完全重新排列的转录本。因此，有一个生产的contig是很好的证据，从一个目标细胞转录存在于GEMs.然而，也有可能转录本是背景的——存在于细胞之间的液体中，而不是在一个完整的细胞中。对于这种情况，需要多个UMI可以提供一些x信息支持。
• 必须至少有三个过滤的umi，每个umi至少有两个read pairs(参见组装算法)。这降低了仅仅基于背景转录本来调用目标细胞的可能性。
• 计算所有条形码上每个UMI的读对数的N50值。如果对于给定的条形码，经过过滤的UMIs的最大读对计数小于N50的3%，则不要将条形码称为cell。这提供了一些保护措施，防止Illumina流式细胞仪上的索引跳变和其他形式的交叉文库污染引起的转录本。

除了以上列出的三个要求外，Cell Ranger 3.1还引入了一个新的过滤器来处理由浆细胞和含有大量RNA的B细胞引入的噪声(如Cell Ranger 3.1发布说明中记录的那样)。1)对与高频率或大型克隆共享一条链的低频率克隆收紧is_cell过滤器，2)缩小高频率克隆，以消除样品处理(例如，并非由于真正的生物克隆扩展)造成的mRNA泄漏带来的噪音。

Annotation Algorithm

V(D)J contig注释的目的是定义V、D、J片段对一个contig的比对，识别CDR3序列，从这些数据判断一个contig是否具有生产性，这意味着它可能对应于一个功能T或B细胞受体。

对于给定的数据集，管道首先确定数据是TCR还是BCR，然后相应地将所有的contigs对齐到TCR或BCR引用序列。在罕见的(混合的)情况下，contig都是对齐的。在12-mer的完美匹配上seeded对齐，然后进行启发式扩展;我们还从C段比对中反向搜索J段比对中不存在12-mer完全匹配的情况，因为这些情况偶尔会出现在体细胞超突变中。

重要的是要理解在比对中V(D)J参考序列的选择可以是任意的，这取决于参考序列彼此之间的相似程度。对于既短又突变较多的D段，通常不可能找到可靠的比对，而且可能没有显示比对。

如果满足下列条件，该条件被称为“有生产力productive ”:

• 完整的长度要求。重叠部分与V基因的起始部分匹配。该基因继续延伸，最终与J基因的末端相匹配。
• 起始要求。V的起始部分匹配contig上的起始密码子。注意，在10x提供的人类和小鼠参考序列中，每个V段都以一个起始密码子开始。
• 连续性。在V开始和J结束之间没有终止密码子。
• 位置。J停止减去V开始等于1模3。这就是说V和J段上的密码子在坐标系中。
• CDR3上要求。有一个带注释的CDR3序列(见下面)。 +结构要求。设VJ为V段和J段长度之和。让len表示J停止减去V开始，在contig上测量，那么VJ - len在-25和+25之间，除了IGH，它必须在-55和+25之间。这个条件是为了防止不可能与功能蛋白相对应的异常结构变化。

对于每一个contig，我们利用CDR3s的侧翼序列是保守的这一事实来搜索CDR3序列。我们将来自V和J参考片段的motifs与人类和小鼠进行比较，如下图所示。这里一个字母代表一种特定的氨基酸，一个点代表任何氨基酸。

left flank   CDR3   right flank
LQPEDSAVYY   C...   LTFG.GTRVTV
VEASQTGTYF          LIWG.GSKLSI
ATSGQASLYL

我们要求CDR3序列长度在5到27个氨基酸之间，以C开头，不包含终止密码子。候选CDR3的侧翼序列与上面的基序匹配，每匹配一个列中的一个条目的位置得分+1。

LTY.... 

前三个氨基酸得分2分。(L匹配第一列中的一个条目，因此为得分贡献1。T匹配第二列中的一个条目，因此为得分贡献1。Y和第三列不匹配，所以对分数没有影响。)

要将候选CDR3声明为CDR3序列，它必须得到至少10分。此外，左翼必须至少贡献3名，右翼必须至少贡献4名。

接下来，我们找到了在叠架上V段末端的隐含停止位置。这是V段在叠架上的起始位置，加上V段的长度。然后，我们要求CDR3序列在停止之前最多启动10个碱基，在V.停止之后最多启动20个碱基(这一段的条件不适用于denovo的情况)。

如果有多个CDR3序列，我们选择得分最高的那个。如果有平局，我们将选择在叠上较晚开始的那一个。如果仍然有一个平局，我们选择较长的CDR3。

如果通过精确匹配共享相同的CDR3核苷酸序列，则将细胞条形码分组为克隆型。请注意，对于B细胞，CDR3内的体细胞突变将破坏实际上与之相关的克隆型。在CDR3外发生体细胞突变的细胞将被认为具有克隆型。

对于每个克隆型和每个CDR3，所有细胞中的contigs组装在一起，产生克隆型一致序列。

因为这个序列是用多个细胞构建的，所以它的准确性比用单个细胞构建的序列还要高。

名词解释

• CDR3 (Complementarity-Determining Region 3)

三个决定互补的区域是T或B细胞受体氨基酸序列的部分，它们被预测与抗原结合。编码CDR3的核苷酸区域跨越了V(D)J连接，使得它比其他cdr更加多样化。这是一种识别独特链的有用方法。

• Cell Barcode (10x Barcode)

这是一个已知的核苷酸序列，它是单个宝石液滴的唯一标识符。每个条形码通常包含对单个cell的reads。

• Clonotype

通过精确的核苷酸匹配，收集共享一组CDR3生产序列的细胞。

• Consensus

对于一个单一的克隆型和链，在所有的细胞与该克隆型之间建立的共识为该链。这个共识是通过从克隆型细胞中重组相应的contigs而建立的。

• Contig

Contiguous sequence of bases produced by assembly.

• Full-length

A contig is full-length if it matches the initial part of a V gene, continues on, and ultimately matches the terminal part of a J gene.

• Productive

See here.

• Gem Group

当将不同组的gem库合并到一个分析中时，我们在每个读取的条形码上附加一个小整数in silico，以识别读取的来自哪个库。这可以防止条形码冲突，否则会在虚拟双重态的形式中造成混乱。

• N50

The N50 of a sorted list of numbers is the midway point by weight. Example:

image

There are implementation differences for exactly how this is computed but they matter little when the list is long. Unlike the mean and median, the N50 discounts the contribution of many small numbers. That is why people use it!

• N-statistic

n -统计量，如N50或N99，是基因组学中常用的中心性度量，因为它们对大量低价值元素的污染具有一定的健壮性。特别是NXX的值是最小的元素子集包括最少,最大的成员,这样的值的总和子集至少是XX %的总额的值数据集。N-statistic表明一个更大的一个大的值可以占总数的比例大的个人价值,对于一个给定的数据集和YY大于XX, NYY wi的值。

• UMI (Unique Molecular Identifier)

Each first-strand cDNA synthesis from a transcript molecule incorporates a random 10 bp nucleotide sequence next to the cell barcode called the UMI. The UMI sequence in each read allows the pipeline to determine which reads came from the same transcript molecule. In other words, the cell barcode distinguishes between cells, and the UMI distinguishes between molecules (for example, RNA fragments) within a cell.

References

[1] 浅谈BCR及TCR基因的克隆重排: https://www.jianshu.com/p/1d95d9a2561f [2] 10× Genomics单细胞测序在免疫组库研究中的应用: https://links.jianshu.com/go?to=https%3A%2F%2Fzhuanlan.zhihu.com%2Fp%2F44788585 [3] BCR: https://www.jianshu.com/p/5df85f7a3524 [4] Difference Between B Cell Receptor and T Cell Receptor: https://links.jianshu.com/go?to=https%3A%2F%2Fwww.differencebetween.com%2Fdifference-between-b-cell-receptor-and-vs-t-cell-receptor%2F [5] 免疫细胞基因重组与免疫组测序: https://links.jianshu.com/go?to=https%3A%2F%2Fabego.cn%2F2019%2F06%2F25%2Fthe-knowlege-of-immunogenic-recombination-and-ir-seq%2F [6] single-cell-vdj: https://links.jianshu.com/go?to=https%3A%2F%2Fsupport.10xgenomics.com%2Fsingle-cell-vdj [7] 免疫组库高通量分析工具：IGoR——更精确剖析免疫组库: https://links.jianshu.com/go?to=https%3A%2F%2Fblog.csdn.net%2Fwoodcorpse%2Farticle%2Fdetails%2F86762325