基於 MaxCompute 的極速的基因測序分析

基因、測序、分析

基因，生命的基本因素，是人類和其他生物的基礎遺傳物質。人有 23 對染色體，總共記錄了大約 3Gb 個鹼基（這裡的 b 是 base，即鹼基，可不是 bit，參考這裡），每個位置上的鹼基可能是 ATCG 中的一個。簡單理解起來，就是有了這 3Gb 長的字元串，就能克隆一個你。基因測序，就是用化學和物理的方法，把你身體里這 3Gb 字元串檢測出來。

當然，由於受當前測序技術的限制，我們並不能一次性測得一個完整的 3Gb 字元串，而是無數個 150bp 左右長度的小碎片。把這無數個小碎片重新組合還原成 3Gb 的長字元串的過程，叫全基因組組裝。人類基因組計劃乾的就是這個組裝拼圖的事情，到了 2003 年，基本上算是拼完了。於是就有了一個標準的 3Gb 人類字元串，業界稱其為『人類基因組參考序列』，也常被簡稱為『參考序列』。

基因測序分析，就是檢測每個個體身上的這 3Gb 字元串，然後再跟標準的字元串（參考序列）做比對，來看是不是哪裡發生了變異，是否有已知的遺傳疾病風險。

現在隨著基因測序技術的成熟，成本在飛速下降，所謂『舊時王謝堂前燕，飛入尋常百姓家』簡直就是指日可待的事情。

於是越來越多的基因數據需要計算分析。

原來的分析計算流程

測序出來的數據是一大堆長度 150bp 的小碎片，但由於已經有了完整的人類參考序列的拼圖，那麼在這個拼圖上尋找位置要比還原拼圖的過程容易很多。多年來，從輸入碎片數據比對到標準的 3Gb 參考序列，再到變異的檢測，形成了如 BWA、Samtools/Picard、GATK 等為單機運行優化良好的業界公認軟體。

將這些軟體組合使用起來可以形成一個基因組數據分析的流程。華大基因基於阿里雲 ECS 的 BGI Online 平台、七橋在 Amazon 和 Google 雲上架構的 IGOR 分析平台，都有這樣的業界金標準流程（下圖為七橋，1.7G 左右的輸入數據跑了 1 小時 48 分）。

但是上面的分析只是外顯子測序分析（WES），所謂外顯子就是那些會直接影響我們外在性狀和健康與否的序列，它的長度只佔了全基因組總長的 1%。因此需要測序的數據量比全基因組測序分析（WGS）少得多。使用 WES 的原因和意義，一方面是因為我們對這 3Gb 的人類基因組所蘊含的全面信息所知甚少，這 1% 的序列是我們目前能夠很有效進行解讀的位點，同時，也是由於它直接影響著我們的性狀和健康，所以在很多疾病的研究中，通常只檢測和分析這些位點和區域；另一方面也是因為全基因組測序目前成本相對較高，數據量巨大，單機計算時間基本要幾天。

基於 MaxCompute 的方案

數據表達

來看基因數據的幾個重要表達。

測序出來的原始輸入數據 FastQ 格式（截取了兩個 record），原本就是四行一個 record 的文本：

@ERR194147.1 HSQ1004:134:C0D8DACXX:1:1104:3874:86238/1GGTTCCTACTTCAGGGTCATAAAGCCTAAATAGCCCACACGTTCCCCTTAAATAAGACATCACGATGGATCACAGGTCTATCACCCTATTAACCACTCACG+CC@FFFFFHHHHHJJJFHIIJJJJJJIHJIIJJJJJJJJIIGIJJIJJJIJJJIJIJJJJJJJJJJIJHHHHFFFDEEEEEEEEDDDCDDEEDDDDDDDDD@ERR194147.2 HSQ1004:134:C0D8DACXX:2:2104:2852:75174/1ACTTCAGGGTCATAAAGCCTAAATAGCCCACACGTTCCCCTTAAATAAGACATCACGATGGATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTC+@BBFDFFFHFFHHIHIJJJJFIHHFHFHJCIHFHIJJJJJJJIJIJIJJIIHIJJJJJJJBEGIGHIHGHHHEFCDFFEDEEDEEDDD?CCCDDDDDDDDC

測序產生的中間及最終結果的 VCF 格式（截取兩個 record），也是一個多行表頭，然後每行一個數據 record 的形式。

##fileformat=VCFv4.1##...##other headers#CHROM POS ID REF ALT QUAL FILTER INFO1 10019 rs376643643 TA T . . OTHERKG;R5;RS=376643643;RSPOS=10020;SAO=0;SSR=0;VC=DIV;VP=0x050000020001000002000200;WGT=1;dbSNPBuildID=1381 10054 rs373328635 CAA C,CA . . NOC;OTHERKG;R5;RS=373328635;RSPOS=10055;SAO=0;SSR=0;VC=DIV;VP=0x050000020001000002000210;WGT=1;dbSNPBuildID=138

用表來表達描述完全沒問題。當然，PairEnd 的 FastQ 總是成對出現，那麼就編排成 8 列的表就好了。VCF 格式的多行表頭也要稍微處理一下，不過這都不是事兒。

問題可分割

這塊涉及到生物學的專業知識，要描述清楚比較困難。基本的想法，是將 3Gb 標準數據切成若干份，然後分散式執行每一份的計算。理論上，切割肯定會對計算的精度帶來影響。但是影響多大要試過才知道。實踐中，也有些辦法來彌補這些影響：比如適當擴大 reference 的 interval，或者通過數據重疊來緩解。

簡言之，跟單機計算結果對比，一致率 98.8%，這個數字可以接受。

MaxCompute Streaming 作業

如果數據如何在分散式系統里流動都想得很清楚的話，那麼最大的障礙其實已經不在，剩下的就都是工程問題。

MaxCompute 支持 Streaming 作業，允許用戶使用任何語言開發 MR 作業，開啟隔離環境之後甚至可以用在公共雲環境。用 Streaming 作業，可以非常方便的將 BWA、samtools、GATK 集成在一個 MR 作業當中：

這個作業，運行時間不到 3 個小時。

結語

或許很快，癌症患者可以當天拿到靶向藥物是否有療效的根據；新生病兒的早期診斷當天就可以排除是否遺傳疾病；基因測序分析就像今天驗血一樣成為普通人司空見慣的體檢手段。

或許將來，MaxCompute 可以更加發揮計算的能力，輔助挖掘 DNA 上目前 99% 還未知的領域。

本來選自阿里雲大數據產品專家「隱林」，擅長MaxCompute、機器學習、分散式、可視化、人工智慧等大數據領域。