黨 孝，孫宇輝，蔣廷亞，周 陽，連超群

(1.美國費(fèi)城兒童醫(yī)院，賓夕法尼亞州 費(fèi)城 19146；2.蘇州奧根診斷，江蘇 蘇州 215000；3.江蘇大學(xué) 生命科學(xué)研究院，江蘇 鎮(zhèn)江 212013；4.蚌埠醫(yī)學(xué)院 醫(yī)學(xué)臨床檢驗(yàn)診斷中心，安徽 蚌埠 230036)

新一代高通量測序技術(shù)的出現(xiàn)及其快速發(fā)展促進(jìn)了遺傳病致病基因的發(fā)現(xiàn)與臨床分子診斷。在基因檢測中，對患者基因組測序產(chǎn)生的大量候選突變位點(diǎn)中發(fā)現(xiàn)致病突變是科研以及臨床分子診斷需要解決的關(guān)鍵科學(xué)問題。準(zhǔn)確區(qū)分中性突變與有害突變對遺傳病的臨床檢測有著重要的意義，目前的研究方法主要是參考突變的功能影響，群體突變頻率信息，基因功能信息以及家系或多樣本的位點(diǎn)驗(yàn)證信息[1]。研究表明，對于單個(gè)樣本的外顯子測序數(shù)據(jù)，即使過濾了群體常見等位基因頻率(minor allele frequency(MAF)>1%)與位點(diǎn)功能，最終仍然有大約400個(gè)左右的罕見非同義候選突變位點(diǎn)[2，3]。

目前國內(nèi)外對于突變的危害性預(yù)測已經(jīng)開發(fā)出了多個(gè)不同的軟件，從預(yù)測原理及預(yù)測方法上區(qū)分，其主要基于：1)蛋白質(zhì)功能的改變：主要是突變引起蛋白質(zhì)功能發(fā)生變化，如PolyPhen-2、SIFT、MutationTaster、FATHMM和 MutationAssessor等；2)進(jìn)化保守性：主要是對多個(gè)物種核酸序列或蛋白序列進(jìn)行多序列比對，分析同源序列的多態(tài)性，如GERP++、SiPhy、fitCons和PhyloP等；3)整合型軟件：主要是利用機(jī)器學(xué)習(xí)等算法整合多個(gè)其他單獨(dú)的預(yù)測軟件為一整體分析預(yù)測，如CADD、fathmm-MKL、MetaLR、MetaSVM、VEST3、Eigen、DANN、GenoCanyon、REVEL和M-CAP等。

有研究表明整合型軟件的預(yù)測效果整體優(yōu)于單個(gè)軟件[4]，因此此類軟件已廣泛應(yīng)用于位點(diǎn)的危害性預(yù)測。對于這些整合型突變危害性預(yù)測軟件，由于其建立在不同的理念與算法基礎(chǔ)上，基于不同的訓(xùn)練集，因此預(yù)測結(jié)果的準(zhǔn)確性與特異性一直是需要評估的重點(diǎn)；評估結(jié)果也能更好地指導(dǎo)研究人員選擇軟件進(jìn)行位點(diǎn)危害性預(yù)測。本研究收集兩份獨(dú)立的測試數(shù)據(jù)集和一份真實(shí)樣本測序數(shù)據(jù)集對整合型預(yù)測軟件(CADD v1.3[5],fathmm-MKL[6],MetaLR,MetaSVM[4],VEST3 v3.0[7],Eigen v1.1[8],DANN[9],GenoCanyon v1.0.3[10],REVEL[11]和M-CAP v1.0[12])進(jìn)行評估，由于REVEL和M-CAP主要針對罕見錯(cuò)義突變的預(yù)測，為了公平評估，本文主要評測不同危害性預(yù)測軟件對罕見錯(cuò)義突變的預(yù)測；在收集的致病突變位點(diǎn)中也發(fā)現(xiàn)錯(cuò)義突變占了致病SNP位點(diǎn)的絕大部分。

1 材料與方法

1.1 致病位點(diǎn)和中性位點(diǎn)測試數(shù)據(jù)集

收集兩份獨(dú)立測試數(shù)據(jù)集用于評測整合型突變危害性預(yù)測軟件的效果，測試集1：ClinVar數(shù)據(jù)集，致病和中性的突變位點(diǎn)來源于ClinVar數(shù)據(jù)庫[13，14](variant_summary.txt.gz文件，下載于ClinVar FTP，數(shù)據(jù)更新時(shí)間：09/10/2018)，選取基因組為GRCh37的種系(germline)SNP位點(diǎn)，其中致病的SNP選取已報(bào)道臨床意義為“Pathogenic”或“Likely pathogenic”的位點(diǎn)，同時(shí)為了保證致病位點(diǎn)的可靠性，過濾掉審查狀態(tài)為沒有明確證據(jù)顯示是否致病的位點(diǎn)(no assertion criteria provided)；中性的SNP選取報(bào)道臨床意義為“Benign”或“Likely benign”的位點(diǎn)。測試集2：UniProt數(shù)據(jù)集，致病和中性的突變位點(diǎn)來源于UniProt/Swiss-Prot[15](數(shù)據(jù)更新時(shí)間10/10/2018)，其中致病的SNP選取報(bào)道標(biāo)簽為“Disease”的位點(diǎn)，中性的SNP選取標(biāo)簽為“Polymorphism”的位點(diǎn)。由于UniProt數(shù)據(jù)庫上記錄的位點(diǎn)是氨基酸突變形式，因此使用TransVar[16]將氨基酸突變形式轉(zhuǎn)化為GRCh37基因組坐標(biāo)；如果測試數(shù)據(jù)集2中有出現(xiàn)數(shù)據(jù)集1的位點(diǎn)，則直接從數(shù)據(jù)集2中過濾掉此位點(diǎn)以保證兩個(gè)數(shù)據(jù)集的獨(dú)立性。

1.2 整合型預(yù)測軟件對罕見錯(cuò)義突變預(yù)測評估

對于上述兩份測試數(shù)據(jù)集，過濾只留下在公共群體數(shù)據(jù)庫(1000 Genomes Project[2]，Exome Sequencing Project(ESP)[17]，Exome Aggregation Consortium(ExAC)[3]和UK10K[18])中次等位基因頻率(MAF)<1%的罕見錯(cuò)義突變，然后分別使用CADD、fathmm-MKL、MetaLR、MetaSVM、VEST3、Eigen、DANN、GenoCanyon、REVEL和M-CAP對這些位點(diǎn)進(jìn)行預(yù)測，預(yù)測效能的評估利用受試者工作特征曲線ROC(Receiver Operating Characteristic)以及ROC曲線下面積AUC進(jìn)行比較，為了保證評估的公平性，選取了在所有要評估的預(yù)測軟件上都有預(yù)測分值的突變位點(diǎn)進(jìn)行ROC評估。ROC曲線的生成以及不同預(yù)測軟件AUC值的計(jì)算使用R語言程序包“ROCR”[19]。

1.3 真實(shí)測序數(shù)據(jù)評估

除了上述兩份獨(dú)立的測試數(shù)據(jù)集，又收集一份真實(shí)的測序數(shù)據(jù)以評估這些軟件在臨床分子診斷上的效果。在EGA(https://www.ebi.ac.uk/ega/home)上申請來自于NIHR BioResource Rare Disease Consortium對遺傳性眼病樣本的測序數(shù)據(jù)(EGA號：EGAD00001002656，數(shù)據(jù)的bam文件來自于全外顯子測序，CRAM文件來自于全基因組測序)，這些樣本中404位病人通過測序以及臨床分子診斷已經(jīng)找到其致病位點(diǎn)[20]。對申請下載的BAM或CRAM數(shù)據(jù)使用GATK v3.6 HaplotypeCaller call變異[21](只選擇外顯子區(qū)域Agilent SureSelect Human All Exon V5的變異位點(diǎn))，對call出的變異過濾“LowQual”以及深度小于5的SNP位點(diǎn)，然后對剩余的候選SNP位點(diǎn)使用M-CAP、fathmm-MKL、CADD、VEST3、REVEL、MetaLR和MetaSVM分別進(jìn)行危害性預(yù)測。在這404份明確致病突變的樣本中，篩選出致病錯(cuò)義突變位點(diǎn)，為保證評估的公平性，挑選所有致病的錯(cuò)義突變在上述要評估的預(yù)測軟件中都有預(yù)測分值，最終得到118個(gè)致病的罕見錯(cuò)義突變，來源于123個(gè)眼科樣本。

2 結(jié)果

2.1 致病位點(diǎn)的頻率以及功能分布

致病突變的收集來源于ClinVar和UniProt數(shù)據(jù)庫，最終在ClinVar中收集了30821個(gè)以及在UniProt中收集了14322個(gè)致病的SNP，致病突變在群體中的發(fā)生率低，因此大多為罕見突變。對這些致病突變在千人基因組1000 Genomes Project，Exome Sequencing Project(ESP)以及Exome Aggregation Consortium(ExAC)中次等位基因頻率的分布也印證了這些致病位點(diǎn)中約99%以上都為罕見突變(MAF<1%)(圖1)；同時(shí)統(tǒng)計(jì)這些致病突變的突變類型，突變類型的注釋基于Variant Effect Predictor(VEP)[22]，發(fā)現(xiàn)錯(cuò)義突變占所有致病SNP的46%以上(表1)。錯(cuò)義突變占了單核苷酸突變類型的大部分，因此難以實(shí)現(xiàn)對疾病的致病性判斷以及變異解讀，根據(jù)美國醫(yī)學(xué)遺傳學(xué)與基因組學(xué)學(xué)會(ACMG)制定的序列變異解讀指南，大多數(shù)錯(cuò)義突變不能明確分為致病或良性，只能判定為意義未明變異(VUS)[23]，這對遺傳病的臨床診斷造成一定的困難，危害性預(yù)測軟件的開發(fā)與實(shí)現(xiàn)可以輔助對位點(diǎn)致病性的判斷，從而幫助研究人員發(fā)現(xiàn)疾病致病位點(diǎn)或者明確臨床分子診斷。

圖1 ClinVar和UniProt數(shù)據(jù)庫中致病SNP頻率分布圖

致病SNP類型ClinVarUniProtSplice acceptor variant21688Splice donor variant27238Stop gained1043121Stop lost435Initiator codon variant31921Missense variant1416014250Splice region variant702-Synonymous variant7095_prime_UTR variant23-3_prime_UTR variant11-non_coding transcript variant26-Intron variant145-

圖2 測試數(shù)據(jù)集ROC曲線圖((a)為ClinVar測試數(shù)據(jù)集結(jié)果；(b)為UniProt測試數(shù)據(jù)集結(jié)果，兩個(gè)測試數(shù)據(jù)集相互獨(dú)立)

2.2 整合型預(yù)測軟件對ClinVar和UniProt數(shù)據(jù)庫罕見錯(cuò)義突變預(yù)測

對收集的ClinVar和UniProt測試數(shù)據(jù)集，過濾群體頻率與突變類型留下罕見錯(cuò)義突變位點(diǎn)，并確保所有位點(diǎn)在要評估的預(yù)測軟件上都有預(yù)測分值。最終測試數(shù)據(jù)集ClinVar剩余12316個(gè)罕見錯(cuò)義致病突變和7988個(gè)罕見中性錯(cuò)義突變；測試數(shù)據(jù)集UniProt剩余11818個(gè)罕見錯(cuò)義致病突變和10427個(gè)罕見中性錯(cuò)義突變(附件1)。將上述兩個(gè)測試集分別進(jìn)行ROC曲線分析，通過ROC曲線以及曲線下面積AUC值的比較，綜合兩個(gè)測試集結(jié)果(圖2(a) ClinVar測試數(shù)據(jù)集結(jié)果；圖2(b) UniProt測試數(shù)據(jù)集結(jié)果)，發(fā)現(xiàn)在罕見錯(cuò)義致病突變與中性突變的識別上，REVEL效果最好，準(zhǔn)確性以及特異性明顯優(yōu)于其他軟件，其他效果較好有M-CAP、MetaSVM、MetaLR和VEST3，其準(zhǔn)確性與特性性高于剩余的5個(gè)軟件。兩個(gè)測試數(shù)據(jù)集是獨(dú)立的，因此評估效果更能代表無偏性，避免評估軟件在特定數(shù)據(jù)集上的過優(yōu)或過差表現(xiàn)。

2.3 來源于123個(gè)眼科樣本的118個(gè)罕見錯(cuò)義致病突變不同軟件的預(yù)測

通過比較M-CAP、fathmm-MKL、CADD、VEST3、REVEL、MetaLR和MetaSVM在真實(shí)123個(gè)眼科測序樣本的基因組罕見SNP危害性預(yù)測效果，同時(shí)重點(diǎn)關(guān)注這些樣本的118個(gè)致病突變的預(yù)測結(jié)果(附件2)。統(tǒng)計(jì)這些軟件對123個(gè)樣本總共預(yù)測的罕見唯一的致病位點(diǎn)數(shù)目，同時(shí)這118個(gè)致病位點(diǎn)不同軟件的預(yù)測結(jié)果分為兩類：“Damaging”和“Tolerated”，比較發(fā)現(xiàn)盡管CADD和FATHMM-MKL對118個(gè)致病位點(diǎn)預(yù)測的Damaging數(shù)量多，但也將更多的中性突變預(yù)測為有害，因此假陽性率相對較高。REVEL、MetaLR和MetaSVM預(yù)測的Damaging數(shù)目少，但同時(shí)對其他中性突變預(yù)測為有害的數(shù)量低，因此在選擇危害性預(yù)測軟件時(shí)需要根據(jù)自己的研究需求，權(quán)衡假陽性率與假陰性率(圖3，表2)。一個(gè)典型的單基因病患者一般攜帶1～2個(gè)致病突變，其余的大部分突變?yōu)橹行酝蛔僛12]。

圖3 不同軟件對118個(gè)罕見錯(cuò)義致病突變的預(yù)測圖(注：左邊坐標(biāo)表示致病突變數(shù)目，柱形圖表示預(yù)測“Damaging”和“Tolerated”的數(shù)量；右坐標(biāo)表示118個(gè)致病位點(diǎn)預(yù)測 Damaging的數(shù)量與123個(gè)樣本中總共預(yù)測的致病的唯一突變數(shù)目的百分比值，由折線圖展示)

軟件有害閾值118個(gè)致病位點(diǎn)預(yù)測Damaging數(shù)目118個(gè)致病位點(diǎn)預(yù)測 Tolerated數(shù)目123個(gè)樣本中總共預(yù)測的致病的唯一突變數(shù)目REVEL>0.589295,304MetaSVM>080385,420MetaLR>0.582366,160M-CAP>0.025115316,636CADD>201071125,394VEST3>0.5952313,695fathmm-MKL>0.5111728,462

3 討論

在已發(fā)現(xiàn)的遺傳病致病位點(diǎn)中很大一部分是罕見錯(cuò)義突變，區(qū)分突變的致病性對疾病致病基因的發(fā)現(xiàn)與分子診斷起著重要作用，本研究通過對基于機(jī)器學(xué)習(xí)開發(fā)的不同危害性預(yù)測軟件使用兩個(gè)獨(dú)立的測試數(shù)據(jù)集進(jìn)行預(yù)測評估，結(jié)果表明，REVEL效果最好，準(zhǔn)確性以及特異性優(yōu)于其他軟件，其他效果較好的軟件有M-CAP、MetaSVM、MetaLR和VEST3，盡管CADD、fathmm-MKL和DANN在罕見錯(cuò)義突變的預(yù)測上效果不突出，但他們對于非編碼以及其他類型突變的危害性預(yù)測起著重要的作用，REVEL和M-CAP是最新開發(fā)的專門針對罕見錯(cuò)義突變的軟件，因此在使用不同軟件時(shí)需結(jié)合具體需求選擇針對性的軟件，輔助位點(diǎn)危害性的判斷。本文也評估了部分軟件在已明確致病突變的真實(shí)測序樣本上的表現(xiàn)，結(jié)果提示在使用危害性預(yù)測軟件時(shí)需要權(quán)衡假陽性率與假陰性率，而且預(yù)測結(jié)果只能作為對位點(diǎn)致病性的輔助判斷，不能作為篩選或確定信息予以使用。

危害性預(yù)測軟件的效果評估一直以來受到特定數(shù)據(jù)集以及突變類型的影響，研究表明不同的預(yù)測軟件預(yù)測結(jié)果不一致率高，在不同的測試數(shù)據(jù)集上準(zhǔn)確性與特異性波動大。不同的預(yù)測軟件建立在不同的理念與算法基礎(chǔ)上，基于不同的訓(xùn)練集，即使針對同一個(gè)位點(diǎn)，不同的軟件對突變的危害性預(yù)測結(jié)果很不一致[24]。而且相同的軟件對于不同的位點(diǎn)，不同的測試數(shù)據(jù)集預(yù)測結(jié)果的準(zhǔn)確性與特異性波動大[25]。多個(gè)預(yù)測軟件不同的預(yù)測結(jié)果會影響對突變危害性的判斷，為了提升危害性預(yù)測軟件的準(zhǔn)確性，目前以更為具體的疾病，基因或通路信息研究是危害性預(yù)測軟件提升的一個(gè)方向，如基于部分有突變熱點(diǎn)或重要功能的基因的危害性預(yù)測軟件開發(fā)[26]；對于特定基因家族的危害性預(yù)測軟件開發(fā)[27]；基于不同的生物通路信息開發(fā)的危害性預(yù)測軟件等。

4 致謝

感謝NIHRBioResource， University of Cambridge和NIHR BioResource Rare Diseases BRIDGE consortium允許下載使用他們測序的遺傳性眼病樣本數(shù)據(jù)。感謝費(fèi)城兒童醫(yī)院Center for Data Driven Discovery of Biomedicine提供訪問交流機(jī)會。