劉渝嬌，陳鄭珊，王美榮，張冠英，范鵬飛，房婷，遲象陽，于長明

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

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

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

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

1 材料與方法

1.1 抗體基因序列來源

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

1.2 基因序列的注釋和標準數(shù)據(jù)集的生成

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

1.3 基于抗體基因標準數(shù)據(jù)集的分析方法

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

2 結(jié)果

2.1 抗體基因序列分析平臺和流程

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

2.2 V（D）J基因使用

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

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

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

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

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

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個基因的組合情況，連接條帶越寬則表明在統(tǒng)計結(jié)果中該組合情況出現(xiàn)的頻率越高。Heatmap圖的橫軸和縱軸分別表示2類節(jié)段基因，圖中的方格點表征特定基因的組合情況，方格點的顏色越深（偏藍色）表明在統(tǒng)計結(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)計Circos圖；B：V基因和D基因的組合統(tǒng)計Circos圖；C：D基因和J基因的組合統(tǒng)計Circos圖；D：D基因和J基因的組合統(tǒng)計Heatmap圖；E：V基因和D基因的組合統(tǒng)計Heatmap圖；F：D基因和J基因的組合統(tǒng)計Heatmap圖

2.3 Clone聚類分析

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

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

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

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

2.5 體細胞高頻突變模型

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

圖6 第239號“克隆組”的抗體譜系樹

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

3 討論

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

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

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

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

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

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