「生資無價系列專題」- 在缺少正常組織的檢體時，該如何有效分析癌症病人的體細胞變異 (Somatic Mutation)？ （以 Mutect2 為例）
緣起
本文的動念來自於三年前所發生「美麗的錯誤」。
針對遺傳性疾病，利用 NGS 來找尋致病的基因變異時，要特別注意 False Positive 的問題（當然也會有 False Negative，不過，這歸類為「未檢出」的議題，不在本文探討範圍內。）。不然，就會導致後續用 Sanger Sequencing 來驗證時卻發現是假的，導致最後 PI 或醫師又會直接下出這個結論：「NGS 沒用！」
雖然本人在癌症資料分析的經驗尚不足，但仍想藉由自己遇到的實際案例跟大家分享，希望提供大家在分析癌症資料時一些參考和注意事項。也希望藉由此文章跟大家多多交流，更歡迎有經驗的高手多給我回饋，謝謝。
體細胞變異 (Somatic Mutation)
何謂體細胞變異 (Somatic Mutation)? 首先，必須先從生殖細胞變異 ( Germline Mutation) 講起，每個人約有 3～5百萬個 Germline Mutation，這些變異來自個人與 Human Reference Genome [1] 之間的差異，大部分變異來自爸爸媽媽遺傳給我們，少數變異是新生變異 (de Novo mutation) 。 然而，Somatic mutation 是指在癌細胞上產生的變異，如下圖一的最左邊的變異點，在 Normal tissue 上沒看到，在 Tumor tissue 上就發現有部分變異產生，這類變異才會稱為 Somatic Mutation。
當分析癌症的體細胞變異 (Somatic Mutation) 時，其困難度與上述的遺傳性疾病的 Germline Mutation 相比，就像是籃球場上的投 Logo shot 和罰球的比較一樣，難度是數十倍以上（除非你是勇士的 Curry 或拓荒者的 Lillard) ，原因有二：
- 腫癌異質性 (Tumor heterogeneity) [2] : 簡單的說就是「一團亂」，也且是隨著時間越來越亂。這會造成治療時的難度大增，一個簡單的比喻是癌症治療如同三峽大霸同時有 10 個地方潰堤一樣，只補住 9 個缺口是沒用的，因為只要漏一個缺口，就會造成下游地區淹水 (癌細胞擴散或轉移）。如何有效找出所有變異點（潰堤點）是所有研究者努力的方向，而且目前很多標靶用藥得靠 NGS 找出這些變異標靶 (Biomarker) 才能對症下藥。例如非小細胞肺癌病患若在 EGFR 的 exon 18–21 有發生突變，則用艾瑞莎 (Iressa or Gefitinib) 的反應較好；然而若 EGFR 上有 T790M 突變發生時，則艾瑞莎有效機率就會變得極低 [3]。
- 無法用 Sanger 驗證：也是因為 Tumor heterogeneity 的特性，無法如生殖細胞變異 (Germline Mutation) 可以用 Sanger 來進行驗證，導致確認這些體細胞變異的真實性及完整性更加困難。
基於以上原因，如何有效找出所有的 Somatic Mutation 是一個重要的題議，而生物資訊界也有很多團隊投入許多心力在開發更好的偵測工具，今年初在 Scientific Reports 上有篇論文 [4] 對目前市面上大家常用的工具進行評比，發現 Stelka2 和 Mutect2 穎脫而出 (參考自教育部成語典 [5] )，Stelka2 在 Higher mutation frequency (≥20%) 比 Mutect2 好一點，而 Mutect2 在 Lower mutation frequency (≤10%) 比 Stelka2 好一點。雖然 Strelka2 比 Mutect2 快近 20 倍，但利用本公司開發的分散式平行加速系統，只需幾小時就可以得到結果，所以執行時間不會是問題。另外，我們對 Lower mutation frequency 比較有興趣，所以選擇採用 Mutect2 進行 Somatic calling。
Somatic Mutation Caller — Mutect2
Somatic Mutation Caller 有二個重點：第一、把所有可能的 variants 盡可以偵測出來；第二、有效去除 Germline mutation 和系統性的假訊號。
在撰寫此文時，GATK 的 Somatic Mutation Pipeline 仍持續更新中，本文將以 v4.1.0.0 為例來進行說明，雖然它己經被歸類為 (Deprecated)，但並不影響本文想傳達的觀點。
首先，讓我們來看看 Mutect2 有提供那些實用的功能，圖二是擷取自 GATK 官網上的一個範例 [6]，可以清楚看到 Mutect2 所提供的呼叫界面，除了提供 Tumor 和 Normal sample 的 BAM files 外，另外還有提供 PoN 和 Germline Resource 的功能。但何謂 PoN 及 Germline Resource 呢？應該在什麼情況下使用它們呢？應如何準備合適的 PoN 和 Germline Resource 呢？接下來會一一為大家說明。
最有效去除 Germline mutation 的方式就是用 Normal tissue 來進行 Somatic-mode (Tumor/Normal）的分析，但常常在樣本收集或經費限制下而缺少 Normal Tissue，進而無法使用 Somatic mode 來去除 Germline Mutation，在此稱之為 Germline Effect。如何有效避免 Germline Effect，Mutect2 提供以下 2 個方式供大家使用 :
Germline Resource
當缺少了 Normal tissue 的資訊來去除 Germline Effect 時，通常會用 Population database 將常見的 common germline mutation 來做過濾。以目前世界上最大的 population database gnomAD 為例，它在 v2.1 版中包含 ~200k exomes 和 ~16k genomes 的資料。另外，因為全部資料太大，所以有特別整理出只包含 AF 欄位的檔案(af-only-gnomad_grch38.vcf.gz)來提升資料載入時間和減少所需的儲存空間。同時，Mutect2 提供另一個參數 — af-of-alleles-not-in-resource 來指定符合 Germline Mutation 的 AF 值應該多高，以此即可過濾掉同一族群常見的 Common germline variants。gnomAD 是統合了多個不同人種的資訊，若有針對欲進行分析的人種更相近的 cohort data 會更好，前提是人數要夠多。例如，Taiwan BioBank 目前提供近 2000 人的 WGS，即可拿這些資料進行 Germline Joint Genotype Pipeline 來得到屬於台灣人種的 Germline Resource。（注意：建議處理這些資料應從 gVCF 開始，而非單純將每個人的 VCF 直接 merge 起來，避免一些 rare mutation 沒被偵測出來）
PoN
PoN 是 a Panel of Normals 的縮寫，主要用意是在去除系統性的假訊號，包括 mapping artifacts 或 other somewhat random but systematic artifacts of sequencing。GATK 提供的文件建議至少要有 40 normal samples 來建立 PoN才會比較合適 [6]。另外，它提供 — min-sample-count (default: 2) 來指定多少以上才視為屬於 PoN 的假訊號。當然，它也可以用來去除 Common germline variants sites，但我會建議使用上述的 Germline Resource 來實作會比較合適，原因有：
- 變異資訊的精細度上的差異：PoN 主要是提供 Site-only 的過濾而 Germline Resource 是 Allele-specific 的過濾資訊。相對而言，Germline Resource 提供比較細緻的過濾。
- 變異普遍性門檻設定的彈性：PoN 必須在建立時就明確指定，當建立好後，並沒有其它方式來調整不同的 Threshold 做過濾。所以若完全仰賴 PoN 來處理 Germline Effect 和 Systematic Artifact 的話，彈性會被局限住。
基於以上的原因，建議把 PoN 定位成 Systematic Artifact 過濾器，Germline Resource 來過濾 Germline mutation 並可微調 Allele Frequency。
一個美麗的錯誤
話說三年前跟台灣一間專門治療癌症病患的醫院進行臨床研究的合作案時，為了探討為什麼有些人的癌症會復發，同時有些人會在用藥後 5 年內引發另一種癌症，希望可以提供以後的病患一些術後回診頻率及用藥上的建議。
所以我們提了一個研究案去申請了台北市研發補助，在經費有了著落，就可以在院內申請 IRB 進行收案。但因經費被評審委員砍掉一大半，導致最後實際參與研究的病例數大減，也讓後續的分析難度變高以及可信度變低。
在這個研究案中，我們會從病人切除下來的細胞檢體進行冷凍，由病理醫師分別採集出癌細胞和正常細胞的樣本，再去做 WES (24 G raw data)，病患包含前期到晚期，有用藥併發症，也有復發的案例，從臨床資料看來可說是非常豐富及多樣性，相信之後的分析會有很多有趣的發現，但最後總共收集了 16 位病患進行研究，只因為經費的限制。（連 PoN 要求最少的數量一半都不到）
先講結果，最後這個臨床研究案跟醫院端醫師們報告完後，結論就是每個想問的問題我都能回答並提供證據，但證據都不夠具說服力（只因為 number of supports 通常只有 2 ～4），必須要再有更多的資料來證明。但那時公司己經連薪水都快發不出來了，怎麼還有資金來提供後續的研究，所以跟醫院端這邊的合作案就先暫停了，真是可惜。
不過，我們在台北市研發補助進行結案報告時，把我們根據這一點點資料所找出的證據和新發現，都讓委員們讚嘆不己（或是說委員們也到了該擔心自己也會得這個癌症的時候了），還直接推薦我們進台北市研發補助的亮點企業的選拔，但我最後還是忍不住跟委員小酸了一下：「要是當初委員不要砍這麼多經費，我們一定可以有更多的發現。」
在整個資料分析的過程中，一個讓我想不透的就是第 616 號病患，她是歸類為晚期病患，同時這個癌症並不是她的第一個癌症，照理說，應該預期會看到很多 Somatic Variants，但分析的結果卻是在所有病患中偵測出最少的 Somatic variants（如圖三最下面一筆資料所示）。正常的晚期病患或有得過其它癌症的人都會有數千到數萬的 Somatic variants (未經任何過濾），但第616號病患卻只有 94 個。這一點讓我想不透，懷疑是不是資料處理出錯導致，經一再重覆確認每一步資料處理和分析過程，還是得到相同的結果，這個奇怪現象讓我在上台跟醫生們報告前好幾個夜晚都失眠。
最後，我還是硬著頭皮上場，特別提出這個奇怪現象跟醫生們討論，大家在場下也議論紛紛，除了講一些不相干的經驗分享，也開始跨不同部門間的溝通（或指責）。我就只能乖乖站在台上等醫師們討論完，最後，有位負責 BioBank 的技術人員終於釐清原因並說道：「因為當初第 616 病患取下的二個組織切片都是 Tumor 的，沒有合適的 Normal tissue，所以最後給出來的 Normal sample 實際上也是 Tumor tissue ! 」
晴天霹靂！聽到這個「美麗的錯誤」，我當場瞪大眼睛呆在台上！
造成我這失眠好幾個夜晚的原因竟然是這種問題！後來，我只能自己打圓場說：「這位技術人員應該是在替醫師們測試我們資料分析的能力，這次應該有過關吧？」
不過，也因為這次的考驗過關，讓醫師們願意跟我們多分享他們的經驗，同時，後續的會議也有越來越多醫師的參與。
對這位第 616 號病患也產生了感情，後來才得知這次的癌症檢體是這位病患的第4種癌症，而且這4種癌症都判定是原位癌，並非轉移。所以主任醫師就提議應該可以去分析她的 Germline Variant，找出她容易得癌症體質的可能致病變異。
另外，針對這次拿 2 個 Tumor Tissue 來跑 Mutect2 的 Somatic mode pipeline 的結果剛好可以讓我們來評估 Germline Effect 的影響程度及𧗠生出其它可能的分析手法，這也是本文想分享的主題，希望對大家的研究有所幫助。
實驗結果
簡單的說，我們目前手上有一位病人取下 2 份癌組織進行 WES 的定序資料，針對這 2 份定序資料進行 GATK HaplotypingCaller 和 Mutect2 (Tumor-Normal mode) 的分析流程分別可以得到其 Somatic variant 和 Germline variant（圖四），我們可以先來比較看看其相似性和不同點。
首先，我們先看 Germline variant 的個數，2個樣本皆偵測出 200 K 以上的變異，但彼此之間相同的變異約只有 8 成，其它 2 成左右是不同的。當然，這只是提供一個初步的數字，還沒有進行任何過濾（例如只看 PASS或進行其它 Quality check）。
其次，雖然是同一個人的癌組織，但因為 Tumor Heterogeneity，所以一定會看到不同的 Somatic mutation，由圖四顯示，分別得到 92 個和 296 個，而且這些 variant 是不重覆，也就是因為把另一個 tumor tissue 當作 normal sample，所以它們的交集才會是 0。然而這樣分析就會少掉在 2 個 samples 同時出現的 Somatic variants，接下來的問題是：「硬使用 Tumor-Normal mode 來分析，到底有多少同時出現的 Somatic Variant 會不見呢？」假設 Mutect2 與 HaplotypeCaller 一樣的能力，那可以估約有 8 成共同的 Somatic Variant，所以數量應該介於 500 ～1000 之間。
再者，約有 1/4 ~ 1/3 的Somatic variants 也會出現在 Germline 中，這點是蠻有趣的，因為 Somatic calling 時因為使用 Tumor — Normal 的分析，照理說己經濾掉所有的 Germline variant 才是，所以說 HaplotypCaller 也會偵測出一些 Somatic variants。
最後，當我們以為 Somatic variant 能被找出來是因為只出現在某一個樣本中而沒出現在另一個樣本，但實際上的結果並不是那樣單純，同時，我們之前擔心會有很多 2 個樣本的出現的 Somatic Variant 都被濾掉了，圖五顯示的結果讓我們可以不用這麼擔心這個問題。可以看到圖五是顯示一個 Somatic variants 在 S616_1 as Tumor 和 S616_2 as Normal 的情況下有被偵測出來，但在角色互換（S616_2 as Tumor 和 S616_1 as Normal）後沒被偵測出的實際 Alignment 狀況，可以發現這個變異應該是我們關切的 Somatic Variant，也實際上也出現在 2 個樣本中，但藉此本文提出跑 2 次 Tumor-Normal mode 分析後再把結果 Union 起來的作法，就可以大大減少因為缺 Normal Tissue 而無法有效降低 Germline Effect的問題，又不會因此而讓所有共同出現的 Somatic Variant 消失不見。
此外，Mutect2 提供 2 種不同分析方式，取決於是否有 Normal sample。也就是當有 Matched normal 的話，就可以跑 Somatic mode；若沒有 Matched normal，則建議可以跑 Tumor-only mode。當然，在沒有 Matched normal 的情況下，可以使用 germline resource 來濾掉 Germline variants。在這裡，感謝第 616 號病人的資料，讓我們可以進行 Tumor-normal 與 Tumor-only 的比較，由圖六可知，Tumor-only 偵測出的 Somatic Variants 個數會比 Tumor-normal 多 2 個 order 以上； 若只看 PASS 的話，也是多一個 order 以上的數量級。當然，這個第 616 號病人的 Somatic Variant 在 Tumor-normal mode 是會低估，但在 tumor-only mode 情況下是會高估。
大家可能注意到一點：「為什麼 Tumor-Normal 偵測到的 Somatic Variant 沒有全部都出現在 Tumor-only 的結果中呢？」
進一步觀察發現，這些缺少大部分出現在 low coverage regions，如圖七所示：
但這些線索只是所謂的 Phenotype，還沒找到這個問題的 Genotype。花了一整個下午追這個問題，後來終於在 GATK 的 source code 中找到一些線索：
原來 Mutect2 在 Tumor-Normal mode 和 Tumor-only mode 在一些參數設定上會有不同的 default value。例如圖八顯示 — af-of-alleles-not-in-resource 在 tumor-only mode 時的 default value 是 5e-8，而在 tumor-normal mode 時卻是 4e-3，因此造成 Tumor-only mode 沒有包含所有 Tumor-normal mode 偵測到的 Somatic Variants。這也再次提醒我們分析細節的重要性，不是每次都使用 Default setting 就可以得到好結果。
結論
總結一下本文的重點，想推薦大家在缺少 Matched Normal 的情況下（如 Bilateral Breast Cancer Patient 各取出左右乳房的癌組織；亦或定序深度足夠，則可以將單一樣本 Split 成 2 組資料），可以進行 Two-round Joint Somatic Calling pipeline：首先，各自將 2 組資料跑 Tumor-Normal 的 Mutect2 分析，之後再將 2 次的結果 Union 起來當作 First-round somatic candidate set來進行後續分析；若沒找到合適的，再使用 Tumor-only 的方式重新跑一次 Mutect2，一樣將 2 次的結果 Union 起來當作 Second-round somatic candidate set來進行後續分析。換句話說，First-round 是小而巧；Second-round 買一送一，提供我們針對後續不同的應該有更多的彈性選擇。
看到這裡，大家可能心中還是有些疑問：「到底第 616 號病人有多少真的 Somatic variants 呀？」
老實說，我也很想知道到底，但目前缺少 Normal Tissue的情況下，任何分析手法都只能去推估，沒有人有正確答案。除非，我們再提出另一個 IRB 申請，但經費呢？想到又要被委員羞辱一番，經費又要被砍一大半，就沒有動力去追出這個問題的答案了。
最後，癌症這個疾病一點都不美麗，但相信每次的資料分析是可以產生美麗的結果，即早讓病人可以得到更合適的治療，在此將本文獻給在天堂的淑珍學妹。
[1] https://en.wikipedia.org/wiki/Reference_genome
[2] https://en.wikipedia.org/wiki/Tumour_heterogeneity
[3] http://www.labmed.org.tw/knowledge_1.ASP?mno=90
[4] https://www.nature.com/articles/s41598-020-60559-5
[5] http://dict.idioms2.moe.edu.tw/idiomView.jsp?ID=-13978&webMd=2&powerMode=1&la=0
[6] https://gatk.broadinstitute.org/hc/en-us/articles/360035889791--How-to-Call-somatic-mutations-using-GATK4-Mutect2-Deprecated-#1
[7] https://github.com/broadinstitute/gatk/blob/3021e6924aeb84d9f3b333e5298abb1ec27d350a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2.java
