基于第二代測(cè)序數(shù)據(jù)識(shí)別腫瘤基因突變的工具比較

李文杰，孫之榮

(清華大學(xué)生命科學(xué)學(xué)院，教育部生物信息學(xué)重點(diǎn)實(shí)驗(yàn)室，北京100084)

癌癥(即惡性腫瘤)一直是導(dǎo)致人類(lèi)死亡的重要原因之一，也一直是世界范圍內(nèi)的科學(xué)家正在攻克的焦點(diǎn)，科學(xué)家們一直在探尋其發(fā)病原因和有效治療方法。雖然癌癥的病因一般可歸納為原癌基因和抑癌基因中的基因組突變，但具體的突變因病人不同和發(fā)生組織不同而不同，所以發(fā)現(xiàn)某個(gè)病人的致癌突變很有挑戰(zhàn)性。

腫瘤發(fā)生過(guò)程中出現(xiàn)的幾類(lèi)基因組突變包括:點(diǎn)突變(如單堿基替換、堿基插入、堿基刪除)，拷貝數(shù)變化(Copy number variation，CNV)，染色體結(jié)構(gòu)變化(如基因融合)。第二代測(cè)序技術(shù)(Next Genera-tion Sequencing，NGS)作為一種新的實(shí)驗(yàn)技術(shù)，給癌癥研究帶來(lái)了很多重要的發(fā)現(xiàn)。在比較腫瘤和正常樣品的測(cè)序數(shù)據(jù)過(guò)程中，最重要的是高正確率和高特異性地找到點(diǎn)突變。

因其低成本性和高通量性，第二代測(cè)序技術(shù)正被研究者使用，以期找到癌癥的起因。每次測(cè)序?qū)嶒?yàn)產(chǎn)生數(shù)以百萬(wàn)計(jì)的短read，每個(gè)位點(diǎn)的堿基都有測(cè)序質(zhì)量。通過(guò)比較一個(gè)病人的正常和腫瘤組織樣品的基因組測(cè)序數(shù)據(jù)，我們可較容易地確定腫瘤組織中發(fā)生了哪些點(diǎn)突變。但是樣品中細(xì)胞純度、測(cè)序錯(cuò)誤、堿基測(cè)序質(zhì)量、read比對(duì)(與參考基因組比對(duì))，都給這個(gè)任務(wù)帶來(lái)一定的挑戰(zhàn)，特別是對(duì)設(shè)計(jì)插入和刪除的點(diǎn)突變。

許多現(xiàn)有的工具先過(guò)濾掉一些read和堿基，然后計(jì)算報(bào)導(dǎo)位點(diǎn)各allele序列的read數(shù)目，之后比較位點(diǎn)在正常和腫瘤樣品中各allele的read頻率。但是這些工具的驗(yàn)證正確率一般都只在54%左右［1］，不同工具間的一致性都比較低［2］。

最多被使用的工具有 VarScan［3］，SAMtools［4］和MuTect［5］。此文將這些工具應(yīng)用在一個(gè)肺癌病人的血液和癌轉(zhuǎn)移組織的基因組測(cè)序數(shù)據(jù)上，分析它們結(jié)果的合理性。同時(shí)，此文還獲取了該病例的癌轉(zhuǎn)移組織的轉(zhuǎn)錄組測(cè)序數(shù)據(jù)，并據(jù)此通過(guò)它來(lái)驗(yàn)證各個(gè)工具發(fā)現(xiàn)的各個(gè)突變位點(diǎn)，以獲得工具的驗(yàn)證正確率。

1 數(shù)據(jù)和方法

此文選擇一個(gè)肺癌病人的血液組織作為正常組織，癌轉(zhuǎn)移作為腫瘤組織(從肺轉(zhuǎn)移至肝，肺部的原癌組織數(shù)據(jù)質(zhì)量很差)［6］。此文應(yīng)用Bowtie2工具，將測(cè)序產(chǎn)生的paired-end read比對(duì)到人類(lèi)參考基因組上(hg19)，然后利用SAMtools工具的部分命令，生成按比對(duì)位置排序的BAM文件［7］。

對(duì)于各個(gè)工具，此文簡(jiǎn)要描述它使用的方法，然后在分析上述樣品測(cè)序數(shù)據(jù)時(shí)，使用默認(rèn)的各閾值和推薦的步驟。各工具比較兩個(gè)樣品的全基因組測(cè)序數(shù)據(jù)，以確定發(fā)生了點(diǎn)突變的基因組位點(diǎn)。

本文對(duì)各突變位點(diǎn)判斷腫瘤特異的等位序列，并用病人癌轉(zhuǎn)移組織的轉(zhuǎn)錄組測(cè)序(RNA-Seq)數(shù)據(jù)判斷這個(gè)序列的真實(shí)存在性。此方法只驗(yàn)證了被轉(zhuǎn)錄了的腫瘤特異等位序列，但也可以作為工具正確率的一種衡量。

2 結(jié)果

2.1 VarScan2 工具分析

VarScan2使用各種過(guò)濾方法，以得到各位點(diǎn)在樣品中的基因型，并計(jì)算各等位序列的read頻率。此工具只支持最多具有兩種等位序列的基因型，一個(gè)與參考基因組相同，一個(gè)為變異序列;同時(shí)這個(gè)變異序列在兩個(gè)樣品中需要相同。然后此工具比較兩個(gè)樣品在該位點(diǎn)各等位序列的read頻率，并通過(guò)Fisher’s Exact檢驗(yàn)得到差異顯著性水平 p值。VarScan2工具將突變位點(diǎn)分為三種突變狀態(tài):Somatic突變、Germline突變、LOH突變;并從三種狀態(tài)的位點(diǎn)，得到’High-Confidence’位點(diǎn)，各種狀態(tài)的High-Confidence位點(diǎn)數(shù)目見(jiàn)表1。

表1 VarScan發(fā)現(xiàn)的各狀態(tài)突變位點(diǎn)數(shù)目Table 1 Number of mutation sites identified by VarScan

但是此工具識(shí)別出過(guò)多的單堿基替換突變位點(diǎn)，結(jié)果中存在比較高的假陽(yáng)性率，也讓工具使用者比較難判斷真正的突變位點(diǎn)。對(duì)于Somatic狀態(tài)的突變位點(diǎn)，它在正常樣品中的基因型需要是純合的跟參考基因組一樣(R/R，見(jiàn)表2和表3)。對(duì)于14 087個(gè)位點(diǎn)，在兩個(gè)樣品中有相同的雜合基因型，但是allele的read頻率有差異，這些位點(diǎn)并認(rèn)為是‘LOH’突變，而不是‘Somatic’突變狀態(tài)(見(jiàn)表2)。

此工具發(fā)現(xiàn)的單堿基替換突變位點(diǎn)中，22 207個(gè)位點(diǎn)有腫瘤特異的等位序列，4 108個(gè)被RNA-Seq數(shù)據(jù)覆蓋，1 359個(gè)位點(diǎn)的腫瘤特異等位序列在RNA-Seq數(shù)據(jù)中存在(驗(yàn)證正確率33%);此工具發(fā)現(xiàn)的涉及插入/刪除突變位點(diǎn)中，2 710個(gè)位點(diǎn)有腫瘤特異的等位序列，422個(gè)被RNA-Seq數(shù)據(jù)覆蓋，159個(gè)位點(diǎn)的腫瘤特異等位序列在RNA-Seq數(shù)據(jù)中存在，驗(yàn)證正確率38%。

表2 VarScan發(fā)現(xiàn)的單堿基替換突變位點(diǎn)在樣品中的基因型統(tǒng)計(jì)Table 2 The genotype of SNV sites in both samples

表3 VarScan發(fā)現(xiàn)的插入/刪除突變位點(diǎn)在樣品中的基因型統(tǒng)計(jì)Table 3 The genotype of indel sites in both samples

2.2 SAMtools工具分析

SAMtools能利用正常和腫瘤樣品的BAM文件，來(lái)確定點(diǎn)突變位點(diǎn)。此工具對(duì)每個(gè)突變位點(diǎn)計(jì)算一個(gè)CLR值(在輸出的VCF文件中)，以表明位點(diǎn)在兩個(gè)樣品間的差異顯著性大小，值越大表明越顯著，范圍在0～255(見(jiàn)圖1)。此文使用70為閾值，只保留差異顯著性大的突變位點(diǎn)，1 012個(gè)單堿基替換突變位點(diǎn)，2 578個(gè)涉及插入/刪除的突變位點(diǎn)。對(duì)于涉及插入/刪除的突變，此工具在相鄰的位點(diǎn)輸出相似的序列，所以導(dǎo)致過(guò)多的榮譽(yù)突變位點(diǎn)。雖然此工具支持一個(gè)位點(diǎn)有多種等位序列，但位點(diǎn)的基因型只支持兩種等位序列。

此工具發(fā)現(xiàn)的位點(diǎn)中，627個(gè)單堿基替換突變位點(diǎn)有腫瘤特異等位序列，95個(gè)被RNA-Seq數(shù)據(jù)覆蓋，43個(gè)位點(diǎn)的腫瘤特異等位序列在RNA-Seq數(shù)據(jù)中存在，驗(yàn)證正確率45%;1 575個(gè)涉及插入/刪除的突變位點(diǎn)有腫瘤特異等位序列，186個(gè)被RNASeq數(shù)據(jù)覆蓋，105個(gè)位點(diǎn)的特異等位序列在RNASeq數(shù)據(jù)中存在，驗(yàn)證正確率56%。

圖1 SAMtools計(jì)算的CLR值在突變位點(diǎn)中的分布Fig.1 Distribution of CLR value among mutation sites calculated by SAMtools

2.3 MuTect工具分析

MuTect不支持涉及插入/刪除的突變，它過(guò)濾掉比對(duì)質(zhì)量低的read，過(guò)濾掉覆蓋read數(shù)低的位點(diǎn)，最后過(guò)濾掉有read比對(duì)時(shí)鏈偏好性的位點(diǎn)。此工具對(duì)每個(gè)位點(diǎn)推斷發(fā)生了突變的可能性:腫瘤樣品在此位點(diǎn)與參考基因組序列不同，正常樣品沒(méi)有變異等位序列。對(duì)于每個(gè)突變位點(diǎn)，此工具還計(jì)算其為真實(shí)突變的可能性:變異等位序列是真實(shí)的，而不是測(cè)序錯(cuò)誤，然后經(jīng)過(guò)log10轉(zhuǎn)換得到的值。大多數(shù)的突變位點(diǎn)的這個(gè)可能性值不太高(見(jiàn)圖2)。本文以15為閾值，得到輸出6 679個(gè)單堿基替換突變位點(diǎn)。

圖2 MuTect計(jì)算的突變真實(shí)可能性在突變位點(diǎn)中的分布Fig.2 Distribution of t_log_fstar among mutation sites calculated by MuTect

此工具發(fā)現(xiàn)的單堿基替換突變位點(diǎn)中，5 397個(gè)位點(diǎn)有腫瘤特異等位序列，851個(gè)被RNA-Seq數(shù)據(jù)覆蓋，305個(gè)位點(diǎn)的腫瘤特異等位序列在RNA-Seq數(shù)據(jù)中存在，驗(yàn)證正確率36%。

3 討論及結(jié)論

現(xiàn)有的工具假定read中測(cè)序質(zhì)量較低的堿基是測(cè)序錯(cuò)誤，然后將它們過(guò)濾掉。這樣的方法可能使得結(jié)果對(duì)選擇的閾值敏感。同時(shí)，上述工具只在一個(gè)位點(diǎn)只支持最多兩種等位序列，這可能丟失了腫瘤樣品中出現(xiàn)的很多突變位點(diǎn)，原因在于腫瘤細(xì)胞中，基因組一個(gè)片段可能會(huì)出現(xiàn)多個(gè)拷貝，每個(gè)拷貝可能經(jīng)過(guò)不同的突變過(guò)程，這樣對(duì)應(yīng)到參考基因組上的某個(gè)位點(diǎn)后，可能存在多種等位序列。

各個(gè)工具的結(jié)果一致性很低，很少比例的位點(diǎn)同時(shí)被兩個(gè)工具發(fā)現(xiàn)。這個(gè)低一致性可能來(lái)自于:(1)它們過(guò)濾read和堿基的方式;(2)它們比較兩個(gè)樣品的方法。建議用戶(hù)嘗試多個(gè)方法，并用更精確的測(cè)序方法驗(yàn)證各自的結(jié)果，然后挑選一個(gè)有最高驗(yàn)證正確率的工具來(lái)使用。

除了此文中討論的點(diǎn)突變之外，腫瘤組織中還發(fā)生了其它類(lèi)型的基因組突變，如染色體重組。這些突變可能涉及大范圍內(nèi)的堿基，也更難發(fā)現(xiàn)和驗(yàn)證。完全研究這些基因組突變依然是比較困難的，更何況確定各突變的功能影響。

當(dāng)前，ANNOVAR是注釋突變功能影響的重要工具之一［8］，這個(gè)工具(及研究者使用的其它工具)都會(huì)忽略不會(huì)引起蛋白質(zhì)序列變化的突變。但是改變蛋白質(zhì)序列并不是基因組突變產(chǎn)生功能影響的唯一途徑:改變序列對(duì)轉(zhuǎn)錄因子或miRNA的親和性都會(huì)顯著得影響功能［9］。因此，需要更多的研究來(lái)完全注釋基因組突變的功能影響。

References)

［1］ KENICHI Y，MASASHI S，YUICHI S，et al.Frequent pathway mutations of splicing machinery in myelodysplasia［J］.Nature，2012，478:64－69.

［2］ MARTIN L，BERNHARD Y，JOS DE G，et al.Confidence-based somatic mutation evaluation and pioritization［J］.PLoS Comput Biol，2012，8(9):e1002714.

［3］ DANIEL C，QUNYUAN Z，DAVID E，et al.VarScan 2:Somatic mutation and copy number alteration discovery in cancer by exome sequencing［J］.Genome Research，2012，450:65.

［4］ HENG L，BOB H，ALEC W，et al.The sequence alignment/map(SAM)format and SAMtools［J］.Bioinformatics，2009，25:2078 －9.

［5］ KRISTIAN C，MICHAEL S，SCOTT L，et al.Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples［J］.Nature Biotechnology，2013，410:60.

［6］ YOUNG S，WON-CHUL L，THOMAS B，et al.A transforming KIF5B and REF gene fusion in lung adenocarcinoma revealed from whole-genome and transcriptome sequencing［J］.Genome Research，2012，22(3):436 －445.

［7］ LANGMEAD B，SALZBERG S.Fast gapped-read alignment with Bowtie 2［J］.Nature Methods，2012，9:357 －359.

［8］ KAI W，MINGYAO L，HAKON H，et al.ANNOVAR:functional annotation ofgenetic variantsfrom highthroughput sequencing data［J］.Nucl.Acids Res.，2010，38(16):e164.

［9］ EMANUELA S，CLAUDIA B，Beatrice P，et al.A somatic mutation in the 5＇UTR of BRCA1 gene in sporadic breastcancercauses down-modulation oftranslation efficiency［J］.Oncogene，2012，s:4596 －4600.