基于機器學習開發(fā)的危害性預測軟件對罕見錯義突變的預測評估

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

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

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

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

有研究表明整合型軟件的預測效果整體優(yōu)于單個軟件[4]，因此此類軟件已廣泛應用于位點的危害性預測。對于這些整合型突變危害性預測軟件，由于其建立在不同的理念與算法基礎上，基于不同的訓練集，因此預測結(jié)果的準確性與特異性一直是需要評估的重點；評估結(jié)果也能更好地指導研究人員選擇軟件進行位點危害性預測。本研究收集兩份獨立的測試數(shù)據(jù)集和一份真實樣本測序數(shù)據(jù)集對整合型預測軟件(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])進行評估，由于REVEL和M-CAP主要針對罕見錯義突變的預測，為了公平評估，本文主要評測不同危害性預測軟件對罕見錯義突變的預測；在收集的致病突變位點中也發(fā)現(xiàn)錯義突變占了致病SNP位點的絕大部分。

1 材料與方法

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

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

1.2 整合型預測軟件對罕見錯義突變預測評估

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

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

除了上述兩份獨立的測試數(shù)據(jù)集，又收集一份真實的測序數(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)找到其致病位點[20]。對申請下載的BAM或CRAM數(shù)據(jù)使用GATK v3.6 HaplotypeCaller call變異[21](只選擇外顯子區(qū)域Agilent SureSelect Human All Exon V5的變異位點)，對call出的變異過濾“LowQual”以及深度小于5的SNP位點，然后對剩余的候選SNP位點使用M-CAP、fathmm-MKL、CADD、VEST3、REVEL、MetaLR和MetaSVM分別進行危害性預測。在這404份明確致病突變的樣本中，篩選出致病錯義突變位點，為保證評估的公平性，挑選所有致病的錯義突變在上述要評估的預測軟件中都有預測分值，最終得到118個致病的罕見錯義突變，來源于123個眼科樣本。

2 結(jié)果

2.1 致病位點的頻率以及功能分布

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

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

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

2.3 來源于123個眼科樣本的118個罕見錯義致病突變不同軟件的預測

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

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

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

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

4 致謝

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