免疫组库高通量分析工具：IGoR

Marcou Q , Mora T , Walczak A M . High-throughput immune repertoire analysis with IGoR[J]. Nature Communications, 2018, 9(1):561.

撰文：王庆开 萨萨里大学

责编：刘永鑫 中科院遗传发育所

题记

个人认为, 当前对于免疫组库的研究还处在上升期，如下图是一份文章发表统计，热度似乎没有微生物组高，但是这并不意味着免疫组库研究的落后，而是存在难度，个体的免疫组库多样性非常之高，举个简单例子，我们在用Hi-seq 2500或者Mi-seq对微生物16S测序时,一个lane可以放置200-500个样本，但是如果用PE-150测免疫组库（人的），如果要得到好的数据，目前来看放置20个就比较多了（也根所用试剂盒有关），有文章放了大大小小100个样本，但我读完后没感觉。我在诺禾致源测序时，他们甚至不接受免疫组库测序，回复是：难度大！raw data里一堆无用数据，影响产出。如果有科学家破解了MHC分子和TCR之间的关系，无异于拿到圣火令，练成乾坤大挪移，登上教主之位！而且作者文末也提到了这点，可惜尚未有人做到。

Citation: “Scott C , Walter S , Eddie S , et al. VDJServer: A Cloud-Based Analysis Portal and Data Commons for Immune Repertoire Sequences and Rearrangements[J]. Frontiers in Immunology, 2018, 9:976-.”

摘要

在医学和生物学领域，免疫组库的高通量测序可以看做很有前途的新统计诊断工具。如果要成功的运用这些工具，就需要正确的描述，分析和解释这些数据。我们发布一款全面的分析软件：IGoR (inference and Generation Of Repertoires), 该软件能接受B细胞和T细胞的受体序列（应当是BCR/TCR CDR3）,并能定量地对以cDNA或gDNA为模板的所获得的受体蓄力作出统计分析。它能从概率角度对序列作出注释，而且它的模块化结构针对不同生物可以用来研究增长的生物复杂性模型。对于B细胞，IGoR能返回超突变（SHM）统计结果，沿着序列，揭示超突变位点的共定位。对比当前的类似软件，IGoR具有更好的准确性，能对获得可靠免疫组库描述所需要的样本容量作出估计。（这里可能有个语法错误，estimate应该是estimates, 主语是IGoR,而并非We。）

背景

适应性免疫系统可以通过表达在B细胞和T细胞表面的特异性受体结合对应的抗原。高通量免疫组库测序的出现让我们更好的了解BCR和TCR的多样性，更有可能在诊断、治疗和预防免疫系统紊乱相关疾病方面改变我们的方式和方法。当前，有诸多算法和软件得以解决这个问题，尤其是对序列的分析、胚系基因分配（#我自己的命名，理解为根据测序序列找到原始的基因组序列，类似于基因组测序中，通过测序序列去mapping完整的基因组，作出定位，反过来通过这种定位，我们知道这个片段最初来自于哪个片段！）、以及克隆重构。无论如何，每个受体序列的产生可以依赖很多方式，称之为“方案”，主要涉及前受体的选择：胚系基因片段的重组、插入和删除以及超突变（其实叫体细胞超突变，因为这个时候跟受精卵半毛钱关系没有，完全是单个B或T细胞自己的事情）。不同的胚系基因片段通过不同的频率各自相互组合，而且插入和删除都是随机的，以至于我们无法决断的对整个受体产生过程作出描述。当我们试图描述这种固有的随机过程时，标准的分配会引入系统性错误（我没能理解这句话）。通过定量描述这些机制的多样性和偏好性对于理解适应性免疫和免疫组库测序的诊断应用都是一个挑战。

IGoR 能处理任意来源的原始序列，不管是以cDNA为模板还是以gDNA为模板，它还能学习V(D)J重排和体细胞超突变的无偏统计。通过这些统计，IGoR能对每一个序列输出一份完整的潜力在组合和超突变方案，以及对应的可能性。IGoR识别正确的VDJ重排方案的能力是当前最好的类似软件/方案的三倍。它还可以作为序列发生器，在指定相同统计情况下，产生任意数量的重排序列（应该是模拟学习功能）。该软件应用到BCR上，它能学习一种依赖环境的超突变模型，以此发现突变热点，这种功能兼顾了BCR突变情景的整体分析。

结果

重组方案的概率分配

V(D)J 重排会从胚系基因库中选择2-3个胚系片段（Variable-V 和 Joining-J 片段形成TCR的链和BCR的轻链；V，Diversity-D和 J片段形成TCR链和BCR的重链），在连接处，通过删除碱基或者插入无模板的碱基将这些胚系片段组装起来（fig.1a）（其实1a并没有很好地说明重排过程，仅是大概的过程）。

B细胞在亲和力成熟的过程中可以通过体细胞超突变继续增加多样性。重组过程会变得简化，因为同样的序列可以由多种方式产生（#这是免疫组库的特点，不同的基因片段最后翻译出来的蛋白链可能是一样的，DNA可能一样，也可能不一样）。IGoR可以读入原始的预处理测序数据，监测读长质量，对唯一序列或相似序列分组（这点其他的软件也可以完成同样的处理，比如MiXCR）。IGoR对数据集中已观测到的序列开始列出可能的重组和突变方案。接着，它拟定出可行方案的概率权重。正如Fig.1a所示，探寻/发掘所获得的方案各不相同，但可能具有形成序列的可行性。因为探求所有的可行性会造成计算机资源的浪费，所以IGoR会对这种功能有所限制（补充材料注释5）。根据受体链的不同，探寻方案在一个单核计算机上所耗费的时间从1ms到1s不等（表1，运行时间的所有分布见补充材料1）。通过指定始发事件间的依赖性（比如基因选择，删除，插入和超突变），IGoR可以形成不同的重组方案，而这些始发事件。可通过贝叶斯网络展示，如Fig.1b所示BCR重链的形成方式（其他形成架构见方法部分）。注意一点，IGoR实际上在长读长（longer reads）上运行的更快, 它倾向于消除V基因拟定/分配时的模糊性，从而检少方案探索的总时长（#长读长信息更多更准确不用浪费时间去匹配，但是相对而言最长能有多长呢？一般PE150足够了，CDR3不长）。

IGoR有三种工作模块：学习，分析和生成（fig1c）。 在学习模块中，IGoR通过使用稀有期望最大化算法（该方法常用在机器学习和聚类中）推断大数据集中的重组统计（见方法部分）。在分析模块中，IGoR根据概率拟定序列的重组方案，进而输出最有可能的重组方案排行和整体上序列产生的可能性。在生产模块中，IGoR输出通过指定的统计产生序列例如通过学习真实数据集获得的训练。

V(D)J重组的推断

我们利用IGoR的学习模块对4个数据集的V(D)J重组进行了推测，四个数据集包含了三条链的唯一的非可编码序列（#在IMDB中将序列划分成了几类，一般分析过程中都是用的可以编码产生多肽链的functional sequences，而在任何的一个过程都可能造成移码，导致序列没有编码功能），要么以mRNA为模板扩增获得（TCRα chain or “TRA,” and TCRβ chain or “TRB”) 要么以 DNA为模板 (TRB, and BCR immunoglobulin heavy chain or or “IGH” from naive B cells），也因此形成了早起的方法。将序列限制成不可编码的唯一序列可以允许我们避免由功能性选择引入的偏好。对于BCR，我们只考虑原始免疫组库中的移码/非可编码序列，这些序列没有超突变，插入和碱基删除，增加了它们是非可编码序列的可能。期望最大化算法在经历了几次迭代后开始收敛（见补充材料Fig.2）。

摒弃所有，不管序列来自哪个个体，数据来源，测序方法，或者模板来源D（DNA，浅蓝色分布，Fig.2；mRNA深蓝色分布），只看我们所采用的用最终序列，通过推断得到TRB的碱基插入和删除分布是相同的。相反，V和J基因的使用却有一定的差异，但在个体间甚至测序方法见表现出重大差异，表明存在引物偏好（16S中也存在同样的问题，不能全面覆盖所有细菌，或者引物扩增偏好）（补充材料Fig.3的TRB以及补充材料Fig.4的IGH D-J usage）.正如之前报道，在 TRA V-J,TRB V-D和 D-J见的插入有相似分布（Fig.2a）。和之前的观察一直，IGH序列的碱基插入比TCR的更显著。碱基删除的统计（补充材料Fig.5），特别是negative删除（按照惯例，对应于回文序列依赖的插入）（negative deletion我不理解这个词组，我在国外教材里也没见过），依赖于基因片段，在BCR和TCR中都是如此（补充材料Fig.6），重新认识了此前对IGH的估计。

（辅助理解N-，P- addition）

我们还在“人为的”的数据集上验证了学习算法。 采用对60bp DNA 为模板的TRB数据得到的统计，IGoR产生了10E3-10E5不等的序列。“人为的”序列在长度上匹配实验中所分析的序列长度。IGoR的学习算法在这些原始序列上运行后，所获得的统计和已知的真实数据相比较。我们发现IGoR对一个10E5条序列的推断准确性非常高，在验证过程中我们把容错率设定为10E-3(Fig.3b)。无论如何，不是所有的高通量测序数据都能达到这样的测序深度（#我的测序就都达到了，但是样本是血液，如果是组织肯定无法达到，只有血液能做到，因为任何研究者都无法把整个组织的mRNA都当做模板），特别是把序列限制到非可编码水平（#虽然作者是用非可编码序列进行验证，验证结果也能看到VDJ的组合情况，但是体内的真实情况是不会选择这种序列，现实中是没有意义的，不编码意味着这些克隆会凋亡，会被淘汰，而在B细胞和T细胞进化成为成熟的过程中，被淘汰的这两种细胞占了非常大的比例，无用的细胞种类比最终成熟的要多得多，所以作者上面所说的很难达到10E5可能指的是这个。）。为了评价这些限制因素对准确性的影响，我们计算了在不同大小数据集和容错率下真正的TRB分布和推断的分布之间的Kullback-Leibler散度（一种针对不同概率分布的非参数方法，见方法部分）。在容错率10E-3水平下，~5000个唯一的移码序列（可从2mL血液中获得）经过充分学习建立一个准确的TRB模型（Fig.3c,同见补充材料7的插入删除分布），多数的估计偏差是由于删除模式造成（删除模式用来解释多数参数）。T典型的KL散度是很小的，表明较少的过度拟合。Increasing the error rate has little effect up to rates of 10−2, but significantly degrades accuracy for larger error rates, 10−1 (Fig. 3d), with the gene usage distribution affected the most (Supplementary Fig. 8). 把容错率增加到10E-2不会有太大影响，但是允许更大的容错率（10E-1）却能显著降低准确性，基因使用分布受影响最大（补充材料Fig.8）。此外，IGoR会把BCR的超突变同样当做出错处理，达到1-10%，而且表现出长的拖尾，约占5%的序列，这些序列都有30%甚至更高的突变率（补充材料Fig.9）。这表明BCR的重组统计可以通过使用初始的、非超突变的B细胞得到更好的推论（正如Fig.2所示）。我们也使用“人为的”数据集的移码序列验证了IGoR的学习模块，显示不会对推断结果造成影响：模型是通过训练“人为的”TRB移码序列获得，相比于真实模型仅有0.4bits的差异，而同时在可编码和非编码序列上训练时，也只有0.3bits。

方案的退化/退步分析

考虑每条序列的所有重组方案，我们的方法和现存的多数方法都显著不同，那些方法的目标是发现最可能的重组方案（但是，Partis基于概率形成方案。）。为了评判哪个序列产生方案是正确的，我们分析了“人为的”序列，这些序列的构建方案在一个模型中是已知的，并且没有超突变。 For each generated sequence, we used IGoR’s analysis mode to enumerate the set of scenarios that were consistent with the nucleotide sequence, and ranked them according to their likelihood. Figure 4a shows the distribution of the rank of the true recombination scenario for TRB and unmutated IGH synthetic data.对于每条产生的序列，我们使用IGoR的分析模块列举出方案组合，这些方案组合与核酸序列对应一致，根据可能性大小排列。Figure 4a 显示的是“人为的”TRB和未突变的IGH数据的重组方案排列分布。

最可能的方案却不是正确的那个，这在IGH序列中占比72%，在TRB序列中占比85%。 两种分布都存在长的拖尾，意味着有相当一部分序列的重组方案存在重组退化。这种退化并不是IGoR不能学习正确的模型，而是由于VDJ重排过程中内在的随机性所导致：插入和删除核酸的不同组合能精确地产生同样的序列。这种内在的限制因素似乎不能被改变，不管有多少数据量，也不管更长的读长，或者不同的重构模型。
接着我们评估，要解释所有序列可能性的一个既定分数f需要多少方案，方案从最可能到最小可能排列。IGH的100000条序列的f分数的分布如Fig.4b所示（补充材料Fig.10是TRB的f分数分布）。如果要列出正确方案中的f在95%置信水平，一般需要至少30-50中方案。这表明，对于TCR和BCR而言，许多方案需要正确地考虑描述出产生过程。
IGoR通过对所有方案的可能性求和，输出要处理的序列产生的概率，那些确定性的分配方法无法做到这一点。这种序列产生的可能可以预测健康个体间的共有特性。在一些试图发现抗原特异性或者自身免疫相关序列的临床研究数据中，这种功能可以被当做重组方案收敛的指示标识。

与其他方法比较

我们将我们的方法和当前两个具有代表性的算法相比较：一个是MiXCR，该软件能很好地对胚系基因作出匹配；另一个是Partis,该软件只针对BCR，通过搜寻超突变利用最大似然法找到最可能的方案。借助IGoR的生产模块，产生出“人为的”长度为130bp的无超突变的IGH序列数据集。将这个数据集应用到MiXCR，Partis和IGoR三个软件，比较他们在序列产生中真实方案的输出。在IGoR和Partis中，模型参数是通过学习已产生数据集来模仿真对真实数据的分析。Figure 4c显示了三种方法分配正确的重组方案的表现。在预测完整的重组方案方面（针对无超突变的IGH），IGoR的表现要好过其他两款软件两倍有余，每一个模块皆是如此。值得注意的是Partis的处理并不包含回文序列为模板的插入，其余两款软件通过在每个胚系基因的末端引入一小段回文序列进行处理；这点限制导致Partis的表现下降到和MiXCR一样（补充材料Fig.11）（#反观，如果引入这中分析功能，Partis的表象将好于MiXCR）。当指定V基因时，由于更长的读长增加了序列信息，所有的三种方法都表现的更好，但是这也没有改变其余两款软件在发现插入和删除图谱的能力，或者正确的召集真正的序列产生方案的能力。因为IGoR能列举所有的可能的方案，运行速度要比其他方法慢一些（见table 1）。

为了核实IGoR在序列组装表现上的改变，尤其是重组后的选择效果，但IGoR针对重组后选择并不创建模型，我们人为选择性地产生了初始IGH数据集，这一数据集根据CDR3长度选自“人为的”序列，以至于在拒绝抽样后，最终的长度分布精确的匹配初始可编码IGH序列。把三种方法应用到这些序列上和那些未经选择的序列上，三种方法的表现很相似，但是IGoR的模型参数来自未经选择的序列（反应我们的学习模型方法的来自未经选择的序列）还是来自经过选择的序列数据（补充材料Fig.12）。然而IGoR的注释功能很好，推断的重组统计可能在经过选择和未经选择的序列上表现差异很大。
接着，我们比较了三种方法通过真实统计学习得到重组统计。对MiXCR和Partis，我们构建了对应到每条序列的重组方案分布，而IGoR仅使用期望最大化算法构建分布。三种方法都产生了相似的统计：V和J基因的使用、碱基删除图谱（补充材料Fig.13）。无论如何，D和J基因之间的相互依赖性被IGoR所捕获，但MiXCR没有（而Partis没有这种功能）。TRB D和J基因在基因组上是分别两个基因簇中（#如下图所示），D1对J1，D2对J2。由于这样的组织结构，D2不能和J1家族所结合进行重排。而IGoR对这样的组合直接指定为0概率事件。同样的结果也可以直接从真实数据中获得（补充材料Fig.14）。最后，IGoR准确的对插入分布重构，而其他两款软件系统性的高估了0插入事件的可能性（补充材料Fig.13a,b for IGH）。

IGoR也具有一些repgenHMM17的特性，repgenHMM17这款软件利用隐马尔科夫模型从概率角度给出重组方案。无论如何，比起IGoR， repgenHMM17运行较慢（表 1），而且不能召集最有可能的重组方案，也不涵盖超突变模型，还不能形成重组方案间的依赖性关系。为了阐述最后这点，我们通过一个合成的TRA数据集学习了一种模型，该数据集是模拟人为设定的V基因选择和碱基插入。但是这种依赖性在TRA中是不存在的，只是和其他受体相关（比如TCRδ and TCRκ），抑或针对指定的受体。（#这里可能存在知识性错误，在人体内TCR不存在κ链，只有，和，，但是IgG倒是存在一条κ轻链，审稿不严格，其实从一开始就感觉到这不是一个专业的免疫出身的人写的文章，或者没有经过免疫研究人员的review，在多数的遗传学书上，重组叫recombination，说的是染色体之间的事件，或者大段的DNA。Rearrangement才是说的V,D,J 之间的组合情况）。简言之，任意一种插入分布（真实的，或者呈几何分布的）的V基因都是随机分配的。此外，在学习正确的V基因插入分布方面，这点repgenHMM无法通过构建实现，IGoR学习到了删除分布，这点基因依赖结构的特性和repgenHMM采用同样的算法，但准确性比repgenHMM高很多（补充材料表格 1）。

体细胞超突变

为了研究记忆性B细胞上BCR的体细胞超突变发生模式，我们给IGoR加入概率机能，用来推断序列依赖的超突变率。在核酸序列上某个既定位置，出错或者突变的概率假定是依赖序列于相邻的nmer环境（Fig.5a），通过额外赏分转换得来，这一赏分通过位置权重矩阵计算而，类似于用来描述DNA结合位点的结合能基序。用IGoR运行文献14 中的移码IGH序列去学习7-mer的PWMs,也学习整体的突变率（在所有7-mers上的突变率的几何平均数），但要固定住此前通过期望最大化算法从初始序列学习得到的重组统计量（方法部分），IGoR的概率架构可以应对基因片段选择或者超突变的收敛组合引起的序列起始的退化。这种学习过程显著不同于文献16的方法，在文献16的方法中，超突变率是统一的。这点也和Partis不同，Partis不通过学习PWM模型，却硬生生地推断出超突变图谱,直接附有每个基因片段上位点的功能（类似于IGoR的posterior mutabilities, Fig.5c）。鉴于此，Partis通过利用更多参数，学习了一个更完整的超突变图谱。其他环境依赖的模型也被提出，这些模型没有附加的假设条件，因此大量的参数都是从其他无义突变统计或者突变小鼠的非功能序列学习获得。三个不同的PWMs是为V,D,J的模板区域而来（Fig.5b）。为了证实我们的PWM和突变率学习算法，基于之前的真实数据及学校得到的模型，我们生产了带有超突变的人工数据，让IGoR再次学习这些参数，发现一致性很好（补充材料Fig.15）。

超突变的位点依赖概率的PWM预测和在序列中实际观测到关联的很好（r=0.7 V基因中，Fig.5c，以及补充材料Fig.16）。在两个测试的个体中，PWM完全可重复（r=0.98,补充材料Fig.17），表明对个体的感染史的推断程序是可行的，也指明了体细胞超突变的普遍性。相反，两个体间推断得到的整体突变率存在两倍的差异，可能由于个体的年龄，既往感染史或者生活方式不同所造成（补充材料Fig.17）。对于每一个基因，我们所发现的基序囊括了之前所报道的热点基序（PWM的正值所示），包括WRCY（或者WRCH）以及WA（W=A或T，Y=C或T，R=G或A；突变位点被标出），也发现了非热点基序，及时较少（SYC，S=C，G）。在上述三种基序中，除了在V和D基因中的突变位点外，C和G位点普遍不出现，T也很少发生突变。通过学习3-9 mers的PWMs，我们对模型的强度进行了评估（补充材料Fig.18）。作为一种context length 功能，每个相对位置的贡献并没有实质性改变。距离突变位点至少达到4个碱基才能对继续有所贡献。这可能意味着前后依赖是广泛的，或者基序模型间接捕获了“非前程后继效应”。整体上，推断的PWMs给出关于规则更详尽和更细致的描述，这些规掌控了热点位置，且不能被简化为几个简单的可描述基序。

Figure 5b 展示了V，D，J基因上基序实质性差异。基于V基因学习的PWMs仅能适度地预测J基因的超突变率（r=0.5 对比r=0.7，V基因上的突变率），基于J基因学习的PWMs预测V基因上突变率的效果更差（r=0.24, 补充材料Fig.16）。为了评估是否这种差的普遍性归结于模型中附加的假设，我们学习了一种“非附加5-mer”模型，这种模型试图分配一个特定的超突变率给1024种可能的5-mer的任意一个（补充材料注释2）。实际上，我们仅能估计一少部分5-mer的突变能力，这一少部分也只是占所有可能的1024种5-mer的不到一半（在160-498个之间）。虽然基于V和J片段学习的模型在个体间是可以通用的（补充材料Fig.19a,c），但并不是很好地相互匹配（补充材料Fig.19d,e）。尽管作为位置功能，模型可以很好地预测突变能力，也只针对少部分5-mer才能行得通（补充材料Fig.20a,b），但是基于V片段学习的模型却不能应用在J片段的超突变率预测上，这与附加模型结果一致。这种不一致表明纯粹基于前后依赖继续的预测不能够解释所有的超突变可能性，而且其他机制存在这样问题。在胚系基因之间，整体的突变率是不相同的，与所报道的染色体状态影响超突变率是一致的。

接下来我们用推断的PWM在IGoR中从概率角度招募序列中假定的超突变。首先，我们检查一条序列中突变数量的分布（Fig.5d）。实际分布（红色）变得更扭曲，相比于通过假定有独立的超突变的预期有更长的拖尾， 正如通过推断的PWM（蓝色）所预测随机产生的超突变序列一样。这一观测和实际的B细胞在亲和力成熟过程中所经历的多次循环相符合。 其次，通过计算两个位置间的超突变富集情况，我们不禁要问，是否超突变在相同的序列上共同存在（Fig.5e）。然而这种富集情况在合成序列中只有一个（因为我们的模型假定超突变之间是相互独立的），真是的数据显示超突变在临近位置有高达4倍的富集。这一不同点与一种实际情况一致，就是AID能引起大片段DNA的修复。这种经典的距离可以对超突变相关区域的长度作出预估，大概15bp，在这一距离上共定位富集指数衰退。

原则上，IGoR能计算任意序列产生的概率。但是，面对高突变序列却存在困难，因为胚系重组序列本身有时候也是不确定的。为了克服这点，IGoR对每条序列都探索出可能的重组和超突变方案，并且计算出每条潜在胚系序列的产生概率。利用合成的数据，我们监测到个体序列的产生概率可以被这种方法很好的预测（r=0.97,补充材料Fig.21和方法部分），而且概率分布是可重复的（补充材料Fig.22）。

讨论

从概率角度，通过对免疫受体比对到胚系基因上，IGoR可以纠正在V(D)J重排统计过程中的系统性偏好，还可以比之前报道的方法更准确的预测重组方案。它详细的重排方案分析能力进一步揭示了由于重排过程固有的随机性，有超过70%的序列重排方案是不正确的，表明在解释确定好的分配方案是要谨慎。

尽管我们只在人类的TRA，TRB和IGH上进行了验证，但IGoR的灵活架构让它可以在任何可变的淋巴细胞受体（TCR或Ig）和物种（前提是这些物种有现存的胚系数据库可用）上运行。和基于隐马尔科夫模型的方法不同（文献10和17），IGoR涵盖更多可能的重组事件中的依赖关系。它还适合处理特殊的或者不完整的重排（D-J重排，TCR δ链的DD2/DD3重排，杂合TRA/TRD重组等）。IGoR还能借助自身合成的对照序列检测特殊的重排。例如，串联D片段的重排已经被报道，但是要想辨识出随机插入还是有些困难。为了测试这点，我们用两个长度大于等于10 nt的D片段作为标记提取这些序列，并用其中一个D片段以此和IGoR预测合成的序列相比较（方法部分）。我们发现IGH数据中的双D片段分配是对照中的五倍，证实了文献14中的发现。相反，在TRB中验证却没有发现串联D片段的存在。未来要发布的IGoR版本应该可以支持多D片段检索的功能。但我们发现IGoR不能在IGH中发现反向的D片段（补充材料23）。

IGoR可以在一个5000个非编码序列的数据中推断重组统计量。一旦对一个既定基因座学习的到了重组模型，IGoR可以以相同的统计量产生一份专断的合成序列，这些序列可以在疾病关系研究中被用作对照，从统计角度以较高的收敛重组频率对公开的序列中识别抗原特异的克隆性，从而免除招募健康个体对照的必要。这种方法基于个体将受体产生过程高度的可重复性，一个个体产生的模型可以应用到不同个体上。基因的分配图谱可以说是最具个性化的（特别是是V基因的选择，被认为和HLA类型相关），但是，他们对受体序列的产生整体贡献并不大。 为了控制不同基因分配时带来的偏好，可以对V-J类型进行族群的分析。要么，在使用IGoR分析时，固定插入和删除模式，而对感兴趣的V和J基因的使用进行推断，因为这些似乎不依赖来于个体或者方法。

我们对超突变的分析使得我们能推断人IGH V,D,J片段上不同基序上的突变靶点，和之前的其他方法有所不同，那些方法假定一个普遍的“上下联动”模型（#自主命名，因为作者表达的意思，突变也依赖周围碱基的状况）。在我们的工作中，我们主要集中在附加模型，因为非附加上下联动模型会受到n-mers数量的限制，突变能力也因此被可靠的估计。我们核查到我们的结果并不是简单地附加假设的结果，当学习一个非附加5-mer模型是他们就会受到约束（补充材料Fig.19）。尽管非附加5-mer模型比附加模型效果要好，但也只是在5-mer水平上超突变率能被可靠估计，这也解释了在所有最好的案例中所有的5-mers不到一半的水准（补充材料Fig.20）。此外，任意既定的5-mer环境是很稀少的，基于此它可以发现独特的基因片段和基因上的位置，漂亮的结果让它看上去是过度拟合的结果。

我们把我们的IGH突变率和“S5F”模型的预测结果做了对。和我们对附加5-mer模型的结果一样，“S5F”模型的预测也不好（Pearson’s r(0.2-0.37),补充材料Fig.24）。但这点在非附加5-mer模型上却有所改变（Peason’s r（0.45-0.47），尽管这种改善需是基于抽取部分较好5-mer。S5F模型的预测能力和附加模型的预测能力相当（补充材料Fig.26）。尤其是，该模型对J片段突变能力的预测较好，尽管只是在V片段上进行了训练。相比于IGoR，值得注意的是S5F模型是在更长的读长和更多的数据及上训练。尽管我们对合成序列的分析标明继续主要在短序列上准确学校得到（补充材料Fig.15），但如果把S5F模型用到的长序列给到IGoR上训练，也会产生很夯实的模型，预测能力更好。

我们还发现超突变倾向于在序列上集中出现。总的来说，至少三个效应决定了超突变热点：超突变邻近的DNA，正如我们的序列基序模型给出的，染色质构型，组蛋白修饰等介导的位置特异效应，还有就是邻近突变的出现。未来在超突变靶点预测的改进可以涵盖以上三方面，而且依赖于对AID介入的更好的量化理解。

除点突变以外，超突变过程也包含indel. 我们估计大概5-12%的记忆性B细胞的IGH序列是非编码的（由于移码框或者终止密码子在CDR3中出现），这些序列在V区有indel（补充材料Fig.27）。IGoR会丢掉大概30%这样的序列，这是由于较大的罚分会把这些序列降低到最小可能阈值以下。尽管当前IGoR的架构还不能对indel建模，但我们检测到这并不影响我们的结果，那就是：增加可能的阈值以此去掉所有包含indel但不影响推断输出的序列。如果IGoR能对indel建模，这将更好的帮助分析那些由于超突变率增加引起的包含高indel率的免疫组库，比如说那些HIV的携带者。

IGoR描绘了VDJ重组过程中基本元件，并给出产生一个既定的TCR或者BCR序列的整体可能性。为了描述基因使用，插入，删除图谱等这些过程的统计量，模型需要在那些尚未经历任何选择的序列上进行学习。在此，我们将移码序列当做这种序列，但是还能正确的筛选出双阴性TCR或者如果能利用，原B细胞也是不错的选择。在那些未经筛选的序列上训练过，IGoR的分析模块也能对筛选过的序列进行注释，还能计算它们产生的可能性。因为选择性只影响序列本身，不对重组方案有直接影响，所以IGoR应该可以在经过选择的和未经选择的序列上运行（补充材料Fig.12）。

此外，IGoR计算得到的序列产生概率可以被用来消除选择带来的效应，还能发现功能序列。IGoR的模型也可以直接在既定免疫组库的可编码序列上训练得到，即使这个免疫组库尚未经历过功能选择（补充材料Fig.12）。然而模型架构可能不适用于描述作用在蛋白序列上选择压力，这样的模型需要能同时捕获序列产生和选择压力，而且可以用来估计经选择过的免疫组库中特殊序列的泛性。一般而言，集中所有工具，重要的一点是要理解数据的局限性，获得和预处理的局限性，这都影响最终的结果解释。IGoR可以对任何平台上获得数据进行处理，甚至来源于不同的是实验，而且它的表现并不依赖于实验技术，值依赖预实验数据的可信度和正确的处理方式。

此外，当前的IGoR版本还需要在运行时输入胚系基因数据。但是这样数据库通常是不完整的，因为他们并不包含一个族群可能的全部多态性。原则上，频繁地被IGoR或者其他软件发现的“错误”可以被用来检测多态性和推断大型不完整数据库的信息丢失，比如一些了解不彻底的物种的免疫组库测序。这种方法需要充实的基因组知识，在没有完整基因组的情况下以此来构建引物和实验设计。这也是IGoR的新的应用。
总结，因为特定的序列更有可能被产生，所以IGoR可以通过设定基线，让我们看到感兴趣的序列。它可以被当做工具从功能上经选择过的序列中识别收敛的重组，也可以和其他更准确的突变模型组合起来联用。

方法

IGoR方法概述

IGoR有三种工作模块：VDJ统计量学习，序列分析，和序列产生。爱国有的模块依赖重组和超突变事件的明确的、随机的描述。在分析和学习模块，每条序列都会被给出所有可能的重组和超突变方案。学习模块会反复的了解分析模块，并持续更新自身模块由期望最大化算法获得的参数。

重组建模。在三个模块中，IGoR假定受体序列经一个重组方案得来，这个方案包括几个随机元件：胚系基因片段的选择，插入，删除。这些特性是随机的，相互之间共有相同的统计依赖关系。为了方便，我们假定这些依赖是可以通过非循环图表示的，也称之为贝叶斯网络（补充材料注释1）。这种结构可以在IGoR的设置文件中进行配置。该研究中，我们使用如下依赖结构针对TRA：

对于TRB和IGH：

V,D,J表示胚系基因的选择，delV,delJ表示在V，J片段中删除的碱基对数量，delDl,delDr表示D片段两侧删除的碱基对数，insVJ，insVJ，insDJ表示插入的碱基对数量（V-J，V-D，D-J之间），ni,mi表示插入的碱基类型。在TRB中，基因使用进一步分解为P（V,D,J）=P（V）P（D,J）。

超突变联动建模：（#自己命名，因为context难翻译，兼有上下文，环境的意思）。当处理TCR或者初始BCR时，假定整个序列存在一个很定的容错率。当处理记忆性B细胞BCR时，联动突变模型假定在V,D,J基因的每个位置上，每个超突变的发生都以Pmut的概率发生，使用公式如下

在公示中（π−m, …, πm）是（2m+1）-mer序列的表示，围绕突变位置。在PWM中，ei(π)是对基序的额外贡献，μ是整体的超突变率。

胚系基因比对和方案枚举

在学习和分析模块中，每条序列首先要和所有可能的胚系基因通过Smith–Waterman算法比对（比如IMGT数据库的胚系基因列表）（#算是当前的主流胚系基因数据库，多数软件可能都以此为参照），只有那些评高于阈值的胚系基因才会被用来进一步分析（补充材料注释5）。可能的重组方案可以通过选取打分高的胚系基因列出，并选取较多的碱基对从比对的末端删除掉。在砍掉多余胚系基因后，以此种方式在胚系基因片段之间的碱基对被称作插入，和胚系基因比对得出的错配成为出错或者超突变。当胚系基因的回文末端被不完全删除时，剩余回文序列碱基对称之为负删除。考虑到在BCR中两个方向上都有D片段插入的可能，我们添加每一个D胚系基因的反向互补序列到所列举出的胚系基因。

序列分析

对数据集中每条序列，每个可能重组方案的概率都会使用重组概率Eq.(1) 或（2）乘以出错概率或者超突变概率Perr: Pscenario = Precomb × P err 计算得来。重组方案按照重组概率递减的方式列出来。所有重组方案和超突变的概率之和既是所观察到的这条序列的概率，Pread。突变前序列的概率通过重组产生，Pgen,定义为所有可能形成这条序列方案的概率（Precomb）之和。因为突变前序列是确定未知的，所以我们计算一个大概的产生概率Pgen作为所有和此读长一致的未突变序列Pgen的几何平均数，而且给予权重，Pgen × Perr/Pread。要么，估计出所有可能未突变的概率Pgen(补充材料注释4)。

为了减少计算时间，IGoR只列出那些近似合理的重组方案。重组方案以分层决策树节点的方式列举出，每个深度对应方案特性的选择。如果它们对序列概率整体的贡献的上限低于特定阈值，这个树枝就会被丢弃。具体步骤见补充材料5.

学习算法

学习算法在一个大的唯一序列的数据中推断Eq. (1) or (2)以及出错或超突变模型参数Eq. (3)。它依赖于序列分析模块，而且遵从期望最大化算法的过程。伪对数似然率定义为所有序列可能方案的对数似然率的权重之和，权重是既定序列方案的条件概率，Precomb/Pread (期望步骤)。当固定权重，这一伪对数似然率较之于对数似然率(Eqs. (1), (2), and (3))就被最大化。参数更新，重复此前的步骤，直到收敛。期望最大化算法数学衍生规则和细节见补充材料注释2。

模型推论的验证。为了比较由合成数据推断得到的θ1和已知数据模型参数θ2，我们计算了两个概率分布间的Kullback-Leibler散度。D(θ1||θ2)=E P（E, θ1）log[P（E, θ1）/ P（E, θ2）],总和是所有方案的E. P(E, θ)，由Eq. (1) or (2)计算而来。KL散度可以被分解为附加贡献，这些贡献源自所有方案的性质，如补充材料注释3所示。

超突变间的关联性

为了证实BCR序列上邻近位置间超突变发生的相关性，我们计算了径向分布函数，定义为

f(i, V) 和 f(i, j, V)分别是在i位置和i,j位置的超突变频率，由个体的重组方案统计量加以权重计算得来。CV(r)是成对位置的集合由r分开，这些成对的位置在V基因中被多次观测到，Nr=V | CV(r)|。

串联D片段的使用。为了评价IGH和TRB的VDJ重组过程中双D片段插入频率，我们通过Smith-Waterman算法计算出这个频率，单次可以至少在10核苷酸上比对到两个非重叠D片段，这些片段都在最好V和J上。接着我们比较了合成序列和真实序列间的频率。

译者简介

王庆开，博士研究生，2009年毕业于山东理工大学，2013年硕士毕业于复旦大学生物医学研究院，主要从事结直肠癌和补体系统的研究（上海市肿瘤医院肿瘤研究所）；后期参加工作从事：焦磷酸测序（肺癌，结直肠癌，乳腺癌靶点检测指导临床用药），小分子药物的高通量筛选（Envision，Echo等高通量工具半自动化）和转基因小鼠销售工作。2017年进入UNISS博士生研究计划，目前从事前列腺癌免疫组库研究和微生物16S的研究以及新发现微生物C.M girerdii 的特异性抗体制备。

猜你喜欢

• 10000+: 菌群分析
  宝宝与猫狗 提DNA发Nature 实验分析谁对结果影响大  Cell微生物专刊 肠道指挥大脑
• 系列教程：微生物组入门 Biostar 微生物组  宏基因组
• 专业技能：生信宝典 学术图表 高分文章 不可或缺的人
• 一文读懂：宏基因组 寄生虫益处 进化树
• 必备技能：提问 搜索  Endnote
• 文献阅读 热心肠 SemanticScholar Geenmedical
• 扩增子分析：图表解读 分析流程 统计绘图
• 16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun
• 在线工具：16S预测培养基 生信绘图
• 科研经验：云笔记  云协作 公众号
• 编程模板: Shell  R Perl
• 生物科普: 肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA