重組炭疽疫苗免疫人體產(chǎn)生抗體基因的生物信息學(xué)分析

劉渝嬌，陳鄭珊，王美榮，張冠英，范鵬飛，房婷，遲象陽(yáng)，于長(zhǎng)明

軍事科學(xué)院 軍事醫(yī)學(xué)研究院 生物工程研究所，北京 100071

抗體是獲得性免疫的重要組成部分，在體液免疫反應(yīng)中起關(guān)鍵作用。其中單克隆抗體藥物因具有高度靶向特異性、親和力和安全性，已經(jīng)在許多重大疾病的治療、診斷和預(yù)防中得到廣泛應(yīng)用，是當(dāng)今新藥研發(fā)的重點(diǎn)?？贵w基因序列的測(cè)定和分析是一個(gè)高速發(fā)展的科學(xué)領(lǐng)域，已經(jīng)應(yīng)用于多種疾病的診斷和防治中，包括疾病微小殘留病變檢測(cè)、移植的預(yù)后檢測(cè)、疫苗系統(tǒng)性評(píng)價(jià)、中和抗體鑒定及B細(xì)胞運(yùn)輸模式追蹤等[1-4]。隨著高通量測(cè)序技術(shù)的發(fā)展和成熟，大規(guī)模、系統(tǒng)性分析抗體基因序列已成為可能。但是，抗體基因的胚系細(xì)胞多樣性和體細(xì)胞多樣性增加了分析的難度，特別是針對(duì)大容量、多樣性高的抗體基因序列庫(kù)，需要借助相應(yīng)的生物信息學(xué)手段以及特定的計(jì)算方法對(duì)抗體基因序列進(jìn)行分析和數(shù)據(jù)挖掘。

炭疽芽胞桿菌是一種致死率極高的人畜共患烈性傳染病原。炭疽疫苗和抗體藥物作為重要的炭疽防治手段，其研究一直是世界范圍內(nèi)關(guān)注的重點(diǎn)。本實(shí)驗(yàn)室研發(fā)的基因工程重組炭疽保護(hù)性抗原（protective antigen，PA）疫苗能夠有效激發(fā)人體產(chǎn)生中和抗體。我們前期通過(guò)流式細(xì)胞分選-單細(xì)胞PCR技術(shù)，從重組PA疫苗免疫志愿者體內(nèi)分選并擴(kuò)增出大量重鏈和輕鏈天然配對(duì)的全人源單克隆抗體基因，包括34株P(guān)A結(jié)合抗體基因和4株中和抗體基因。在本研究中，我們?cè)贚inux操作系統(tǒng)（Ubuntu）上，基于NumPy、SciPy、pandas、Biopython、presto、Change-O[5]等構(gòu)建了抗體基因序列分析平臺(tái)，并利用alakazam、shazam、tigger等R語(yǔ)言包提供的函數(shù)進(jìn)行編程實(shí)現(xiàn)相應(yīng)功能，利用R、Perl等語(yǔ)言編寫了一系列數(shù)據(jù)處理和接口程序，對(duì)前期獲得的炭疽抗體可變區(qū)序列進(jìn)行了一系列高效的生物信息學(xué)分析，包括抗體重鏈的VDJ基因使用和統(tǒng)計(jì)特征分析、基于聚類分析的抗體譜系樹分析、中和抗體和結(jié)合抗體的重鏈CDR3區(qū)物理化學(xué)特性統(tǒng)計(jì)和分析、5-mer motifs體細(xì)胞高頻突變模型構(gòu)建。研究結(jié)果表明，前期獲得的抗體基因序列中IGHV3、IGHD3和 IGHJ4的使用頻率最高，IGHV3和IGHD5、IGHD6和IGHJ4、IGHV3和IGHJ4的組合是序列中出現(xiàn)頻率最高的組合方式；23.3%的PA結(jié)合抗體都出現(xiàn)在同一個(gè)抗體聚類（239號(hào)克隆組）中；中和抗體CDR3區(qū)較其他抗體對(duì)堿性和酸性氨基酸的使用頻率較低，偏向于使用疏水性氨基酸；抗體的體細(xì)胞高頻突變多發(fā)生在AAAAC、GTTAA、AGCTC、AAGTA等熱點(diǎn)motif。另外，運(yùn)用ARGalaxy Web Sever[6-7]的抗體基因序列分析模塊得到相應(yīng)的數(shù)據(jù)可視化結(jié)果。

本研究搭建的抗體基因序列分析平臺(tái)整合安裝了完整的序列處理和分析程序，可離線使用。該平臺(tái)實(shí)現(xiàn)了對(duì)高通量抗體基因序列的快速處理分析，而且具有較好的拓展性，可自由組合使用各個(gè)功能模塊以及編程開發(fā)滿足實(shí)際需求的新功能。本研究結(jié)果對(duì)于重組炭疽疫苗評(píng)價(jià)和生物信息學(xué)輔助的抗體鑒定和改造具有一定的借鑒意義，提供的方法也可為其他相關(guān)疫苗激發(fā)免疫系統(tǒng)應(yīng)答產(chǎn)生抗體譜的研究及中和抗體的篩選和鑒定提供幫助。

1 材料與方法

1.1 抗體基因序列來(lái)源

采集重組PA疫苗免疫志愿者的血液樣品，利用密度梯度離心法分離外周血單個(gè)核細(xì)胞（PBMC），進(jìn)行熒光染色，通過(guò)流式細(xì)胞分選儀從PBMC中分選單個(gè)抗體分泌細(xì)胞克隆（CD19+/ CD3-/CD20-/CD27high/CD38high）。采用反轉(zhuǎn)錄PCR從單個(gè)B細(xì)胞中擴(kuò)增抗體基因，再以該P(yáng)CR產(chǎn)物為模板進(jìn)行單細(xì)胞巢式PCR反應(yīng)，經(jīng)瓊脂糖凝膠電泳鑒定為陽(yáng)性的巢式PCR產(chǎn)物送上海生工生物工程股份有限公司測(cè)序，批量整理測(cè)序結(jié)果即得到抗體可變區(qū)基因序列，保存為單個(gè)fasta格式文件。共計(jì)346條H鏈可變區(qū)序列（其中4條序列為有中和活性的抗體，30條序列為僅有結(jié)合活性的抗體）、268條κ鏈可變區(qū)序列、173條λ鏈可變區(qū)序列。考慮到抗體的特性主要由H鏈可變區(qū)決定，故未對(duì)κ和λ鏈可變區(qū)基因序列進(jìn)行類似分析。

1.2 基因序列的注釋和標(biāo)準(zhǔn)數(shù)據(jù)集的生成

由于使用反向測(cè)序引物，因而需要將測(cè)序所得的抗體基因序列進(jìn)行反向互補(bǔ)后再進(jìn)行后續(xù)分析。利用BioEdit[8]軟件批量操作，得到全部原始基因序列的反向互補(bǔ)序列。利用平臺(tái)中本地化的IgBLAST程序?qū)贵w可變區(qū)基因序列進(jìn)行注釋，IgBLAST程序具有識(shí)別序列的FR/CDR邊界、最佳匹配胚系V（D）J基因、顯示V（D）J基因連接細(xì)節(jié)、提供核苷酸的重排連接信息等功能。進(jìn)而運(yùn)用Change-O程序集提供的MakeDb預(yù)處理程序（Python語(yǔ)言編寫）對(duì)IgBLAST程序的輸出文件進(jìn)行整理和格式轉(zhuǎn)化，生成可供平臺(tái)內(nèi)其他應(yīng)用程序使用的標(biāo)準(zhǔn)tab-delimited數(shù)據(jù)集文件。

1.3 基于抗體基因標(biāo)準(zhǔn)數(shù)據(jù)集的分析方法

對(duì)得到的標(biāo)準(zhǔn)數(shù)據(jù)集文件在平臺(tái)內(nèi)進(jìn)行處理分析，使用R語(yǔ)言包alakazam[10]提供的count?Gene函數(shù)對(duì)抗體重鏈基因數(shù)據(jù)集進(jìn)行VDJ基因使用情況的統(tǒng)計(jì)，可得到VDJ基因節(jié)段數(shù)量、節(jié)段家族數(shù)量等統(tǒng)計(jì)信息。R語(yǔ)言包shazam[5]提供的distToNearest函數(shù)用于計(jì)算調(diào)整抗體基因序列分層聚類的合適閾值[11]，并將該閾值作為Defines?Clone程序（由Change-O程序集提供）的dist參數(shù)值。進(jìn)而運(yùn)行DefineClones以定義克隆組（clone group），對(duì)抗體基因序列進(jìn)行聚類。使用Perl語(yǔ)言編寫的程序?qū)σ丫垲惖男蛄羞M(jìn)行克隆組選擇和數(shù)據(jù)處理，使用PHYLIP[12]軟件包中的dnapars程序?yàn)橥豢寺〗M內(nèi)的抗體序列建立譜系樹，并通過(guò)R語(yǔ)言包igraph調(diào)整和繪制譜系樹網(wǎng)絡(luò)圖。al?akazam包提供的countClones函數(shù)可用于統(tǒng)計(jì)克隆組的情況。alakazam包中還提供了一系列用于分析抗體氨基酸序列理化性質(zhì)的函數(shù)，如利用ami?noAcidProperties函數(shù)分析序列CDR3區(qū)的長(zhǎng)度、親疏水性和氨基酸的酸堿性等特征。shazam包提供的createTargetingModel函數(shù)[13]可通過(guò)序列信息構(gòu)建出體細(xì)胞高頻突變模式模型（SHM targeting models），并通過(guò)plotMutability函數(shù)將模型可視化為直觀的猬形圖（Hedgehog）。在ARGalaxy Web Sever上傳序列文件后選擇ARGalaxy模塊下Im?mune Repertoire pipeline功能，根據(jù)實(shí)際情況選擇參數(shù)，提交運(yùn)行后得到分析結(jié)果。另外，Galaxian Web Sever的可視化模塊提供Circos Builder功能，提交相關(guān)數(shù)據(jù)文件和配置文件可在線繪制Circos圖[14]。

2 結(jié)果

2.1 抗體基因序列分析平臺(tái)和流程

將測(cè)序結(jié)果用BioEidt軟件處理得到反向互補(bǔ)序列，在搭建的平臺(tái)上對(duì)反向互補(bǔ)序列文件進(jìn)行一系列數(shù)據(jù)處理和序列分析。具體流程為，用本地IgBLAST程序注釋后再用MakeDb.py進(jìn)行格式轉(zhuǎn)換，生成標(biāo)準(zhǔn)tab-delimited數(shù)據(jù)集文件，用countGene函數(shù)進(jìn)行VDJ基因使用統(tǒng)計(jì)分析、De?fineClones程序進(jìn)行抗體序列聚類分析、dnapars程序繪制抗體譜系樹、createTargetingModel函數(shù)構(gòu)建體細(xì)胞高頻突變模式模型、aminoAcidProperties函數(shù)分析CDR3區(qū)氨基酸理化性質(zhì)，完成序列分析。圖1顯示了所搭建的抗體基因序列分析平臺(tái)和抗體基因序列分析的流程。

2.2 V（D）J基因使用

通過(guò)R語(yǔ)言包alakazam中的countGene函數(shù)統(tǒng)計(jì)抗體重鏈VDJ基因使用情況（圖2），包括節(jié)段基因家族、節(jié)段基因和等位基因的統(tǒng)計(jì)結(jié)果。結(jié)果顯示VDJ基因使用頻率最高的胚系基因節(jié)段家族分別為IGHV3、IGHD3、IGHJ4，使用頻率最高的胚系基因分別為IGHV3-23、IGHD3-10，使用頻率最高的等位基因分別為IGHV3-23*01、IGHD6-13*01、IGHJ4*02。

圖1 抗體基因序列分析流程

圖2 VDJ基因的部分統(tǒng)計(jì)結(jié)果A：VDJ基因家族使用統(tǒng)計(jì)；B：VDJ基因使用統(tǒng)計(jì)；C：VDJ等位基因使用統(tǒng)計(jì)

另外，ARGalaxy Web Sever生成的VDJ基因使用情況的柱狀統(tǒng)計(jì)圖（圖3）顯示使用頻率最高的胚系基因節(jié)段家族分別為IGHV3、IGHD3、IGHJ4，使用頻率最高的胚系基因分別為IGHV3-23、IGHD3-10，與平臺(tái)分析得到的結(jié)果一致。

圖3 VDJ基因使用頻率統(tǒng)計(jì)圖A：V基因家族使用統(tǒng)計(jì)；B：D基因家族使用統(tǒng)計(jì)；C：J基因家族使用統(tǒng)計(jì)；D：V基因使用統(tǒng)計(jì)；E：D基因使用統(tǒng)計(jì)

ARGalaxy Web Sever生成的Circos圖和Heat?map圖（圖4）反映抗體重鏈VDJ基因的組合情況。圖4A、D表現(xiàn)D基因和J基因的組合，B、E表現(xiàn)V基因和D基因的組合，C、F表現(xiàn)V基因和J基因的組合。Circos圖圓周上的每段弧代表特定基因，2段弧之間的連線（連接條帶）即表示2個(gè)基因的組合情況，連接條帶越寬則表明在統(tǒng)計(jì)結(jié)果中該組合情況出現(xiàn)的頻率越高。Heatmap圖的橫軸和縱軸分別表示2類節(jié)段基因，圖中的方格點(diǎn)表征特定基因的組合情況，方格點(diǎn)的顏色越深（偏藍(lán)色）表明在統(tǒng)計(jì)結(jié)果中該組合出現(xiàn)的頻率越高。從Circos和Heatmap圖可以看出，DJ基因組合出現(xiàn)頻率較高的是IGHD6-13和IGHJ4、IGHD-10和IGHJ4、IGHD3-10和IGHJ6；VD基因組合出現(xiàn)頻率較高的是IGHV3-64和IGHD6-13、IGHV3-23和IGHD4-24、IGHV3-23和IGHD3-10；VJ基因組合出現(xiàn)頻率較高的是IGHV3-23和IGHJ4、IGHV3-64和IGHJ4、IGHV3-23和IGHJ6。

圖4 VDJ基因組合情況A：D基因和J基因的組合統(tǒng)計(jì)Circos圖；B：V基因和D基因的組合統(tǒng)計(jì)Circos圖；C：D基因和J基因的組合統(tǒng)計(jì)Circos圖；D：D基因和J基因的組合統(tǒng)計(jì)Heatmap圖；E：V基因和D基因的組合統(tǒng)計(jì)Heatmap圖；F：D基因和J基因的組合統(tǒng)計(jì)Heatmap圖

2.3 Clone聚類分析

2.3.1 克隆組的定義 基于抗體基因序列間的漢明距離（Hamming distance）對(duì)序列進(jìn)行聚類，定義克隆組，將每條序列都?xì)w入特定的“克隆組”。根據(jù)該模型的定義，“克隆組”內(nèi)的抗體序列可能源自相同的B細(xì)胞。利用DefineClones程序?qū)λ芯康男蛄袛?shù)據(jù)集劃分了254個(gè)克隆組，同一克隆組內(nèi)的序列數(shù)目從1條至26條不等。定義只包含1條序列的克隆組為單序列克隆組，則有中和活性抗體的4條重鏈序列均分布在單序列克隆組中；僅有結(jié)合活性抗體的總共30條重鏈序列中，22條分布在單序列克隆組中，7條集中分布在239號(hào)克隆組中，1條分布在36號(hào)克隆組中。圖5顯示了克隆組內(nèi)序列數(shù)量最多的前15個(gè)克隆組。

2.3.2 生成B細(xì)胞譜系樹（lineage trees） PHY?LIP軟件包中的dnapars程序?yàn)?39號(hào)“克隆組”內(nèi)的抗體基因序列生成抗體譜系樹（圖6），結(jié)合抗體基因所在的節(jié)點(diǎn)位置用紅圈標(biāo)記，根節(jié)點(diǎn)是胚系基因。dnapars程序基于最大簡(jiǎn)約法（maximum parsimony）建樹，該算法常用于系統(tǒng)發(fā)生學(xué)分析，根據(jù)離散型性狀（如基因序列）的變異程度，構(gòu)建系統(tǒng)發(fā)育樹和分析演化關(guān)系。

2.4 重鏈CDR3區(qū)物理化學(xué)性特性分析

考慮到抗體重鏈的CDR3區(qū)在抗原抗體相互作用中的重要性，統(tǒng)計(jì)分析了其氨基酸序列的理化特性，包括CDR3區(qū)長(zhǎng)度、氨基酸使用情況等。圖7顯示中和抗體CDR3區(qū)較其他抗體對(duì)堿性和酸性氨基酸的使用頻率較低，偏向于使用疏水性氨基酸。

2.5 體細(xì)胞高頻突變模型

圖5 包含序列數(shù)量最多的15個(gè)克隆組

圖6 第239號(hào)“克隆組”的抗體譜系樹

針對(duì)給定的抗體基因序列數(shù)據(jù)集，構(gòu)建出5個(gè)堿基長(zhǎng)度motif模式的體細(xì)胞高頻突變模型（SHM targeting models），反映特定突變發(fā)生的可能性。模型結(jié)果以猬形圖的形式表現(xiàn)，猬形圖主要由5層圓環(huán)和外圍放射狀線條構(gòu)成，第3層圓環(huán)表示特定的堿基（A、T、G、C中的某一種），內(nèi)外側(cè)4個(gè)圓環(huán)上依次循環(huán)排布著ATGC等4種堿基（顏色區(qū)別）。因此，A圖表示NNANN模式motifs，B圖表示NNTNN模式motifs，C圖表示NNCNN模式motifs，D圖表示NNGNN模式motifs，其中N表示4種堿基中的某一種。沿固定半徑從內(nèi)圈到外圈即可讀取到特定的5個(gè)堿基長(zhǎng)度的motif，其對(duì)應(yīng)的放射狀線條的長(zhǎng)度表示發(fā)生突變的可能性。放射狀線條的不同顏色分別表示傳統(tǒng)的hotspots（熱點(diǎn)）和cold-spots（冷點(diǎn)），hot-spots標(biāo)記為紅色（WRC/GYW）和綠色（WA/TW），cold-spots標(biāo)記為藍(lán)色（SYC/GRS）。從圖8可以看出，抗體的體細(xì)胞發(fā)生高頻突變更可能發(fā)生在AAAAC、GT?TAA、AGCTC、AAGTA等熱點(diǎn)motif。

3 討論

前期通過(guò)采集重組PA疫苗免疫志愿者血液樣品，并進(jìn)行流式細(xì)胞分選和單細(xì)胞PCR快速擴(kuò)增出全人源單克隆抗體基因，對(duì)于測(cè)序后得到大量的抗體可變區(qū)基因序列，采用簡(jiǎn)單的序列比對(duì)和基礎(chǔ)的統(tǒng)計(jì)方法已經(jīng)遠(yuǎn)遠(yuǎn)不能滿足批量序列分析和生物學(xué)信息提取的需要。本研究采取了一系列行之有效的生物信息學(xué)手段處理和研究序列信息，取得較好的分析結(jié)果，這些流程和方法也同樣適用于更高通量的抗體基因序列的處理和分析。

抗體重鏈可變區(qū)基因包括V基因片段（vari?able segment）、D基因片段（diversity segment）和J基因片段（joining segment）。對(duì)VDJ基因使用情況的統(tǒng)計(jì)顯示，對(duì)于IgHV基因，使用頻率最高的胚系基因節(jié)段家族為IgHV3，其次是IgHV4。僅發(fā)現(xiàn)IgHV1、IgHV3、IgHV4、IgHV5等4個(gè)胚系基因節(jié)段家族的使用，這可能與單細(xì)胞PCR引物的偏性或抗體序列有限的樣本量有關(guān)。對(duì)于IgHD基因，使用頻率最高的胚系基因節(jié)段家族為IgHD3。對(duì)于IgHJ基因，使用頻率最高的胚系基?因節(jié)段家族為IgHJ4。這些結(jié)果與預(yù)期及以前的報(bào)道相符[15-16]。VDJ基因組合方式的頻率也有較大差別，IgHV3-23和IgHJ4、IgHV3-64和IgHD6-13、IgHD6-13和IgHJ4組合出現(xiàn)在同一條抗體重鏈基因上的可能性最大?？梢暬慕y(tǒng)計(jì)結(jié)果同時(shí)直觀地表明了抗體基因的組合多樣性，VDJ基因在B細(xì)胞發(fā)育過(guò)程中發(fā)生重排，從而形成多種多樣的組合形式，產(chǎn)生豐富的多樣性。輕鏈基因結(jié)構(gòu)與此類似，只是缺少D基因片段。

圖7 重鏈CDR3區(qū)氨基酸序列理化特性

對(duì)克隆組的定義和序列的聚類分析發(fā)現(xiàn)，23.3%的PA結(jié)合抗體都出現(xiàn)在同一個(gè)抗體聚類（239號(hào)克隆組）中，而具有中和活性的重鏈可變區(qū)基因序列更傾向于分布在單序列克隆組（克隆組內(nèi)只有一條序列樣本），即中和抗體的重鏈基因序列比較“獨(dú)特”，與數(shù)據(jù)集中的其他序列具有演化關(guān)系的可能性較小。同時(shí)，中和抗體CDR3區(qū)較其他抗體對(duì)堿性和酸性氨基酸的使用頻率較低，偏向于使用疏水性氨基酸，在炭疽中和抗體親和力成熟改造和中和活性提高的過(guò)程中，可以參考上述結(jié)果，對(duì)可變區(qū)的改造朝著特定聚類方向，或者使用疏水性氨基酸從而提高抗體的親和力和中和活性。

圖8 SHM模型猬形圖A：NNANN模式motifs；B：NNTNN模式motifs；C：NNCNN模式motifs；D：NNGNN模式motifs

體細(xì)胞高頻突變也是抗體多樣性的重要來(lái)源，其在抗原刺激后作用于已重排過(guò)的V基因上。其實(shí)質(zhì)是抗體基因（尤其是在突變熱點(diǎn)位置，主要是CDR區(qū)）發(fā)生的點(diǎn)突變，而且發(fā)生頻率很高。突變后有些分子的親和力會(huì)優(yōu)于原先的分子，從而發(fā)生親和力成熟的現(xiàn)象。在統(tǒng)計(jì)抗體基因序列數(shù)據(jù)集信息基礎(chǔ)上，構(gòu)建出體細(xì)胞高頻突變模式模型，模型可視化表現(xiàn)為猬形圖的形式，直觀反映5個(gè)堿基長(zhǎng)度的motif的突變可能性以及hot-spots和cold-spots的分布情況，對(duì)于抗體改造靶點(diǎn)的確定有一定的指導(dǎo)意義。如果能使用更高通量的抗體組庫(kù)序列信息，并且做對(duì)照分析（如比較免疫前后的抗體基因變化，比較不同物種之間抗體基因的差異等），則構(gòu)建的細(xì)胞高頻突變模式模型更具科學(xué)性和實(shí)際意義。本研究結(jié)果提供了一種將DNA和蛋白質(zhì)的功能、結(jié)構(gòu)和進(jìn)化信息聯(lián)系起來(lái)的生物信息學(xué)手段，可為其他疾病的相關(guān)研究提供參考。

[1]Robins H.Immunosequencing:applications of immune repertoiredeep sequencing[J].CurrOpin Immunol, 2013,25(5):646-652.

[2] Stern J N,Yaari G,Vander Heiden J A,et al.B cells populating the multiple sclerosis brain mature in the draining cervical lymph nodes[J].Sci Transl Med, 2014,6(248):248ra107.

[3] Khurana S,Fuentes S,Coyle E M,et al.Human anti?body repertoire after VSV-Ebola vaccination identifies noveltargets and virus-neutralizing IgM antibodies [J].Nat Med,2016,22(12):1439-1447.

[4] Wang B,Kluwe C A,Lungu O I,et al.Facile discov?ery of a diverse panel of anti-Ebola virus antibodies by immune repertoire mining[J].SciRep,2015,5: 13926.

[5] Gupta N T, Heiden J A V, Uduman M, et al. Change-O:a toolkit for analyzing large-scale B cell immunoglobulin repertoire sequencing data[J].Bioinfor?matics,2015,31(20):3356.

[6] Ijspeert H,van Schouwenburg P A,Van Z D,et al. Antigen receptor galaxy:a user-friendly,Web-based tool for analysis and visualization of T and B cell re?ceptor repertoire data[J].J Immunol,2017,198(10):415 6-4165.

[7] Moorhouse M J.ImmunoGlobulin galaxy(IGGalaxy)for simple determination and quantitation of immunoglobu?lin heavy chain rearrangements from NGS[J].BMC Im?munol,2014,15(1):59.

[8] Hall T A.BioEdit:a user-friendly biological sequence alignment editor and analysis program for Windows 95/ 98/NT[J].Nucleic Acids Symp Ser,1999,41(41):95-98.

[9] Ye J,Ma N,Madden T L,et al.IgBLAST:an immu?noglobulin variabledomain sequenceanalysistool[J]. Nucleic Acids Res,2013,41(Web Server issue):W34.

[10]Stern J N,Yaari G,Vander Heiden J A,et al.B cells populating the multiple sclerosis brain mature in the draining cervical lymph nodes[J].Sci Transl Med, 2014,6(248):248ra107.

[11]Glanville J,Zhai W,Berka J,et al.Precise determina?tion of the diversity of a combinatorial antibody li?brary givesinsightinto the human immunoglobulin repertoire[J].Proc Natl Acad Sci USA,2009,106(48): 20216-20221.

[12]Felsenstein J.PHYLIP-phylogenetic inference package (Version 3.2)[J].Cladistics,1989,5(4):164-166.

[13]Gur Y,Vander H J A,Mohamed U,et al.Models of somatic hypermutation targeting and substitution based on synonymous mutations from high-throughput immu?noglobulin sequencing data[J].Front Immunol,2013,4: 358.

[14]Krzywinski M,Schein J,Birol I,et al.Circos:an in?formation aesthetic forcomparative genomics[J].Ge?nome Res,2009,19(9):1639-1645.

[15]Laserson U,Vigneault F,Gadalamaria D,et al.Highresolution antibody dynamics of vaccine-induced im?mune responses[J].Proc Natl Acad Sci USA,2014,111 (13):4928-4933.

[16]Mroczek E S,Ippolito G C,Rogosch T,et al.Differ?ences in the composition of the human antibody reper?toire by B Cell subsets in the blood[J].Front Immu?nol,2014,5(5):96.