曾 永,陳建威,崔恩海,胡國(guó)海,王曉光,季巧英,叢 迎,楊 陽(yáng),楊偉斌
(1. 青田縣疾控中心,青田 323900;2. 青島華大基因研究院,青島 266555;3. 湖州市中心醫(yī)院,湖州 313000;4. 深圳國(guó)家基因庫(kù),深圳 518083;5. 麗水市疾控中心,麗水 323000;6. 麗水市中心醫(yī)院,麗水 323000)
2019年12 月,湖北省武漢市出現(xiàn)嚴(yán)重急性呼吸綜合征冠狀病毒2(severe acute respiratory syndrome coronavirus 2,SARS-CoV-2)引起的新型冠狀病毒肺炎(corona virus disease 2019,COVID-19)[1-2]。病毒感染隨后蔓延到全國(guó)各地,導(dǎo)致嚴(yán)重的呼吸道疾病,嚴(yán)重危害人們的生命健康[3]。除中國(guó)外,COVID-19疫情在全球多個(gè)國(guó)家也出現(xiàn)了,成為國(guó)際關(guān)注的突發(fā)公共衛(wèi)生事件[4-5]。同時(shí),中國(guó)開始出現(xiàn)海外輸入性病例,防控病毒的戰(zhàn)斗還在繼續(xù)。
SARS-CoV-2在人群中的持續(xù)傳播可能會(huì)產(chǎn)生適應(yīng)性突變,因此極有可能會(huì)累計(jì)基因突變,進(jìn)而抵抗人類的免疫系統(tǒng)[6]。已有研究表明,基因突變可能會(huì)導(dǎo)致新型冠狀病毒的毒力、傳染性和傳播性的改變[7-8]。因此,對(duì)不同來(lái)源的新冠肺炎病例進(jìn)行基因組分析有助于了解SARS-CoV-2的變異與感染病程的相關(guān)性,從而為疾病預(yù)防和控制提供建議。
鑒于患者咽拭子和痰液中RNA含量的差異,我們采用最新研發(fā)的捕獲測(cè)序[9-10]和宏轉(zhuǎn)錄組測(cè)序方法[10],對(duì)意大利歸國(guó)和國(guó)內(nèi)的新冠肺炎確診病例進(jìn)行基因組突變分析。與宏轉(zhuǎn)錄組測(cè)序相比,捕獲測(cè)序所需數(shù)據(jù)量?jī)H為宏轉(zhuǎn)錄組測(cè)序的1/10或更少,且測(cè)序結(jié)果準(zhǔn)確。對(duì)病毒基因組的突變分析為深入研究意大利歸國(guó)患者SARS-CoV-2的致病力、傳染性以及國(guó)內(nèi)毒株的差異打下了基礎(chǔ)。
本項(xiàng)目涉及的所有臨床標(biāo)本(包括咽拭子和痰液)均來(lái)自浙江省新冠肺炎確診的19例患者,其中7例為國(guó)內(nèi)病例,12例為歐洲輸入性病例(11例為意大利輸入性病例,1例為西班牙輸入性病例)。在11名意大利患者中有7名是意大利貝加莫同一家餐廳的工作人員(詳細(xì)臨床資料見表1)。該項(xiàng)倫理審查已獲麗水市中心醫(yī)院和深圳華大生命科學(xué)研究院倫理委員會(huì)(institutional review board,IRB)批準(zhǔn)。所有參與者在納入這項(xiàng)研究之前都簽署了書面知情同意書。
通過Maxwell?16 Total RNA純化試劑盒進(jìn)行SARS-CoV-2 RNA提取。隨后針對(duì)SARS-CoV-2的RdRp基因和N基因進(jìn)行實(shí)時(shí)逆轉(zhuǎn)錄聚合酶鏈反應(yīng)(quantitative real-time PCR,qRT-PCR)。采用HCoV-19核酸檢測(cè)試劑盒檢測(cè)和定量檢測(cè)臨床標(biāo)本中的病毒RNA的載量(捷諾生物,上海)。
宏轉(zhuǎn)錄組測(cè)序技術(shù)可及時(shí)快速地對(duì)傳染病進(jìn)行檢測(cè)、防控和變異分析[11-12]。首先,我們對(duì)每個(gè)RNA樣本用DNase I去除DNA,并對(duì)純化的產(chǎn)物進(jìn)行定量(Qubit RNA HS Assay Kit, Thermo Fisher Science, Waltham, MA, USA);然后將純化的RNA在87℃緩沖液中孵育6 min進(jìn)行打斷,隨后將RNA片段利用隨機(jī)引物合成雙鏈(ds)cDNA并進(jìn)行末端修復(fù)、接頭連接;最后對(duì)(ds)cDNA進(jìn)行18個(gè)循環(huán)的PCR擴(kuò)增并加上雙端測(cè)序標(biāo)簽(unique dual-indexed adapters,UDI),使用MGIEasy RNA(華大智造,深圳)的試劑盒經(jīng)過環(huán)化、滾環(huán)復(fù)制后構(gòu)建DNA納米球(DNA nanoball,DNB)文庫(kù)[10]。使用Qubit 3.0 熒光定量?jī)x (Invitrogen, Carlsbad, CA, USA)檢測(cè)文庫(kù)濃度,并使用Agilent 2100 生物分析儀檢測(cè)文庫(kù)片段大小,對(duì)合格文庫(kù)采用DNBSEQ T1測(cè)序平臺(tái)(華大智造,深圳)或MGISEQ-2000測(cè)序平臺(tái)(華大智造,深圳)進(jìn)行雙端100 bp 讀長(zhǎng)模式的高通量測(cè)序[4]。
病毒全基因組雜交捕獲技術(shù)是病毒分離鑒定的常用方法之一,可提升檢測(cè)靈敏度和病毒基因組覆蓋度[13-14]。利用MGIEasy外顯子捕獲輔助試劑盒(華大智造,深圳)替代了SARS-CoV-2 DNA/RNA捕獲試劑盒(伯科生物,北京)中的阻滯劑寡核苷酸(blocker oligos)和PCR引物寡核苷酸(PCR primer oligos)后,采用基于混合捕獲的富集方法從宏轉(zhuǎn)錄組逆轉(zhuǎn)錄雙鏈cDNA片段中富集SARS-CoV-2特異性片段,并根據(jù)試劑盒說(shuō)明對(duì)SARS-CoV-2病毒全基因組進(jìn)行擴(kuò)增[10]。給富集擴(kuò)增得到的病毒cDNA片段加上雙端測(cè)序標(biāo)簽,經(jīng)過環(huán)化、滾環(huán)復(fù)制后構(gòu)建DNB文庫(kù),對(duì)每個(gè)庫(kù)檢合格的文庫(kù)使用MGISEQ-2000測(cè)序平臺(tái)(華大智造,深圳)進(jìn)行雙端100 bp 讀長(zhǎng)模式的高通量測(cè)序[4]。
為了獲得高質(zhì)量的SARS-CoV-2病毒測(cè)序序列,我們采取了4個(gè)嚴(yán)格的質(zhì)量控制步驟,從每個(gè)宏轉(zhuǎn)錄組文庫(kù)和雜交捕獲文庫(kù)測(cè)序結(jié)果中獲得高質(zhì)量的SARS-CoV-2片段序列。首先用Kraken(v0.10.5)軟件[15]將原始測(cè)序數(shù)據(jù)同SARS-CoV-2參考基因組(GISAID獲取號(hào):EPI_ISL_402119)[16]數(shù)據(jù)進(jìn)行對(duì)比,篩選得到候選的SARS-CoV-2原始測(cè)序數(shù)據(jù);然后依次用fastp(v0.19.5,參數(shù):-q20-u20-n1-150)、SOAPnuke(v1.5.6,參數(shù):-l20-q0.2-E50-d)和PRINSEQ(v0.20.4,參數(shù):LC_Method Dust-LC_Threshold 7)對(duì)測(cè)序數(shù)據(jù)的低質(zhì)量數(shù)據(jù)、接頭污染、PCR重復(fù)測(cè)序序列和低復(fù)雜度序列進(jìn)行過濾。如果該文庫(kù)最終高質(zhì)量有效測(cè)序序列數(shù)小于1 000條,則舍棄其測(cè)序數(shù)據(jù),不用于下游分析。針對(duì)每個(gè)樣本,本文將兩個(gè)不同建庫(kù)方法獲得的大規(guī)模并行高通量測(cè)序(massively parallel sequencing,MPS)數(shù)據(jù)合并用于變異檢測(cè)、組裝和其他下游分析。使用Samtools (v1.9)計(jì)算SARS-CoV-2參考基因組在所有文庫(kù)測(cè)序數(shù)據(jù)中每個(gè)核苷酸位點(diǎn)的覆蓋深度。對(duì)每個(gè)文庫(kù)的測(cè)序數(shù)據(jù)以200 bp的滑動(dòng)窗口用R(v3.6.1)軟件中的ggplot2包繪制SARS-CoV-2參考基因組的測(cè)序覆蓋率和測(cè)序深度。
對(duì)于參考覆蓋深度大于15 X(X表示全基因組的測(cè)序覆蓋深度的層數(shù))的樣本,截取最多100 X的SARS-CoV-2測(cè)序數(shù)據(jù),使用SPAdes(v3.14.0)軟件[17]默認(rèn)參數(shù)進(jìn)行組裝得到完整的病毒基因序列。對(duì)于其他低深度樣本,使用IDBA_hyrid(v1.1.3,參數(shù)為“-mink 21--maxk 91--step 10”)進(jìn)行病毒基因組組裝[18]。
為了檢測(cè)每個(gè)病毒樣品的單核苷酸變異(single nucleotide variations,SNVs)信息,首先使用BWA aln(v0.7.16)將每個(gè)文庫(kù)的有效測(cè)序數(shù)據(jù)與SARSCoV2參考基因組(EPI_ISL_402119)進(jìn)行比對(duì),然后使用Picard (v2.10.10)去除重復(fù)的比對(duì)結(jié)果。基于去重的比對(duì)結(jié)果使用FreeBayes (v1.3.1,參數(shù):-p 1-q 20-m 60-V)[19]分析得到每個(gè)測(cè)序數(shù)據(jù)的SNVs,最后使用snippy-vcf_filter(參數(shù):--minqual 100--minfrac 0.8;https://github.com/tseemann/snippy)過濾低置信度SNVs。剩余的高置信度SNVs在SARS-CoV-2參考基因組中用SNVeff(v4.3)[20]的默認(rèn)參數(shù)進(jìn)行注釋。如果測(cè)序數(shù)據(jù)的參考基因組覆蓋深度大于15 X,則所有SNVs變異檢測(cè)調(diào)用軟件的最小覆蓋參數(shù)設(shè)置為10 X;如果測(cè)序覆蓋深度小于15 X,則設(shè)置為4 X。為了評(píng)估每個(gè)樣本的SNVs準(zhǔn)確性,使用pysamstats (v1.1.2)軟件(https://github. com/alimanfoo/pysamstats)和Pilon(v1.23)軟件[21]進(jìn)行變異檢測(cè),并保留至少兩種方法檢測(cè)到的SNVs。對(duì)于高深度測(cè)序(>40 X)的樣本,對(duì)BWA軟件比對(duì)結(jié)果使用與Ni等[22]相同方法及參數(shù)進(jìn)行宿主內(nèi)單核苷酸變異(intradividual single nucleotide variations,iSNVs)檢測(cè)。
我們通過qRT-PCR檢測(cè)19個(gè)SARS-CoV-2樣本的Ct值,結(jié)果表明所有樣本的Ct值均在22.43~38.84之間(表1)。其中僅3例樣本(C06、F03和F12)檢測(cè)到高病毒拷貝(Ct<24.5),6例樣本(C05、F05、F04、C01、C07和C02)檢測(cè)到低病毒拷貝(24.5<Ct<28.7),其余10例樣本檢測(cè)到極低病毒拷貝(Ct>28.7)。
表1 19個(gè)新冠病例樣本臨床信息Tab. 1 Detailed clinical information of 19 cases
對(duì)所有SARS-CoV-2 RNA樣本構(gòu)建逆轉(zhuǎn)錄文庫(kù),并對(duì)其全基因組進(jìn)行了測(cè)序。這些樣品獲得的原始測(cè)序數(shù)據(jù)雙端序列條數(shù)在8 995萬(wàn)條到140 824萬(wàn)條之間,平均約為64 968萬(wàn)條(130 Gb)。用Kraken和BWA對(duì)SARS-CoV-2病毒測(cè)序數(shù)據(jù)進(jìn)行分類,并計(jì)算每個(gè)文庫(kù)數(shù)據(jù)的覆蓋深度。雖然每個(gè)文庫(kù)均產(chǎn)出了大量的測(cè)序數(shù)據(jù),但由于病毒拷貝數(shù)較低,在Ct<28.7的6個(gè)樣本中,僅有0.004 5%~0.055 8%的測(cè)序數(shù)據(jù)被歸類為SARSCoV-2病毒序列(比例最低的樣本C01和最高的樣本C05除外),而在Ct>28.7的11個(gè)樣本中,僅有0.001 4%~0.006 2%的測(cè)序數(shù)據(jù)被歸類為SARSCoV-2。在19個(gè)測(cè)序結(jié)果中,只有10個(gè)樣本的有效參考基因組覆蓋率超過95%,且測(cè)序覆蓋深度大于5 X(表2)。
然后,我們構(gòu)建了14個(gè)基于雜交捕獲的文庫(kù)并進(jìn)行測(cè)序以獲得病毒的基因組數(shù)據(jù),其中包括9個(gè)宏轉(zhuǎn)錄測(cè)序數(shù)據(jù)基因組覆蓋度不足的樣本(F05、F10、F01、C04、F02、F06、C03、F09和F11)和5個(gè)宏轉(zhuǎn)錄測(cè)序數(shù)據(jù)基因組高覆蓋度的樣本(C06、F12、C05、F04和F07,包含不同的Ct值)。雜交捕獲測(cè)序的每個(gè)樣本平均得到了約7 000萬(wàn)條病毒測(cè)序序列條數(shù),分類鑒定得到SARS-CoV-2病毒測(cè)序數(shù)據(jù)的平均比率為8.4%(表2),其比例顯著高于宏轉(zhuǎn)錄組測(cè)序所得的結(jié)果,表明捕獲測(cè)序法更適合于高Ct(低病毒載量)樣本的測(cè)序。該方法測(cè)序的所有樣品覆蓋深度均大于5 X,基因組覆蓋率為95%(表2),可滿足用于突變和組裝分析。
為了比較兩種不同建庫(kù)方法大規(guī)模并行MPS結(jié)果的準(zhǔn)確性,我們使用宏轉(zhuǎn)錄組和雜交捕獲方法分別對(duì)5個(gè)樣本(C06、F12、C05、F04和F07)進(jìn)行了文庫(kù)構(gòu)建及測(cè)序。分析表明,每個(gè)文庫(kù)測(cè)序數(shù)據(jù)的參考基因組深度都超過了5 X,且參考基因組覆蓋率均大于96%(圖1a、表2)。使用Pilon和FreeBayes方法對(duì)每個(gè)測(cè)序數(shù)據(jù)進(jìn)行SNVs檢測(cè),取兩個(gè)方法都認(rèn)為是可靠的SNVs作為該測(cè)序數(shù)據(jù)檢測(cè)到的SNVs。結(jié)果表明,不論在高病毒拷貝數(shù)樣本(C06、F12)中,還是在低病毒拷貝數(shù)樣本(C05、F04)和極低病毒拷貝數(shù)樣本(F07)中,兩個(gè)不同建庫(kù)方法測(cè)序數(shù)據(jù)檢測(cè)到的SNVs完全一樣,具有很高的一致性(表3)。上述結(jié)果顯示,相對(duì)直接進(jìn)行宏轉(zhuǎn)錄組測(cè)序,雜交捕獲測(cè)序?qū)Σ煌珻t值的樣本均具有準(zhǔn)確的突變檢測(cè)結(jié)果,且所需樣本量?jī)H為宏轉(zhuǎn)錄組測(cè)序的1/10或更少。高Ct(低病毒載量)的樣品即使測(cè)序超過100 Gb也不能獲得有效的SARS-CoV-2病毒測(cè)序信息,但捕獲測(cè)序可獲得有效的信息,表明雜交捕獲測(cè)序是一種更適合于低病毒載量樣本的檢測(cè)方法。
表2 19個(gè)新冠病例樣本宏轉(zhuǎn)錄組測(cè)序數(shù)據(jù)和雜交捕獲測(cè)序數(shù)據(jù)與參考基因組比對(duì)統(tǒng)計(jì)Tab. 2 Sequencing data and SARS-CoV-2 reference genome mapping statistics
表3 宏轉(zhuǎn)錄組測(cè)序和雜交捕獲測(cè)序檢測(cè)的SNVs變異信息Tab. 3 SNVs detection by metatranscriptomic and hybrid capture
合并兩個(gè)不同建庫(kù)方式的測(cè)序數(shù)據(jù)后,總共有10個(gè)樣品的測(cè)序深度大于15 X,其余9個(gè)樣品測(cè)序深度為5~15 X,全部樣品的參考基因組覆蓋率均大于93%(圖1b)。使用SPAdes或IDBA_hyrid進(jìn)行組裝,共得到10個(gè)完整的SARS-CoV-2病毒基因組序列和9個(gè)SARS-CoV-2草圖基因組序列,其中9個(gè)草圖基因組序列含有12~2 174 bp的gaps(表4)。除F09的覆蓋率為93%外,所有組裝的基因組相對(duì)于參考基因組的覆蓋率均在95%以上。
圖1 SARS-CoV-2參考基因組(EPI_ISL_402119)覆蓋深度情況Fig. 1 The sequencing coverage depth of SARS-CoV-2 reference genome (EPI_ISL_402119)(a)5個(gè)不同Ct樣本宏轉(zhuǎn)錄組測(cè)序和雜交捕獲測(cè)序參考基因組覆蓋深度;(b)19個(gè)樣品全部測(cè)序數(shù)據(jù)統(tǒng)計(jì)參考基因組覆蓋深度。(a) The sequencing depth of 5 different Ct samples by metatranscriptomic sequencing and hybrid capture sequencing; (b) Statistical sequencing depth for all sequencing data of 19 samples.
使用FreeBayes、Pysamstats和Pilon 3種方法對(duì)19個(gè)樣本進(jìn)行SNVs檢測(cè),取至少2種方法檢測(cè)到的高置信SNVs為最終鑒定的SNVs (表5)。在19個(gè)樣品中總共檢測(cè)到了ORF1a:C409S、ORF1a:P1207H、ORF8:L84S和S:D614G 4種突變型,分布在基因ORF1a、ORF8和S上(表4、5)。突變分支S(Clade S,ORF8:L84S)包括4個(gè)中國(guó)樣本(C01、C02、C07和C03)和唯一的西班牙樣本F01,這些樣品都是在基因ORF8的第251位堿基由T變成C,導(dǎo)致病毒第84位編碼氨基酸由亮氨酸(Leu,L)突變?yōu)榻z氨酸(Ser,S)。突變分支G(Clade G,S:D614G)包括來(lái)自同一家意大利餐廳的4個(gè)樣本(F07、F08、F09和F10)以及另外5個(gè)意大利樣本(F05、F02、F03、F04和F12)。這9個(gè)意大利樣本的S蛋白具有相同的突變,都是SRAS-CoV-2病毒第614位氨基酸由天冬氨酸(Asp,D)變成了甘氨酸(Gly,G)。此外,具有相同突變ORF1a:P1207H的樣本包括2個(gè)來(lái)自中國(guó)慶元市的樣本(C05和C06)和2個(gè)來(lái)自同一家意大利餐廳的樣本(F06和F12),具有突變ORF1a:C409S的只有1個(gè)來(lái)自中國(guó)的樣本C04。變異分析表明,在中國(guó)及意大利的新型冠狀病毒肺炎患者檢測(cè)到的SARS-CoV-2病毒基因組均存在不一樣的變異類群。同時(shí)樣本間檢測(cè)具有差異的SNVs,表明了SARS-CoV-2病毒變異的多態(tài)性。
表4 病毒基因組組裝和突變類型統(tǒng)計(jì)Tab. 4 Genome assembly and mutation clade statistics
表5 (續(xù))Tab. 5 (continued)
表5 19個(gè)新冠病例樣本SNVs變異檢測(cè)結(jié)果Tab. 5 SNVs detection of 19 sequencing SARS-CoV-2 samples
表5 (續(xù))Tab. 5 (continued)
我們對(duì)高測(cè)序深度(>40 X)的9個(gè)樣本進(jìn)行了iSNVs分析。結(jié)果表明:在F03、F08和C06 3個(gè)樣本中沒有檢測(cè)到任何iSNVs;在C02、C07、C05和F04這4個(gè)樣本中至少檢測(cè)到一個(gè)iSNVs,但沒有iSNVs位于基因區(qū)或非翻譯區(qū)(untranslated region,UTR)。此外,在F06和F11這2個(gè)樣本中均檢測(cè)到了4個(gè)iSNVs,都具有相同的突變c.1841A>G(p.Asp614Gly,即S:D614G)(表6)。F06和F11樣本基因區(qū)的iSNVs與來(lái)自同一家意大利餐廳的其他5個(gè)病例樣本的SNVs突變結(jié)果意外一致,表明在同一患者中可能存在重復(fù)感染。
表6 高深度測(cè)序樣品iSNVs變異信息統(tǒng)計(jì)Tab. 6 The statistics of iSNVs for high depth sequencing samples
本研究從新冠肺炎病例咽拭子和痰液樣本中提取Total RNA,并采用宏轉(zhuǎn)錄組和雜交捕獲方法進(jìn)行高通量全基因組測(cè)序,最終獲得了19個(gè)樣本有效的高通量測(cè)序數(shù)據(jù),組裝得到了10個(gè)完整的新冠病毒基因組序列和9個(gè)新冠病毒基因組草圖。通過比較宏轉(zhuǎn)錄組測(cè)序和捕獲測(cè)序的參考基因組覆蓋度以及變異檢測(cè)情況,證實(shí)了捕獲測(cè)序?qū)τ诓煌《据d量的樣品均可獲得有效、準(zhǔn)確的組裝和變異檢測(cè)結(jié)果。
此外,通過變異結(jié)果分析,11個(gè)意大利歸國(guó)華人中9個(gè)病例均檢測(cè)到了S蛋白D614G突變,且具有D614G突變的樣本大部分都伴有5'UTR中的C到T突變(相對(duì)于參考基因組EPI_ISL_402119的241位)、3 037位的C到T突變、14 408位的C到T突變。另外2個(gè)病例雖然直接檢測(cè)到的變異型不屬于D614G,但在iSNVs檢測(cè)中也均發(fā)現(xiàn)了D614G突變。具有D614G突變的SARS-CoV-2病毒株2020年2月初開始在歐洲傳播[23]。我們3月初采集的來(lái)自意大利不同地區(qū)的F02、F03、F04、F07、F08、F10和F12 7個(gè)樣品都能檢測(cè)到同時(shí)含有這4種突變的相同突變型(表5),表明在3月初之前具有這4個(gè)遺傳連鎖突變的單倍型已在意大利不同地區(qū)開始傳播流行。根據(jù)目前研究,D614G突變毒株已成為世界范圍內(nèi)的主流毒株,在歐洲和北美的測(cè)序樣本中占近70%,且其具有更強(qiáng)的感染性和傳播性[24],這可能是意大利疫情暴發(fā)的一個(gè)重要因素。
我們對(duì)19個(gè)樣本進(jìn)行SNVs檢測(cè),結(jié)果共發(fā)現(xiàn)4種突變,其中意大利患者中分離的病毒株是具有高傳染性的S:D614G突變病毒株,而ORF1a:C409S、ORF1a:P1207H和ORF8:L84S這3個(gè)突變主要在中國(guó)人群中被檢測(cè)到,可見SARS-CoV-2病毒突變?cè)谝獯罄椭袊?guó)并不相同。鑒于SARS-CoV-2基因組會(huì)不斷的發(fā)生突變,那么對(duì)SARS-CoV-2基因組的測(cè)序就顯得極為重要,而我們新研發(fā)的捕獲測(cè)序的方法與傳統(tǒng)的宏轉(zhuǎn)錄組測(cè)序相比所需的樣本量更少,是一種更適合于低病毒載量樣本的檢測(cè)方法。通過捕獲測(cè)序獲得SARS-CoV-2病毒的突變信息并研究突變的功能,對(duì)于后續(xù)抗體和疫苗設(shè)計(jì)具有重要的指導(dǎo)意義。