李文杰，孫之榮

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

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

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

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

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

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

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

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

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

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

2 結(jié)果

2.1 VarScan2 工具分析

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

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

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

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

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

表3 VarScan發(fā)現(xiàn)的插入/刪除突變位點在樣品中的基因型統(tǒng)計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個涉及插入/刪除的突變位點。對于涉及插入/刪除的突變，此工具在相鄰的位點輸出相似的序列，所以導(dǎo)致過多的榮譽突變位點。雖然此工具支持一個位點有多種等位序列，但位點的基因型只支持兩種等位序列。

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

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

2.3 MuTect工具分析

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

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

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

3 討論及結(jié)論

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

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

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

當(dāng)前，ANNOVAR是注釋突變功能影響的重要工具之一［8］，這個工具(及研究者使用的其它工具)都會忽略不會引起蛋白質(zhì)序列變化的突變。但是改變蛋白質(zhì)序列并不是基因組突變產(chǎn)生功能影響的唯一途徑:改變序列對轉(zhuǎn)錄因子或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.