基于第二代測序數據識別腫瘤基因突變的工具比較

李文杰，孫之榮

(清華大學生命科學學院，教育部生物信息學重點實驗室，北京100084)

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

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

因其低成本性和高通量性，第二代測序技術正被研究者使用，以期找到癌癥的起因。每次測序實驗產生數以百萬計的短read，每個位點的堿基都有測序質量。通過比較一個病人的正常和腫瘤組織樣品的基因組測序數據，我們可較容易地確定腫瘤組織中發生了哪些點突變。但是樣品中細胞純度、測序錯誤、堿基測序質量、read比對(與參考基因組比對)，都給這個任務帶來一定的挑戰，特別是對設計插入和刪除的點突變。

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

最多被使用的工具有 VarScan［3］，SAMtools［4］和MuTect［5］。此文將這些工具應用在一個肺癌病人的血液和癌轉移組織的基因組測序數據上，分析它們結果的合理性。同時，此文還獲取了該病例的癌轉移組織的轉錄組測序數據，并據此通過它來驗證各個工具發現的各個突變位點，以獲得工具的驗證正確率。

1 數據和方法

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

對于各個工具，此文簡要描述它使用的方法，然后在分析上述樣品測序數據時，使用默認的各閾值和推薦的步驟。各工具比較兩個樣品的全基因組測序數據，以確定發生了點突變的基因組位點。

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

2 結果

2.1 VarScan2 工具分析

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

表1 VarScan發現的各狀態突變位點數目Table 1 Number of mutation sites identified by VarScan

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

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

表2 VarScan發現的單堿基替換突變位點在樣品中的基因型統計Table 2 The genotype of SNV sites in both samples

表3 VarScan發現的插入/刪除突變位點在樣品中的基因型統計Table 3 The genotype of indel sites in both samples

2.2 SAMtools工具分析

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

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

圖1 SAMtools計算的CLR值在突變位點中的分布Fig.1 Distribution of CLR value among mutation sites calculated by SAMtools

2.3 MuTect工具分析

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

圖2 MuTect計算的突變真實可能性在突變位點中的分布Fig.2 Distribution of t_log_fstar among mutation sites calculated by MuTect

此工具發現的單堿基替換突變位點中，5 397個位點有腫瘤特異等位序列，851個被RNA-Seq數據覆蓋，305個位點的腫瘤特異等位序列在RNA-Seq數據中存在，驗證正確率36%。

3 討論及結論

現有的工具假定read中測序質量較低的堿基是測序錯誤，然后將它們過濾掉。這樣的方法可能使得結果對選擇的閾值敏感。同時，上述工具只在一個位點只支持最多兩種等位序列，這可能丟失了腫瘤樣品中出現的很多突變位點，原因在于腫瘤細胞中，基因組一個片段可能會出現多個拷貝，每個拷貝可能經過不同的突變過程，這樣對應到參考基因組上的某個位點后，可能存在多種等位序列。

各個工具的結果一致性很低，很少比例的位點同時被兩個工具發現。這個低一致性可能來自于:(1)它們過濾read和堿基的方式;(2)它們比較兩個樣品的方法。建議用戶嘗試多個方法，并用更精確的測序方法驗證各自的結果，然后挑選一個有最高驗證正確率的工具來使用。

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

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

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.