MutPrimerDesign：用于人類基因編碼區(qū)域突變位點的引物設(shè)計程序

曹英豪，彭公信

(中國醫(yī)學(xué)科學(xué)院 基礎(chǔ)醫(yī)學(xué)研究所 & 北京協(xié)和醫(yī)學(xué)院 基礎(chǔ)醫(yī)學(xué)院，北京 100730)

DNA突變不僅是導(dǎo)致生物多樣性和生物進(jìn)化的重要原因之一，同時也與人類的健康和疾病發(fā)生有著密切的關(guān)系。比如引起腫瘤發(fā)生最重要的原癌基因和抑癌基因，其大部分功能與細(xì)胞的生長、分裂、錯配修復(fù)和細(xì)胞凋亡等相關(guān)。如果這兩類基因在蛋白編碼區(qū)發(fā)生突變，則可能導(dǎo)致蛋白質(zhì)功能的改變或者喪失，造成細(xì)胞的非正常生長，從而導(dǎo)致腫瘤的產(chǎn)生。因此，檢測基因突變并尋找可用于診斷和治療的突變位點已經(jīng)成為腫瘤研究中的常規(guī)內(nèi)容，具有重要的臨床和生物學(xué)意義[1]。

發(fā)生在基因編碼區(qū)的DNA突變，根據(jù)其對蛋白產(chǎn)物的影響可以分為移碼突變、非同義突變、同義突變等。其中前兩種突變由于引起了蛋白質(zhì)的序列變化，并可能影響其正常功能的發(fā)揮，因此，對這兩種突變的檢測和研究尤為重要。近年來，以高通量測序技術(shù)為主的檢測技術(shù)得到了迅猛的發(fā)展，比如目標(biāo)區(qū)域捕獲技術(shù)(Panel)、全外顯子組測序技術(shù)(WES)、全基因組測序技術(shù)(WGS)等[2]。尤其是Panel測序技術(shù)，通過特異性捕獲探針對多個基因目標(biāo)區(qū)域DNA片段進(jìn)行捕獲測序，具有成本低、測序深度高的優(yōu)勢，能夠準(zhǔn)確給出重要基因的突變信息，在臨床上得到了廣泛的應(yīng)用[3-5]。

為了應(yīng)用于疾病的診斷和治療，通過測序平臺獲取的或者文獻(xiàn)中已發(fā)表的突變信息，需要進(jìn)一步擴大人群規(guī)模來驗證該突變的可靠性。利用qPCR及衍生技術(shù)檢測突變位點，在降低成本的同時又保證其準(zhǔn)確性，是一種簡單而可行的方法[6-8]。設(shè)計引物是qPCR的首要問題，目前已有許多成熟的引物設(shè)計軟件可以使用。這其中除了Primer3[9]、Oligo[10]等綜合性的引物設(shè)計軟件包之外，還有一些專用的引物設(shè)計軟件，比如整合了Primer3和BLAST的模板特異性引物設(shè)計工具Primer-BLAST[11]，突變引物設(shè)計工具WASP[12]，PrimerX(http://www.bioinformatics.org/primerx/)，MutScreener[13]，Edesign[14]以及適用于克隆表達(dá)載體的PrimerCE[15]等。但由于測序平臺的突變檢測報告以及文獻(xiàn)中給定的突變信息，大部分均是以HGNC(HUGO Gene Nomenclature Committee)[16]的基因名稱(Gene Symbol)命名的基因編碼區(qū)域的相關(guān)坐標(biāo)，因此在設(shè)計引物時，這些軟件或多或少有些不適用。主要有兩個問題：(1)人類基因普遍有多個轉(zhuǎn)錄本，需要通過不同的數(shù)據(jù)庫尋找對應(yīng)的標(biāo)準(zhǔn)轉(zhuǎn)錄本序列。(2)由于內(nèi)含子的影響，對于外顯子編碼區(qū)的坐標(biāo)，需要根據(jù)基因結(jié)構(gòu)進(jìn)行基因組坐標(biāo)的手工轉(zhuǎn)換。另外，為了能夠進(jìn)一步直接從擴增子中檢測序列是否發(fā)生突變，針對突變位點可增加探針序列的設(shè)計。

利用Python語言編寫了MutPrimerDesign，通過對人類基因組序列數(shù)據(jù)庫以及基因注釋信息的解析，可識別多種數(shù)據(jù)庫的基因名稱編號，轉(zhuǎn)換基因編碼區(qū)坐標(biāo)為基因組坐標(biāo)。同時，通過對primer3-py(https://github.com/libnano/primer3-py)程序包的整合，可批量設(shè)計突變位點附近的引物及探針序列。軟件的下載地址為:https://github.com/bioinfo-ibms-pumc/MutPrimerDesign。

1 方 法

MutPrimerDesign的流程軟件主要分為三步，第一步接收輸入信息，并將基因名稱轉(zhuǎn)換為HGNC中標(biāo)準(zhǔn)基因名稱，同時獲取標(biāo)準(zhǔn)轉(zhuǎn)錄本以及相對應(yīng)的RefSeq[17]編號。在檢測輸入數(shù)據(jù)時，MutPrimerDesign可以接收三種不同的基因命名，包括基因名稱、Ensembl編號或者HGNC編號。第二步轉(zhuǎn)換編碼區(qū)坐標(biāo)為基因組坐標(biāo)，并提取突變位點附近基因組序列。在轉(zhuǎn)換坐標(biāo)前，通過對人類基因組的注釋文件進(jìn)行解析，構(gòu)建所有基因?qū)?yīng)的基因結(jié)構(gòu)特征，并根據(jù)該結(jié)構(gòu)，將輸入的標(biāo)準(zhǔn)轉(zhuǎn)錄本的編碼區(qū)坐標(biāo)，轉(zhuǎn)換為對應(yīng)基因組上的坐標(biāo)，并進(jìn)一步利用python的BioSeq程序包對人類基因組序列進(jìn)行解析，提取相應(yīng)的擴增子序列(Amplicon)。第三步調(diào)用python的primer3-py程序包，進(jìn)行引物設(shè)計，輸出報告。為了方便在引物擴增序列的同時能夠區(qū)分野生型和突變型位點，軟件中加入了探針序列的設(shè)計，并默認(rèn)設(shè)定探針序列的Tm值比引物序列的Tm值低5度。具體流程(見圖1)。

1.1 數(shù)據(jù)源下載及安裝

通過git命令可以克隆下載軟件全部代碼，命令如下：

git clone https://github.com/bioinfo-ibms-pumc/MutPrimerDesign.git

下載后的目錄中主要包括4個文件，分別為MutPrimerDesign.py，gffLoader.py，fastaLoc.py和namedb文件。主程序MutPrimerDesign.py用于接收基因突變信息并設(shè)計引物，gffLoader.py用于解析基因注釋文件，fastaLoc.py用于提取基因組序列，namedb文件為基因名稱數(shù)據(jù)庫。MutPrimerDesign依賴BioPython、Pandas、Primer3-Py數(shù)據(jù)包，這些python依賴包均可使用pip命令進(jìn)行標(biāo)準(zhǔn)安裝。依賴包的安裝命令如下：

pip3 install biopython pandas primer3-py

另外，由于HGNC中基因名稱對應(yīng)于RefSeq的標(biāo)準(zhǔn)序列，因此，為了確保序列坐標(biāo)的準(zhǔn)確性，MutPrimerDesign目前只接收由NCBI網(wǎng)址提供的人類基因組序列文件以及相對應(yīng)版本的基因注釋文件。人類GRCh38版本基因組及基因注釋文件下載及解壓縮命令如下：

wget-c ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.39_GRCh38.p13/GCF_000001405.39_GRCh38.p13_genomic.fna.gz

wget-c ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.39_GRCh38.p13/GCF_000001405.39_GRCh38.p13_genomic.gff.gz

gunzip-c GCF_000001405.39_GRCh38.p13_genomic.fna.gz > genome.fa

gunzip-c GCF_000001405.39_GRCh38.p13_genomic.gff.gz > gene.gff

圖1 MutPrimerDesign引物設(shè)計流程圖及原理Fig.1 Flow chart of MutPrimerDesign and its principle

1.2 命令行說明

MutPrimerDesign主要由兩部分命令組成，一部分為程序主要參數(shù)，包括輸入輸出以及各種依賴數(shù)據(jù)庫的定義；另一部分為引物參數(shù)，主要包括引物的Tm值，長度，GC含量等(見圖2)。

圖2 MutPrimerDesign的命令行說明截圖Fig.2 Screen shot for commands of MutPrimerDesign

MutPrimerDesign接收的輸入?yún)?shù)有兩種類型。第一種，是“-i”命令，直接接收基因突變信息。每個基因突變的信息分為三個部分，以冒號分開，內(nèi)容依次分別為基因名稱、突變在該基因的CDS區(qū)域坐標(biāo)、突變長度。不同突變信息以逗號隔開，比如“KRAS:24:2,MET:1124:1,EGFR:2573:1”。另一種是“-f”命令，直接接收文本文件。文件中每一行為一個基因突變信息，內(nèi)容與-i中突變信息格式一致。

MutPrimerDesign除了輸入的基因信息外，還需要三個輸入文件，分別為基因名稱轉(zhuǎn)換文件namedb，基因組序列文件genome.fa以及相對應(yīng)版本的基因注釋文件gene.gff。其中namedb為MutPrimerDesign預(yù)處理后的文件，另外兩個文件可從NCBI網(wǎng)站上下載解壓后得到。

MutPrimerDesign的輸出分為三個文本文件，通過“-o”參數(shù)來控制。主要包括擴增子序列文件，另外還有一個以“.primer”后綴結(jié)尾的引物序列文件，以及一個以“.detail”后綴結(jié)尾的詳細(xì)報告文件。

2 結(jié) 果

2.1 不同數(shù)據(jù)庫來源的引物設(shè)計舉例

經(jīng)統(tǒng)計，50%的黑色素瘤中原癌基因BRAF均發(fā)生了突變，其中90%的突變發(fā)生在第600位纈氨酸的位置上[18]。該位點是一個非常重要的免疫檢查點，位點的突變會影響腫瘤用藥效果。對此，分別使用BRAF基因的基因名稱、Ensembl ID號、HGNC ID號作為輸入，對其V600突變位點進(jìn)行引物設(shè)計，運行命令如下所示：

(1)python3 MutPrimerDesign.py -g genome.fa -G gene.gff -k namedb -o result.fa -i “BRAF:1799:1”

(2)python3 MutPrimerDesign.py -g genome.fa -G gene.gff -k namedb -o result.fa -i “ENSG00000157764:1799:1” -s ensembl_id

(3)python3 MutPrimerDesign.py -g genome.fa -G gene.gff -k namedb -o result.fa -i “1097:1799:1” -s hgnc_id

BRAF基因編碼區(qū)第1 799位引物設(shè)計的部分結(jié)果輸出截圖(見圖3)。輸出內(nèi)容中，每一個擴增子信息以Amplicon開始，后接設(shè)計順序及引物編號，第一行為擴增子編號，第二行包括基因突變位點信息，依次為輸入突變信息、基因、Refseq編號、染色體編號、擴增子在基因組上坐標(biāo)區(qū)間、突變位點在基因組上坐標(biāo)、突變位點在編碼區(qū)坐標(biāo)、參考堿基及長度和引物編號。第三行為擴增子長度。第四行到第七行為primer3-py設(shè)計的引物及探針序列信息。其中，第四行為根據(jù)熱力學(xué)方法計算引物對整體形成二級結(jié)構(gòu)的可能性以及3’末端形成二級結(jié)構(gòu)的可能性。第五行為5’端引物，主要信息包括：引物序列，引物位置，Tm值，GC含量，引物自身形成二級結(jié)構(gòu)的可能性，3’末端形成二級結(jié)構(gòu)的可能性，形成發(fā)夾結(jié)構(gòu)的可能性，3’末端的穩(wěn)定性(delta G, kcal/mol)。同樣，第六行和第七行分別為3’端引物和探針序列信息。第八行為擴增子序列，第九行為引物、探針、突變位點在擴增子中的位置。由圖3中可以看出，無論是采用基因名稱，Ensembl編號或者HGNC編號，軟件均可以完成相同的引物設(shè)計。

2.2 批量引物設(shè)計舉例

Maki-Nevala S等通過對425例非小細(xì)胞型肺癌病人樣本進(jìn)行測序，找到了9個復(fù)發(fā)突變熱點[19]。利用MutPrimerDesign的命令行輸入“-i”和文件輸入“-f”兩種方式進(jìn)行引物設(shè)計，引物及探針設(shè)計結(jié)果(見表1)。運行命令如下所示：

(1)python3 MutPrimerDesign.py -g genome.fa -G gene.gff -k namedb -o lung.fa -i "KRAS:34:2,MET:1124:1,MET:3029:1,EGFR:2573:1,TP53:820:1,PIK3CA:1624:1,EGFR:2235:1,TP53:469:1"

(2)python3 MutPrimerDesign.py -g genome.fa -G gene.gff -k namedb -o lung.fa -f samples.file

2.3 引物重新設(shè)計舉例

對于已經(jīng)完成引物設(shè)計的基因位點，并且擴增子序列已經(jīng)保存在result.fa文件中時，如果需要修改引物參數(shù)進(jìn)行重新設(shè)計，則可以使用-p參數(shù)忽略序列提取過程，僅進(jìn)行引物重新設(shè)計，運行命令如下所示：

(1)python3 MutPrimerDesign.py -g genome.fa -G gene.gff -k namedb -o result.fa -p --primerNum 3 --primer_opt_tm 61 --primer_max_tm 64 --primer_opt_size 23

2.4 軟件運行效率測試

MutPrimerDesign在設(shè)計引物時消耗內(nèi)存及計算資源較小，普通的臺式機或者筆記本環(huán)境下均可以運行。為了測試軟件的運行性能和穩(wěn)定性，以Linux操作系統(tǒng)為例(系統(tǒng)版本：Linux Mint 19.3 Cinnamon；CPU：Intel Core i7-9700；內(nèi)存：24 GB；存儲：1T機械硬盤)，采用MutPrimerDesign軟件默認(rèn)參數(shù)進(jìn)行了性能測試：系統(tǒng)開機后，首次完成8個突變位點的引物設(shè)計耗時為1 m 33 s。而后當(dāng)系統(tǒng)自動加載了基因組信息緩存后，無論是多基因還是單基因的100個突變位點，設(shè)計引物只需要30 s即可完成。

3 結(jié) 論

通過Python，開發(fā)了MutPrimerDesign引物設(shè)計軟件，可專門用于人類基因編碼區(qū)的突變位點引物設(shè)計。通過對基因名稱進(jìn)行解析，MutPrimerDesign將自動選擇HGNC給定的標(biāo)準(zhǔn)轉(zhuǎn)錄本信息，完成基因編碼區(qū)坐標(biāo)到基因組坐標(biāo)的轉(zhuǎn)換以及基因組序列的提取，從而能夠批量、自動化地完成基因編碼區(qū)突變位點的引物及探針序列的設(shè)計。目前軟件的不足之處是輸出格式較為簡單，輸出結(jié)果可以更加豐富全面。

4 討 論

以Primer3, Oligo為首的引物設(shè)計軟件的出現(xiàn)，使得引物設(shè)計相關(guān)工作變得越來越簡單。但是，伴隨著生物技術(shù)的不斷發(fā)展，生物學(xué)問題越來越多，引物設(shè)計的需求也越來越細(xì)化。比如高通量測序平臺的不斷發(fā)展，突變位點檢測等相關(guān)的下游應(yīng)用將會變得越來越多，對此，WASP，MutScreener等特定的引物設(shè)計軟件一定程度下緩解了這些需求。但針對突變位點設(shè)計引物時，用戶輸入基因的名稱可能來自不同數(shù)據(jù)庫，如何自動地轉(zhuǎn)換基因名稱，對含有多個轉(zhuǎn)錄本的基因進(jìn)行標(biāo)準(zhǔn)轉(zhuǎn)錄本的選擇，以及如何根據(jù)基因的編碼區(qū)位置自動定位到基因組上相應(yīng)的位置，這些看似簡單的問題卻沒有完善及考慮。MutPrimerDesign通過對不同數(shù)據(jù)庫的整合以及基因結(jié)構(gòu)的解析解決了這些特定的問題。另外，引物設(shè)計開發(fā)的一個重要指標(biāo)是如何能夠更合理地減少用戶的工作量，輔助用戶完成批量、自動化的引物設(shè)計。MutPrimerDesign試圖通過最簡潔的輸入來完成多種情況下突變位點的引物設(shè)計，也是對該指標(biāo)的一次嘗試。在將來，MutPrimerDesign將會考慮整合HGVS(Human Genome Variation Society)[20]的突變標(biāo)準(zhǔn)規(guī)則，進(jìn)一步直接對接測序平臺的突變檢測報告，豐富并完善突變序列與參考序列的對比信息，更加方便、合理地完成突變位點檢測的引物設(shè)計工作。另外，MutPrimerDesign將會考慮應(yīng)用到更多的物種，實現(xiàn)多物種的突變引物設(shè)計。