基于二代高通量測(cè)序技術(shù)的新型冠狀肺炎病毒S朊突變點(diǎn)分析

徐朋朋， 李 爽 ， 金德才， 萬(wàn)云洋

(1.南京普維康生物科技有限公司，江蘇 南京 210017；2.中國(guó)石油大學(xué)(北京) 油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室地層微生物資源與應(yīng)用研究中心 非常規(guī)油氣科學(xué)技術(shù)研究院，北京 102249；3.華北理工大學(xué) 實(shí)驗(yàn)動(dòng)物中心，河北 唐山 063210；4.中國(guó)科學(xué)院 生態(tài)環(huán)境研究中心，北京 100085)

新型冠狀肺炎病毒或新型冠狀病毒(SARS-CoV-2，簡(jiǎn)稱新冠病毒)導(dǎo)致的新型冠狀病毒肺炎(COVID-19，簡(jiǎn)稱新冠肺炎)短短幾個(gè)月內(nèi)在全球爆發(fā)，截止2020年5月，累計(jì)確診已從400萬(wàn)上升到近2 000萬(wàn)例，死亡已從約30萬(wàn)上升到超70萬(wàn)人，且目前仍未有全球消退趨勢(shì)。相關(guān)疫苗研制是當(dāng)務(wù)之急，而對(duì)于免疫相關(guān)較為重要的是其基因組中的刺突糖蛋白(Spike glycoprotein，簡(jiǎn)稱S朊)編碼基因[1]。該基因已被明確報(bào)道有突變點(diǎn)(D614G)[2]，而隨著病毒繼續(xù)流行，不排除有類似D614G的新突變點(diǎn)產(chǎn)生，這對(duì)疫苗效價(jià)會(huì)有不確定影響。從2020年初中國(guó)科學(xué)家[3]公布第一條新冠病毒全基因組序列開始，幾個(gè)月內(nèi)，該病毒已發(fā)生抗原漂變而出現(xiàn)了至少6個(gè)不同基因型[4]。雖然目前確定的是該病毒全基因組整體序列變異度遠(yuǎn)遠(yuǎn)低于流感病毒[5]，但不能排除其重要的功能位點(diǎn)或片段發(fā)生突變的可能。截止目前，已報(bào)道在全球流行SARS-CoV-2毒株S朊區(qū)已經(jīng)有14個(gè)位點(diǎn)發(fā)生了突變[5-6]。新冠病毒侵入人體后，單個(gè)病毒S基因在宿主體內(nèi)不斷復(fù)制，并在復(fù)制過程中發(fā)生隨機(jī)變異。通過高深度測(cè)序，判讀這些突變并研究其可能的突變規(guī)律，是目前可行的評(píng)估手段，目前尚未見國(guó)內(nèi)外有相關(guān)研究。不同于一代和三代測(cè)序技術(shù)，二代測(cè)序技術(shù)雖然讀長(zhǎng)短，但有著高準(zhǔn)確率Q30值和數(shù)據(jù)量大且測(cè)序覆蓋倍數(shù)高的優(yōu)勢(shì)[7]。經(jīng)過多年發(fā)展，二代高通量測(cè)序和數(shù)據(jù)處理方法已經(jīng)逐漸趨于完善和穩(wěn)定，測(cè)序成本和時(shí)間也已降至合理水平。從目前報(bào)道的新冠肺炎數(shù)據(jù)來看，大多采用此種測(cè)序手段。本研究采用2020年3至4月份美國(guó)華盛頓州爆發(fā)的新冠病毒二代高通量測(cè)序數(shù)據(jù)作為原始數(shù)據(jù)，共對(duì)100例病人數(shù)據(jù)進(jìn)行了糖朊編碼區(qū)突變點(diǎn)分析。

1 材料與方法

1.1 材料

1.1.1 數(shù)據(jù)準(zhǔn)備 從NCBI下載130例美國(guó)華盛頓洲公布的從急性感染新冠肺炎患者采集得到的二代高通量測(cè)序數(shù)據(jù)(SRR11621803～SRR11621932)，作為分析準(zhǔn)備數(shù)據(jù)。

1.1.2 分析相關(guān)軟件 數(shù)據(jù)質(zhì)控采用fastp軟件[8]完成(參數(shù)默認(rèn))，序列組裝使用Spades軟件[9]，基因預(yù)測(cè)使用Prokka軟件[10]，序列比對(duì)和格式轉(zhuǎn)換使用bowtie2，samtools軟件，突變點(diǎn)頻率計(jì)算使用VirVarSeq軟件[11]，統(tǒng)計(jì)和繪圖使用R語(yǔ)言完成。

1.2 方法

1.2.1 序列上游分析 將130份二代原始fastq文件(原始數(shù)據(jù)為已剔除宿主的測(cè)序序列，共92.804 M reads)，使用fastp軟件進(jìn)行質(zhì)控，共得到91.224 M reads。使用spades軟件進(jìn)行序列組裝，共得到93份完整新冠病毒序列。將不能組裝為完成圖的樣本數(shù)據(jù)，采用分層拼接，對(duì)剩余樣本進(jìn)行重新組裝，獲得7個(gè)樣本的完整病毒基因組序列。將成功獲得全基因組序列的100份樣本數(shù)據(jù)作為分析的可用數(shù)據(jù)。將拼接出的序列使用Prokka軟件進(jìn)行基因編碼區(qū)(CDS，Coding Sequence)預(yù)測(cè)，選定WuHan株的S區(qū)片段作為標(biāo)準(zhǔn)參考序列(MN908947.3)進(jìn)行BLAST比對(duì)確認(rèn)，獲得完整的100條新冠病毒S區(qū)核酸及翻譯后的S朊序列。

1.2.2 序列下游分析 將100份樣本質(zhì)控后fastq數(shù)據(jù)，分別同各自獨(dú)立組裝出的S片段全長(zhǎng)序列，使用bowtie軟件進(jìn)行比對(duì)(mapping)，使用Virvarseq軟件以氨基酸為單位，統(tǒng)計(jì)各自樣本的病毒S區(qū)每個(gè)密碼子所有突變類型和突變頻率。同時(shí)，以MN908947.3為參考，將100份樣本中所有突變型與參考序列進(jìn)行比較。將所得數(shù)據(jù)導(dǎo)入R語(yǔ)言中，使用reshape2包，ggplot2包完成統(tǒng)計(jì)和繪圖?？紤]到樣本起始數(shù)據(jù)量差別較大，每個(gè)樣本的病毒平均覆蓋倍數(shù)從30×～1 226×不等。參考流感病毒相關(guān)研究報(bào)道，深度測(cè)序的樣本點(diǎn)突變一般可以選擇1%～5%作為陽(yáng)性突變判定標(biāo)準(zhǔn)[12-13]。由于該病毒突變率遠(yuǎn)低于流感病毒，同時(shí)為了盡可能地舍棄掉測(cè)序或拼接組裝錯(cuò)誤等因素影響，本次分析選取0.15(即15%)作為單個(gè)樣本的單個(gè)密碼子陽(yáng)性突變率的閾值。

2 結(jié)果與分析

2.1 新冠病毒S朊在體內(nèi)隨機(jī)突變結(jié)果統(tǒng)計(jì)

統(tǒng)計(jì)100份新冠病毒S朊各自在其病人體內(nèi)的隨機(jī)突變點(diǎn)，結(jié)果如圖1所示。

圖1 基于美國(guó)華盛頓州公開二代高通量測(cè)序數(shù)據(jù)的S朊突變點(diǎn)分析Fig.1 All spike protein mutations based on public second-generation sequencing data in 100 SRAR-CoV-2 viruses samples collected from Washington, the United StatesS朊全長(zhǎng)1273aa，分為S1和S2區(qū)；重要的結(jié)構(gòu)區(qū)有SP-信號(hào)肽區(qū)(1aa～11aa)，RBD-受體結(jié)合區(qū)( 330aa～583aa)；A：100份樣本整體與參考序列(MN908947.3)每個(gè)氨基酸位相似頻率平均值，該值越低說明100份樣本整體突變率越高，紅色陰影為標(biāo)準(zhǔn)差(SD)；B：100份樣本中所有突變率大于15%的突變型集合(1個(gè)相同類型突變可能出現(xiàn)在不同樣本中，因此可以有多個(gè)記錄，同時(shí)不考慮密碼子搖擺性)，圖中每個(gè)點(diǎn)對(duì)應(yīng)單獨(dú)某個(gè)樣本在此處密碼子的突變率The full length of S protein is 1273aa,which is divided into S1 an S2 regions.The important structural regious are sp signal peptide region (1aa-11aa) and RBD receptor binding region(330aa-583aa);A: Average amino acid mutation frequency of spike protein in 100 samples, compared with reference sequence(MN908947.3). The lower frequency value, the higher mutation rate;B: All the mutation points collections with mutation rates greater than 15% in each of 100 samples (no regard the codon wobble mutation). Each point corresponds to the mutation rate of a single sample at this codon position

從美國(guó)華盛頓州公布二代高通量測(cè)序數(shù)據(jù)分析可見(圖1)，S區(qū)在全長(zhǎng)區(qū)域內(nèi)，整體突變率超過15%的僅有1個(gè)點(diǎn)(為已報(bào)道的D614G點(diǎn)[5])(圖1A)，同時(shí)未發(fā)現(xiàn)D614G在該100例病人中存在有雜合突變情況，這說明對(duì)于這些在3至4月底采集的病例，該突變點(diǎn)可能已在華盛頓疫區(qū)趨向于穩(wěn)定流行。

2.2 新冠病毒S朊低頻率隨機(jī)突變分析

除D614G點(diǎn)外，其余突變點(diǎn)整體突變率雖均<5%，但是有較明顯的連續(xù)分布規(guī)律(圖1B)，為此進(jìn)一步統(tǒng)計(jì)了整體突變率>3%的點(diǎn)(表1)。

表1 與參考序列MN908947.3相比整體突變率>3%的密碼子集合

現(xiàn)有研究已證實(shí)COVID-19同其他冠狀病毒科病毒類似，S朊主要參與介導(dǎo)受體相互作用及免疫識(shí)別，從而刺激機(jī)體產(chǎn)生抗體，因此是開發(fā)疫苗或單抗十分重要的目標(biāo)區(qū)[2]。根據(jù)現(xiàn)有資料[5]，本研究和已報(bào)道的S朊區(qū)預(yù)測(cè)的抗原表位位點(diǎn)做了匹配(圖2)，其中黑色大寫標(biāo)記為本研究找到的突變集中區(qū)域。

由圖2看出，新冠病毒在病人體內(nèi)復(fù)制時(shí)產(chǎn)生的隨機(jī)突變點(diǎn)主要集中在S1區(qū)，且集中在S1區(qū)的前段部分(126aa～204aa)，而RBD區(qū)未見突變位點(diǎn)。從穩(wěn)定性來說，S1中的RBD區(qū)>S2區(qū)>S1其他區(qū)。這與報(bào)道的冠狀病毒科其他病毒特征基本一致。

圖2 S朊突變集中區(qū)預(yù)測(cè)Fig.2 Prediction on high-frequency mutation region of spike protein 綠色-S1區(qū)：預(yù)測(cè)的抗原表位；紫色-S2區(qū)：預(yù)測(cè)的抗原表位；下劃線：預(yù)測(cè)的S1/S2裂解位點(diǎn)；黑色大寫粗體字母：預(yù)測(cè)的突變集中區(qū)(單個(gè)病人體內(nèi)>15%的突變點(diǎn)集合)Green-S1 region: predicted epitope sites; Purple-S2 region: predicted epitope sites; Underline: predicted cleave point of S1/S2 region; Black capital letters in bold: predicted mutation points (frequency > 15% in a single patient)

3 討 論

目前新冠病毒侵染機(jī)理和S朊的三維結(jié)構(gòu)已經(jīng)得到解析[14-15]，研究已取得了較喜人的進(jìn)展。但由于研究時(shí)間有限，其突變的發(fā)生規(guī)律及突變后對(duì)病毒毒力的影響尚未十分清楚。而這些突變對(duì)于疫苗效價(jià)以及核酸診斷試劑準(zhǔn)確率，有時(shí)會(huì)有較大影響。

本研究通過研究100例美國(guó)新冠肺炎病人SARS-CoV-2二代高通量測(cè)序的數(shù)據(jù)，初步確定了S朊的突變區(qū)主要集中在126aa～153aa和194aa～204aa區(qū)域以及S1朊開頭的9個(gè)信號(hào)肽和S2朊末端20aa區(qū)域，屬首次報(bào)道。這些突變集中區(qū)，以單個(gè)樣本散在突變?yōu)橹?，推測(cè)是病毒侵入宿主后復(fù)制過程中隨機(jī)產(chǎn)生的。這些突變點(diǎn)主要集中在預(yù)測(cè)的抗原表位區(qū)，推測(cè)可能和該病毒對(duì)抗機(jī)體免疫有關(guān)。此外，本研究也證實(shí)了S2區(qū)突變少于S1區(qū)，與冠狀病毒科S朊特征一致[16]。

值得注意的是本次組裝結(jié)果顯示，同武漢株(WuHan strain)相比，一共有10個(gè)位點(diǎn)在單個(gè)或多個(gè)樣本中發(fā)生了改變(僅有1個(gè)D614G密碼子突變樣本數(shù)較多，其余變異均為單株散在變異)。與根據(jù)全球流行株上千條序列統(tǒng)計(jì)得到的14突變點(diǎn)(呈一定比例流行的統(tǒng)計(jì)點(diǎn))相比[5]，有兩突變點(diǎn)與本研究一致(D614G，Y145H)。本研究著重于病毒在同一個(gè)病人體內(nèi)復(fù)制的突變規(guī)律研究，這說明該病毒雖然突變率較低，但不排除有類似D614G的新突變點(diǎn)出現(xiàn)并穩(wěn)定流行的可能。

根據(jù)目前的研究顯示，SRA數(shù)據(jù)庫(kù)中編號(hào)為ERR4085412(記錄為英國(guó)2020年3月27號(hào)采集)的新冠病人樣本，其體內(nèi)出現(xiàn)D614G的突變比率為59.6%(數(shù)據(jù)未附)，說明類似D614G這種突變點(diǎn)大流行，很可能是最初來源于某個(gè)或某些群體感染，病毒在其個(gè)體內(nèi)復(fù)制時(shí)，產(chǎn)生了高頻的D614G點(diǎn)突變。

國(guó)外已有學(xué)者將流感病毒研究方法[17]引入新冠病毒研究中[18]，通過構(gòu)建新冠病毒S朊突變庫(kù)質(zhì)粒，整合進(jìn)酵母表達(dá)系統(tǒng)中，測(cè)定表達(dá)的S朊與ACE親和力，通過超高深度測(cè)序，研究其突變點(diǎn)與病毒毒力之間的關(guān)系，從而篩選出免疫原性最佳的突變點(diǎn)組合，同時(shí)發(fā)現(xiàn)了較D614G結(jié)合能力更強(qiáng)的突變點(diǎn)，但目前尚未在流行株中發(fā)現(xiàn)這些突變點(diǎn)。這種方法較為準(zhǔn)確地闡明了突變點(diǎn)與病毒毒力之間的對(duì)應(yīng)關(guān)系。而我們研究重點(diǎn)在于野毒株在單個(gè)病人體內(nèi)復(fù)制時(shí)的隨機(jī)突變規(guī)律，這對(duì)于實(shí)時(shí)觀測(cè)病毒是否向高毒力突變點(diǎn)變異有著較為重要的意義。

已有報(bào)道表明[19]，可能是地域隔離使病毒被鎖定到一個(gè)區(qū)域后造成特有突變。呼吁國(guó)內(nèi)外學(xué)者更多地公開原始數(shù)據(jù)，基于同一種測(cè)序方法分析更大的數(shù)據(jù)集合，這些突變點(diǎn)在不同的地域和時(shí)間里將如何變化，從而更全面地了解和預(yù)測(cè)病毒突變規(guī)律，這將是下一步的主要工作。