基於眾核協處理器的三級流水序列比對方法
2024-01-21 17:53:15 1
基於眾核協處理器的三級流水序列比對方法
【專利摘要】本發明公開了一種基於眾核協處理器的三級流水序列比對方法，目的是提高序列比對軟體的比對速度。技術方案是：採用MIC眾核協處理器使用多線程並行進行序列比對，並將MIC上序列比對過程中從主存讀取序列到MIC、比對序列和將比對結果寫回主存這三個串行步驟採用三級流水方式，即在進行每一輪序列比對的同時，讀取下一組比對所需要的序列，同時將上一組比對的結果寫回主存，使得讀寫操作和比對操作並行進行。本發明實現了序列讀取，序列比對和比對結果寫回的三個主要過程的並行執行，提高了比對效率，降低了比對時間，與兩路八核CPU相比，本發明可加速比對過程2.3倍以上，且避免拷貝大量內存空間，提高了程序的時空效率。
【專利說明】基於眾核協處理器的H級流水序列比對方法
【技術領域】
[0001] 本發明涉及生物信息領域序列比對的方法，尤其指一種基於眾核協處理器的序列 比對方法。
【背景技術】
[0002] 分子生物學是從分子水平上研究生命現象物質基礎的學科，通過研究生物分子的 結構、功能和合成等方面的原理，從而使生物體的功能和性狀在前所未有的分子細節上得 到詳盡的分析和理解，進而更加科學嚴謹地闡明生命現象的本質。
[0003] 在分子生物學研究中，DNA的序列分析是進一步研究和改造目的基因的基礎。 DNA(脫氧核糖核酸）是一種生物大分子，一共分為四種鹼基，記為A、T、C、G，該些大分子的 排列順序決定了某種遺傳指令，該些遺傳指令是建構細胞內其他的化合物，如蛋白質與核 糖核酸的需要。帶有蛋白質編碼的DNA片段稱為基因，即遺傳物質，是DNA分子上具有遺傳 信息的特定核巧酸序列。基因經過轉錄、翻譯，最終產生結構和功能各異的、表現生物體性 狀的蛋白質。
[0004]DNA序列分析的基礎是對DNA分子進行測序，即確定DNA分子中A、T、C、G四種鹼 基的排列順序。當前的DNA測序技術，一次實驗最多只能直接測得不大於5000個鹼基的排 列順序，形成多個DNA短序列（稱為read)。而一般生物的基因組鹼基數目都十分巨大，女口 人類基因組總長約為30億個鹼基對。該樣，絕大多數生物的基因組都不能通過實驗手段一 次性獲得，而必須藉助於計算機技術進行後續拼接得到完整的基因組。
[000引序列比對是目前廣泛使用的DNA序列分析方法，它是將測序得到的read短序列直 接與拼接完成的參考基因組進行比對，確定read在參考基因組中是否出現W及出現的具 體位置。通過序列比對進行DNA序列分析，避免了對目標基因組進行組裝，可W很大程度上 節省序列分析的時間和工作量，提高序列分析的效率。
[0006] 由於比對時read數量都較大，無法一次性全部存放到計算機主存中。所W目前常 用的DNA序列比對方法均按照W下步驟進行：
[0007] 步驟1 ;根據計算機主存可用空間大小，將read平均分為若干組，每一組所佔空間 大小不超過計算機主存容量；
[0008] 步驟2 ;從磁碟上讀取一組read到主存中；
[0009] 步驟3 ;對讀取到主存中的read逐個進行比對；
[0010] 步驟4 ;將read比對結果寫回磁碟；
[0011] 步驟5 ;檢查磁碟中是否還存在未比對的序列，如果存在，返回步驟2 ;如果不存 在，結束比對過程。
[0012]目前比對主要使用的運算器件為計算機中央處理器CPU或者圖形處理器GPU。
[0013] 雖然與進行序列組裝相比，序列比對可W節省大量時間，但是目前廣泛使用的 基於CPU的串行序列比對方法的速度仍比較慢，如在配備兩路八核Intel2.4GHzCPU 的常用伺服器上，米用李恆在論文《Fastandaccurateshortreadalignmentwith Burrows-WheelerTransform》中公布的基於BW(Burrows-Wheeler)變換的序列比對方法， 對長度為100個鹼基的8千萬條序列進行比對，需要花費一天W上時間，很難滿足後序的序 列分析對於時間的要求，更是無法滿足時效性要求較高的臨床需求。
[0014] 基於CPU的並行序列比對方法使得多個線程能夠並行進行序列比對，有效地提高 了序列比對的速度。但是目前絕大部分研究機構使用的是單節點伺服器，CPU計算能力十 分有限。而隨著測序技術的發展，特別是新一代高通量測序技術的出現，單位時間內產生的 read數量翻了幾翻，基於CPU的並行序列比對軟體也已經很難有效處理如此大量的read。
[0015] 基於GPU的序列比對軟體，利用GPU具有大量計算核也的特性，使用其對序列比對 進行加速，有效地增強了伺服器的計算能力，與基於CPU的並行序列比對方法相比，進一步 提高了速度。GPU指令集設計與CPU相比較為簡單，但是能夠快速處理簡單的浮點和整型 計算。而目前廣泛使用的李恆在論文《Fastandac州rateshortreadalignmentwith Burrows-WheelerTransform》中公布的基於BW炬urrows-Wheeler)變換的序列比對方法運 算過程複雜，程序分支多，當GPU中的一個核也遇到分支時，與其同組的其他核也均要等待 該分支處理完畢才能繼續並行執行，很大程度上降低了序列比對的效率。
[0016]MIC(ManyIntegratedCore)是Intel公司開發的一款眾核協處理器，具有與傳統 X86CPU兼容的指令集，能夠快速處理分支判斷等複雜指令。而且，每個MIC協處理器配備 了 50個W上的計算核也，每個計算核也可W啟動4個硬體線程，並行規模可達200個線程 W上，其核也主頻約為1.IGHz，並且包含512位寬度的向量處理單元，配備6GBW上片上存 儲空間，其單卡雙精度浮點計算峰值超過ITFlops。MIC與GPU相比，在峰值計算能力相當 的情況下，能夠有效地處理複雜指令，提高複雜程序的運行效率。因此，本發明採用MIC眾 核協處理器加速序列比對。
[0017] 但是MIC仍然存在一定局限，其計算核也只能訪問它自身的片上存儲空間，無法 直接訪問計算機主板上的主存空間，只能由CPU來負責主機與MIC之間的數據傳輸。一種 解決方案是在運行前將程序和所需的數據事先拷貝到MIC上，然後登錄上MIC上的操作系 統，啟動程序，程序的輸出寫在MIC上，等運行完畢後再將結果手動拷貝到主存。W人類基 因組為例，序列比對過程的輸入、輸出，W及中間文件所佔空間超過10GB，但是目前MIC卡 上的存儲空間無法滿足需求，所W該方案困難非常大。另一種解決方案是從CPU端啟動程 序，在程序運程過程中由CPU將序列比對的輸入數據拷貝到MIC上，MIC完成比對後，再由 CPU將比對結果拷貝到主存，並將下一批序列拷貝到MIC上，如此循環，直到所有序列比對 完畢。但是CPU與MIC之間頻繁的數據傳輸佔用了大量的程序運行時間，很大程度上影響 了程序的運行效率。如何提高計算效率，降低數據傳輸開銷是利用MIC加速序列比對的難 點。目前還沒有利用MIC進行序列比對的技術方案的公開報導。
【發明內容】
[0018] 本發明要解決的技術問題是提出一種基於MIC眾核協處理器的H級流水序列比 對方法，提高序列比對軟體的比對速度。
[0019] 本發明的技術方案是：採用MIC眾核協處理器使用多線程並行進行序列比對，並 將MIC上序列比對過程中從主存讀取序列到MIC、比對序列和將比對結果寫回主存該H個 串行步驟採用H級流水方式，即在進行每一輪序列比對的同時，讀取下一組比對所需要的 序列，同時將上一組比對的結果寫回主存，使得讀寫操作和比對操作並行進行。
[0020] 具體技術方案如下：
[0021] 主要變量定義：
[002引 M_CPU;計算機主存可用空間大小。
[0023] M_DNA;DNA短序列所佔空間大小。
[0024] M_MIC;MIC上可用空間大小。
[002引 M_REF;參考基因組大小。
[002引 M_SEQ;MIC上每塊緩存空間大小。
[0027] 步驟1 ;CPU根據計算機主存可用空間大小M_CPU，W及DNA短序列所佔空間大小 M_DNA，將DNA短序列即read平均分為L組，L為正整數，L=「M_DNA/M_CPU-1，''「xl"表示對 "X"向上取整；
[0028] 步驟2 ;在MIC上聲明H個指針變量；Seqs_ptr、Read_ptr、和Write_pt;r，並根據 MIC上可用空間大小M_MICW及參考基因組大小M_REF，在MIC上分別為H個指針分配同樣 大小的緩存空間，MIC上每塊緩存空間大小M_SEQ= (M_MIC-M_REF) /3,其中Seqs_ptr指向 的空間存儲當前一組正在比對的序列，ReacLptr指向的空間存儲下一組將要比對的序列， Write_pt;r指向的空間存儲上一組序列比對的結果；
[0029] 步驟3 ;CPU初始化循環變量i為零；
[0030] 步驟4 ;CPU將磁碟中L組read中的第i組讀入主存中；
[0031] 步驟5 ;CPU根據MIC上Seqs_ptr指向的空間的大小M_SEQ，將讀入主存中的read 平均分為M小組，M為正整數，M=「(M_DNA/L)/M_SEQ1;
[0032] 步驟6 ;CPU將循環變量m置為零；
[0033] 步驟7 ;CPU將主存中M個小組read中的第m組讀到MIC的Seqs_ptr指向的空間 內；
[0034] 步驟8 ;根據MIC上可用計算核也數Core_MIC，W及MIC上每個計算核也支 持的最大硬體線程數化read_MIC，在MIC上同時啟動化2 (N〉0,為整數，N= (Core_ MIC-1) *T虹eacLMIC，其中MIC卡需要保留一個核也處理主存與MIC卡之間數據調度）個線 程，線程編號為0到化1，化2個線程並行執行W下步驟：
[003引步驟8. 1 ;第0到第N-1號線程並行比對Seqs_ptr對應空間中的read,並將比較結 果寫入Seqs_ptr對應空間中；比對方法採用李恆在論文《Fastandaccurateshortread alignmentwithBurrows-WheelerTransform》中公布的基於BW炬urrows-Wheeler)變換 的方法，該一步與【背景技術】步驟3的區別在於是在MIC上由N個線程並行地對進行比對，而 不是在CPU上由單個線程完成，將Seqs_ptr對應空間中的所有read比對完畢後，轉到步驟 9 ；
[0036] 步驟8. 2 ;第N號線程將循環變量m加1，判斷m是否等於M，如果m不等於M，執行 步驟8. 2. 1，如果m等於M，結束第N號線程，轉到步驟9 ;
[0037] 步驟8. 2. 1 ;第N號線程將主存中M個小組read中的第m組讀到MIC的ReacLptr 對應空間內，讀取完畢後，轉到步驟9 ;
[0038] 步驟8. 3 ;第化1號線程判斷Write_ptr對應空間是否為空，如果Write_ptr對應 空間不為空，執行步驟8. 3. 1，如果Write_ptr對應空間為空，結束第化1號線程，轉到步驟 9 ；
[0039] 步驟8.3. 1 ;第化1號線程將Write_ptr對應空間中的read比對結果寫回主存， 寫完成後轉到步驟9;
[0040] 步驟9 ;同步第0到第化1號線程，同步完成後，在MIC上的多線程部分執行完畢， W下步驟為單線程執行；
[0041] 步驟10 ;MIC進行指針交換，在MIC上聲明臨時指針tmp_ptr，將Seqs_ptr值賦給 tmp_ptr，將Read_ptr值賦給Seqs_ptr，將Write_pt;r值賦給Read_ptr，將tmp_ptr值賦給 Write_pt;r，將tmp_ptr值置為空。
[0042] 步驟11 ;MIC判斷Seqs_ptr對應空間是否為空，如果不為空，轉步驟8,如果為空， 執行步驟12;
[0043] 步驟12 ;MIC將Write_ptr對應空間中的read比對結果寫回主存；
[0044] 步驟13 ;CPU將內存中第i大組read比對的結果寫回磁碟，並清空相應內存空間； [004引步驟14 ;CPU將循環變量i加1;
[0046] 步驟15 ;CPU判斷i是否等於L如果i不等於L轉步驟4,如果i等於L執行步 驟16;
[0047]步驟 16;MIC釋放Seqs_ptr、Read_ptr、Write_pt;r指向的空間；
[004引步驟17;結束比對。
[0049] 採用本發明可W達到W下技術效果：
[0050] 本發明通過多線程並行技術，利用新型運算器件MIC眾核協處理器進行序列比 對，其中
[005。步驟8在MIC上利用化2 (N〉0)個線程並行進行序列比對，提高了比對速度，在多 線程下可W取得接近線性的加速比，並實現了序列讀取，序列比對和比對結果寫回該H個 主要過程的並行執行，提高了比對效率，降低了比對時間，與兩路八核CPU相比，本發明可 加速比對過程2.3倍W上。
[0052] 步驟9在MIC上通過交換指針的方式實現數據交換，而不是直接對變量進行賦值， 避免拷貝大量內存空間，提高了程序的時空效率。
【專利附圖】
【附圖說明】
[0053] 圖1是本發明總體流程圖；
[0054] 圖2是對圖1中步驟8的分解示意。
【具體實施方式】
[00巧]國防科大採用配備兩路八核2. 4GHzCPU和一塊57核1.IGHzMIC卡的伺服器作 為環境，伺服器硬碟大小為43TB，內存大小為132GB，MIC卡片上存儲空間大小為6GB，輸入 數據為人類基因組，參考基因組佔空間大小為3GB，DNA短序列佔空間大小為240GB，包括8 千萬條序列，驗證本發明的效果：
[0056] 如圖1所示，具體實施步驟如下：
[0057] 步驟1 ;CPU根據計算機主存可用空間大小M_CPU= 45GB(作業系統及其他服務佔 用了一定內存，另外還需要預留一部分內存作程序運行時使用，所W可用內存大小為45GB， 小於安裝內存大小132GB)，W及DM短序列所佔空間大小M_DNA= 240GB，將DM短序列即read平均分為L=6組，L為正整數，L=「IVI_DNA/M_CP川=「240GB/45G1=6;
[0058] 步驟2;在MIC上聲明H個指針變量；Seqs_ptr、Read_ptr、和Write_pt;r，並根據 MIC上可用空間大= 4. 5GB(MIC內存大小為6GB，其中1. 5GB空間用於存儲MIC片 上作業系統W及運行時使用）W及參考基因組大小M_REF= 3GB，在MIC上分別為H個指針 分配同樣大小的緩存空間，每塊空間大小1_569 = (M_MIC-M_RE巧/3 = 0. 5GB，其中Seqs_ ptr指向的空間存儲當前一組正在比對的序列，ReacLptr指向的空間存儲下一組將要比對 的序列，Write_pt;r指向的空間存儲上一組序列比對的結果；
[0059] 步驟3;CPU初始化循環變量i為零；
[0060] 步驟4 ;CPU將磁碟中L= 6組read中的第i組讀入主存中；
[0061] 步驟5;CPU根據MIC上Seqs_ptr指向的空間的大小M_SEQ=0. 5GB，將讀入主存 中的read平均分為M=80小組，M為正整數，M=「(M_DNA/L)/M_SEQM(240GB,'巧化5GB1 =80;
[0062] 步驟6 ;CPU將循環變量m置為零；
[0063] 步驟7;CPU將主存中M=80個小組read中的第m組讀到MIC的Seqs_ptr指向 的空間內；
[0064] 步驟8如圖2所示；根據MIC上可用計算核也數Core_MIC= 57,W及MIC上每個 計算核也支持的最大硬體線程數化reacLMIC= 4,在MIC上同時啟動化2 = 226(N〉0,為整 數，N= (Core_MIC-l)*!'虹ead_MIC=巧7-1)*4 = 224,其中MIC卡需要保留一個核也處理 主存與MIC卡之間數據調度）個線程，線程編號為0到化1 = 225,化2 = 226個線程並行 執行W下步驟：
[006引步驟8. 1;第0到第N-1 = 223號線程並行比對Seqs_ptr對應空間中的read,並將 比較結果寫入Seqs_ptr對應空間中；比對方法採用李恆在論文《Fastandaccurateshort readalignmentwithBurrows-WheelerTransform》中公布的基於BW炬urrows-Wheeler) 變換的方法，該一步與【背景技術】步驟3的區別在於是在MIC上由N= 224個線程並行地對 進行比對，而不是在CPU上由單個線程完成，Seqs_ptr對空間中的所有read比較完畢後， 轉到步驟9;
[0066] 步驟8. 2 ;第N= 224號線程將循環變量m加1，判斷m是否等於M= 80,如果m不 等於M= 80,執行步驟8.2. 1，如果m等於M= 80,結束第N= 224號線程，轉到步驟9 ;
[0067] 步驟8. 2. 1;第N=224號線程將主存中M=80個小組read中的第m組讀到MIC 的ReacLptr對應空間內，讀取完畢後，轉到步驟9 ;
[006引步驟8. 3;第化1 = 225號線程判斷Write_ptr對應空間是否為空，如果Write_ptr對應空間不為空，執行步驟8. 3. 1，如果Write_ptr對應空間為空，結束第化1 = 225號 線程，轉到步驟9;
[0069] 步驟8.3. 1;第化1=225號線程將Write_ptr對應空間中的read比對結果寫回 主存，寫完成後轉到步驟9;
[0070] 步驟9;同步第0到第化1=225號線程，同步完成後，在MIC上的多線程部分執 行完畢，W下步驟為單線程執行；
[0071]步驟10 ;MIC進行指針交換，在MIC上聲明臨時指針tmp_ptr，將Seqs_ptr值賦給 tmp_ptr，將Read_ptr值賦給Seqs_ptr，將Write_pt;r值賦給Read_ptr，將tmp_ptr值賦給 Write_pt;r，將tmp_ptr值置為空。
[0072] 步驟11 ;MIC判斷Seqs_ptr對應空間是否為空，如果不為空，轉步驟8,如果為空， 執行步驟12 ;
[0073] 步驟12 ;MIC將Write_ptr對應空間中的read比對結果寫回主存；
[0074] 步驟13 ;CPU將內存中第i大組read比對的結果寫回磁碟，並清空相應內存空間； [007引步驟14 ;CPU將循環變量i加1 ;
[0076] 步驟15 ;CPU判斷i是否等於L= 6,如果i不等於L= 6,轉步驟4,如果i等於L =6,執行步驟16 ;
[0077]步驟 16;MIC釋放Seqs_ptr、Read_ptr、Write_pt;r指向的空間；
[007引步驟17;結束比對。
[0079] 統計運行數據後發現，與僅使用兩路八核CPU並行進行序列比對相比，本發明可 W將序列比對速度提高2. 3倍W上。
【權利要求】
1. 一種基於眾核協處理器的三級流水序列比對方法，其特徵在於包括以下步驟： 步驟I :CPU根據計算機主存可用空間大小M_CPU，以及DNA短序列所佔空間大小M_ DNA，將DNA短序列即read平均分為L組，L為正整數，L=「M_DNA/M_CPUl, "「χΓ表示 對"X"向上取整； 步驟2 :在MIC即眾核協處理器上聲明三個指針變量：Seqs_ptr、Read_ptr和Write_ ptr，並根據MIC上可用空間大小M_MIC以及參考基因組大小M_REF，在MIC上分別為三個 指針分配同樣大小的緩存空間，MIC上每塊緩存空間大小M_SEQ = (M_MIC-M_REF)/3,其中 Seqs_ptr指向的空間存儲當前一組正在比對的序列，Read_ptr指向的空間存儲下一組將 要比對的序列，Write_ptr指向的空間存儲上一組序列比對的結果； 步驟3 :CPU初始化循環變量i為零； 步驟4 :CPU將磁碟中L組read中的第i組讀入主存中； 步驟5 :CPU根據MIC上Seqs_ptr指向的空間的大小M_SEQ，將讀入主存中的read平 均分為M小組，M為正整數，M=「(M_DNA/L),/M_SEQl; 步驟6 :CPU將循環變量m置為零； 步驟7 :CPU將主存中M個小組read中的第m組讀到MIC的Seqs_ptr指向的空間內； 步驟8 :根據MIC上可用計算核心數Core_MIC，以及MIC上每個計算核心支持的最大硬 件線程數Thread_MIC，在MIC上同時啟動N+2個線程，線程編號為O到N+l，N>0,為整數，N =(Core_MIC-l)*Thread_MIC，其中MIC卡需要保留一個核心處理主存與MIC卡之間數據調 度，N+2個線程並行執行以下步驟： 步驟8. 1 :第O到第N-I號線程並行比對Seqs_ptr對應空間中的read，並將比較結果 寫入Seqs_ptr對應空間中；比對方法採用李恆在論文《Fast and accurate short read alignment with Burrows-Wheeler Transform》中公布的基於 BW 變換的方法，將 Seqs_ptr 對應空間中的所有read比對完畢後，轉到步驟9 ; 步驟8. 2 :第N號線程將循環變量m加1，判斷m是否等於M，如果m不等於M，執行步驟 8. 2. 1，如果m等於M，結束第N號線程，轉到步驟9 ; 步驟8. 2. 1 :第N號線程將主存中M個小組read中的第m組讀到MIC的Read_ptr對 應空間內，讀取完畢後，轉到步驟9 ; 步驟8. 3 :第N+1號線程判斷Write_ptr對應空間是否為空，如果Write_ptr對應空間 不為空，執行步驟8. 3. 1，如果Write_ptr對應空間為空，結束第N+1號線程，轉到步驟9 ; 步驟8. 3. 1 :第N+1號線程將Write_ptr對應空間中的read比對結果寫回主存，寫完 成後轉到步驟9 ; 步驟9 :同步第0到第N+1號線程，同步完成後，在MIC上的多線程部分執行完畢，以下 步驟為單線程執行； 步驟10 :MIC進行指針交換，在MIC上聲明臨時指針tmp_ptr，將Seqs_ptr值賦給tmp_ ptr，將 Read_ptr 值賦給 Seqs_ptr，將 Write_ptr 值賦給 Read_ptr，將 tmp_ptr 值賦給 Write_ptr，將 tmp_ptr 值置為空。 步驟11 :MIC判斷Seqs_ptr對應空間是否為空，如果不為空，轉步驟8,如果為空，執行 步驟12 ; 步驟12 :MIC將Write_ptr對應空間中的read比對結果寫回主存； 步驟13 :CPU將內存中第i大組read比對的結果寫回磁碟，並清空相應內存空間； 步驟14 :CPU將循環變量i加1 ; 步驟15 :CPU判斷i是否等於L，如果i不等於L，轉步驟4,如果i等於L，執行步驟16 ; 步驟 16 :MIC 釋放 Seqs_ptr、Read_ptr、Write_ptr 指向的空間； 步驟17 :結束比對。
【文檔編號】G06F9/38GK104375807SQ201410745667
【公開日】2015年2月25日 申請日期:2014年12月9日 優先權日:2014年12月9日
【發明者】廖湘科, 朱小謙, 崔英博, 彭紹亮, 鄒丹, 王恆, 朱敏, 劉欣, 王海強, 高明 申請人:中國人民解放軍國防科學技術大學