簡介這是一款開源的Matlab工具箱面向蛋白質(zhì)多序列比對的協(xié)同進化分析采用簡化的BSD許可證發(fā)布并且不依賴Matlab自帶的生物信息工具箱。工具箱集成了統(tǒng)計耦合分析、直接耦合分析、互信息、顯式似然等多種主流協(xié)同進化計算方法每種方法均以獨立函數(shù)形式封裝用戶可按需調(diào)用便于對不同算法的計算結(jié)果進行橫向比較同時工具內(nèi)還附帶多種生物信息學(xué)矩陣檔案為后續(xù)的深度分析提供基礎(chǔ)。整個壓縮包共包含116個文件其中核心為103個Matlab腳本另有支持Windows、macOS與Linux平臺的動態(tài)鏈接庫、可執(zhí)行文件以及蛋白質(zhì)結(jié)構(gòu)、二級結(jié)構(gòu)、序列等示例數(shù)據(jù)文件并包含C源碼以便有需要的用戶進行底層修改。壓縮包體積僅為264KB輕量小巧下載和部署都非常方便。示例數(shù)據(jù)中提供了蛋白質(zhì)序列與二級結(jié)構(gòu)文件能夠幫助使用者快速驗證整個分析流程。目前已有13958人學(xué)習(xí)或下載特別適合從事蛋白質(zhì)結(jié)構(gòu)預(yù)測、進化保守性研究以及相關(guān)算法開發(fā)的科研人員和Matlab開發(fā)者。1. 從多序列比對里找“一起變”的殘基AKtoolbox 是什么、為什么值得拆我剛接觸多序列比對的時候總覺得分析到保守位點就完事了。后來在某個模擬項目X里做功能殘基篩選被一個問題卡住兩個位點單獨看都不夠保守但它們的變異模式總是成對出現(xiàn)這種信號在傳統(tǒng)的序列Logo圖里根本看不見。AKtoolbox就是為解決這個問題而生的開源軟件一個純MATLAB工具箱專門做多序列比對的協(xié)同進化分析。它把FASTA讀取、位點過濾、互信息/統(tǒng)計耦合計算、置換檢驗、結(jié)果導(dǎo)出串成一條完整鏈路尤其適合手里沒有Python生信環(huán)境、卻已經(jīng)在用MATLAB做數(shù)據(jù)分析的工程師和研究生。如果你手頭已經(jīng)有一份比對好的FASTA從敲第一行命令到拿到殘基配對得分矩陣通常用不了一個上午。2. 協(xié)同進化分析為什么選 MATLAB先想清楚提取的是哪一種“聯(lián)動信號”2.1 協(xié)同進化信號來自哪里MSA 給我們的原始輸入是什么多序列比對本質(zhì)上是一張字符表每一行是一條序列每一列是比對上的一個位點字符可能是氨基酸殘基也可能是-這樣的空位。我們看一個位點是否保守通常會統(tǒng)計這個位點上出現(xiàn)頻率最高的殘基占多大比例但這是單點統(tǒng)計它丟掉了一個關(guān)鍵信息——位點之間是否存在“組合偏好”。舉個最簡單的例子。假設(shè)一個假想MSA里有100條序列位點12和位點80單獨看都各有約一半序列出現(xiàn)兩類殘基單點保守性得分都不高?？僧?dāng)你把兩個位點放在一起看時會發(fā)現(xiàn)規(guī)律非常強只要位點12是K位點80幾乎一定是D只要位點12是E位點80幾乎一定是R。這說明兩個位點根本不是獨立變異的它們很可能在三維結(jié)構(gòu)上靠在一起或者共同受到同一選擇壓力一個變了另一個必須跟著補償。協(xié)同進化分析想提取的就是這種成對非隨機聯(lián)動。它把每個位點當(dāng)作一個離散變量把每一對位點看作一個二維列聯(lián)表問題判斷觀測到的聯(lián)合分布是否顯著偏離兩個邊緣分布的乘積。理解這一點之后后面所有預(yù)處理參數(shù)都有了方向我們希望在MSA里盡可能去掉噪聲留下真正能反映“一起變”的信號。2.2 工具箱里的三類核心算法MI、SCA 與擾動統(tǒng)計AKtoolbox把主流的協(xié)同進化算法封裝成了幾個統(tǒng)一入口你不需要把每個算法從頭推導(dǎo)一遍但必須清楚它們各自的輸出含義否則很容易拿錯結(jié)果去解釋生物學(xué)問題。第一類是互信息也就是 MI。它比較的是每個位點對聯(lián)合分布p(x,y)與兩個位點邊緣分布乘積p(x)p(y)的差異差異越大互信息值越高。MI實現(xiàn)最簡單、計算最快適合做初步篩選。但它的缺點也很明顯天然偏向保守位點。因為保守位點本身只有少數(shù)幾個字符類別概率估計相對穩(wěn)定一旦出現(xiàn)少量共變殘差就會被放大而高度可變的位點類別太多概率被攤薄即使存在真實共變也容易被淹沒。第二類是統(tǒng)計耦合分析也就是 SCA。它不會直接給出二維的“位點A到位點B”的分?jǐn)?shù)而是先對每個位點做熵加權(quán)頻率歸一化然后對加權(quán)頻率矩陣做特征分解再通過譜重構(gòu)得到位點之間的耦合關(guān)系。SCA對保守位點有一定懲罰所以更容易找出結(jié)構(gòu)域內(nèi)部的功能模塊適合分析較大范圍的殘基網(wǎng)絡(luò)。但它比較“吃數(shù)據(jù)”要求MSA序列數(shù)和序列多樣性都足夠如果只拿20條序列去跑SCA出來的特征向量基本是噪聲。第三類可以叫擾動統(tǒng)計法風(fēng)格上和ELSC、McBASC這一類方法接近。它的思路是把每條序列當(dāng)作一次自然擾動用回歸或協(xié)方差模型分離出其他位點對當(dāng)前位點變異的貢獻然后做顯著性檢驗。這類方法比MI更“條件化”能消除一部分背景相關(guān)但計算量最大而且對MSA里的異常序列特別敏感如果里面混進幾條和外類群序列差異極大的序列它們可能主導(dǎo)整個協(xié)方差估計。下面這張表可以幫助選型算法類型輸出形式適合場景主要風(fēng)險互信息 MI位點對之間的非負(fù)分?jǐn)?shù)快速初篩找候選配對偏向保守位點gap處理不當(dāng)容易假高統(tǒng)計耦合 SCA位點譜權(quán)重與重構(gòu)的耦合矩陣找結(jié)構(gòu)域協(xié)同模塊、功能界面相關(guān)殘基需要足夠多樣本小樣本下不穩(wěn)定擾動統(tǒng)計殘差或檢驗統(tǒng)計量中長序列、需要定量顯著性時對異常序列敏感計算耗時長2.3 為什么在 MATLAB 里做這件事以及 AKtoolbox 的定位很多人會先問現(xiàn)在Python生態(tài)里做生物信息學(xué)已經(jīng)很成熟為什么還要用MATLAB我在實際使用中的體會是很多實驗室和工程團隊的老課題代碼都是MATLAB寫的交互界面上也習(xí)慣用MATLAB做熱力圖和散點圖如果為了一個協(xié)變分析再搭一套Python環(huán)境先要解決BioPython版本、numpy版本這些依賴沖突可能比分析本身更耗時。AKtoolbox把從FASTA讀取到統(tǒng)計檢驗的完整流程都收進了一組前綴為ak_的函數(shù)。它不提供多序列比對功能輸入必須是已經(jīng)比對好的MSA這意味著它定位的是一個中段分析組件前面接任何比對軟件后面接MATLAB的可視化和統(tǒng)計分析。作為開源軟件另一個好處是當(dāng)你對某個結(jié)果存疑時可以直接打開ak_covar這類函數(shù)的源碼確認(rèn)里面到底有沒有做gap masking、有沒有去冗余而不是把算法當(dāng)黑匣子用。3. 環(huán)境準(zhǔn)備與數(shù)據(jù)格式把 FASTA 變成 MATLAB 喜歡的形狀3.1 拉源碼并登記路徑我手頭這版 AKtoolbox 是純命令行函數(shù)式架構(gòu)沒有 App Designer 界面目錄下主要是ak_前綴的腳本函數(shù)、一個測試目錄和幾個示例數(shù)據(jù)文件。沒有界面反而是好事命令行接口穩(wěn)定寫批處理也方便。拿到源碼后第一步不是馬上跑demo而是先把工具箱路徑加到 MATLAB 搜索路徑下并確認(rèn)核心函數(shù)能被識別。% 假設(shè)源碼解壓后放在 ~/src/AKtoolbox addpath(genpath(~/src/AKtoolbox)); % 檢查核心函數(shù)是否可見 which ak_readfasta這段代碼里genpath的作用是遞歸添加 AKtoolbox 下所有子目錄避免漏掉嵌套的工具函數(shù)目錄。執(zhí)行which ak_readfasta后如果返回的是完整路徑說明路徑配置正常如果返回“未找到”先檢查是否用了addpath(genpath(...))而不是只addpath外層目錄這是最常見的漏配原因。3.2 從 FASTA 到字符矩陣AKtoolbox 的輸入函數(shù)把 FASTA 文件讀進來后會保留每條序列的頭部信息和序列內(nèi)容再用另一個函數(shù)把變長序列轉(zhuǎn)成統(tǒng)一的字符矩陣。我不建議自己用fastaread之后直接堆字符串因為序列長度不一致時很容易在拼接處產(chǎn)生隱性錯誤。% 讀入已經(jīng)比對好的多序列比對文件 fastaData ak_readfasta(A3L1_aln.fasta); % 轉(zhuǎn)換成字符矩陣行序列列MSA位點 seqMat ak_alignment2matrix(fastaData); % 看一下規(guī)模序列數(shù)和位點數(shù) [ns, np] size(seqMat); fprintf(序列數(shù)%d位點數(shù)%d\n, ns, np);這里的fastaData在函數(shù)內(nèi)部被處理成一個結(jié)構(gòu)數(shù)組每個元素包含頭部和序列體。ak_alignment2matrix做的是字面量轉(zhuǎn)換把每一條序列展開成一行字符不足的長度用-補齊如果輸入的多序列比對里原本就有g(shù)ap這些gap會原樣保留。拿到seqMat之后后續(xù)所有函數(shù)都按這個矩陣接收數(shù)據(jù)這是整個分析鏈路的統(tǒng)一數(shù)據(jù)格式。3.3 關(guān)鍵預(yù)處理過濾高gap列并對近緣序列去冗余多序列比對不是拿來就能直接算協(xié)同進化的。很多MSA里存在這樣的列只有少數(shù)序列有殘基大部分序列在這個位置是-這種列帶來的“變異信號”根本不是進化約束而是比對軟件對插缺事件的占位。另外如果數(shù)據(jù)集里包含大量同源度極高的近緣序列比如同一菌株重復(fù)測序得到的序列它們會把譜系偏好放大讓少數(shù)遠(yuǎn)緣序列的變異信息被淹沒。常見的做法是先做一次位點過濾和去冗余。AKtoolbox里對應(yīng)函數(shù)是ak_filter_alignment% 過濾高gap列并對相似度高于90%的序列去冗余 [alnClean, info] ak_filter_alignment(seqMat, ... max_gap_frac, 0.3, ... min_identity, 0.9); % info里保存了保留下來的是原始MSA哪些列 origPos info.positions; fprintf(過濾后保留位點數(shù)%d\n, size(alnClean, 2));max_gap_frac表示一列中g(shù)ap占比超過30%就刪除min_identity表示序列相似度高于90%時只保留其中一條代表性序列。我自己的習(xí)慣是如果MSA本來序列數(shù)就不多min_identity可以放到0.95甚至1.0優(yōu)先保證樣本量如果序列冗余嚴(yán)重則改成0.8。這里還有一個必須養(yǎng)成的習(xí)慣info.positions一定要保存下來。一旦過濾函數(shù)執(zhí)行完列號就變成了“過濾后新矩陣的列號”而原始MSA的列號存在info.positions里后面做結(jié)構(gòu)映射時全靠它。3.4 數(shù)據(jù)規(guī)模不同時的參數(shù)決策預(yù)處理沒有一組參數(shù)能適配所有數(shù)據(jù)按照我跑過的幾個MSA案例粗估可以按下面的思路選數(shù)據(jù)情況建議參數(shù)理由序列數(shù)很多超過200條min_identity設(shè)為0.8~0.9去掉重復(fù)代表避免近緣序列主導(dǎo)序列數(shù)偏少只有30~50條min_identity設(shè)為0.95以上或不去冗余盡量保留樣本量否則統(tǒng)計不穩(wěn)定比對覆蓋很差gap比例高max_gap_frac設(shè)為0.3先砍掉大量無信息的占位列比對覆蓋很好gap很少max_gap_frac設(shè)為0.5只刪除那些幾乎全是gap的區(qū)域這里不用把這些值當(dāng)成鐵律它們解決的是同一個問題讓矩陣?yán)锏拿恳粋€統(tǒng)計量都盡量來自真實的殘基狀態(tài)而不是來自比對噪聲。4. 跑通一次默認(rèn)流程拿到互信息得分矩陣和顯著性估計4.1 最小運行命令數(shù)據(jù)預(yù)處理完之后就可以直接調(diào)用協(xié)同進化分析入口。AKtoolbox 里最常見的入口函數(shù)是ak_covar你不需要分別調(diào)用十幾個子函數(shù)只要把字符矩陣放進去并指定算法即可。res ak_covar(alnClean, ... method, mi, ... normalize, zscore, ... nperm, 1000);method參數(shù)指定算法mi表示互信息也可以換成sca或perturb。normalize指定標(biāo)準(zhǔn)化方式zscore會基于置換檢驗的隨機分布把原始分?jǐn)?shù)轉(zhuǎn)成z分?jǐn)?shù)。nperm是置換次數(shù)這里設(shè)置1000次意味著對每個位點對都做1000次隨機打亂的統(tǒng)計量估計用于后面計算經(jīng)驗p值。這里的邏輯是先計算觀測數(shù)據(jù)的原始互信息矩陣然后把每一列內(nèi)的殘基字符隨機打亂破壞位點之間的真實聯(lián)動再重新計算一次隨機情況下的互信息。重復(fù)nperm次之后就能得到隨機背景的均值和標(biāo)準(zhǔn)差再用觀測值減去隨機均值、除以標(biāo)準(zhǔn)差得到標(biāo)準(zhǔn)化后的得分。這個得分比原始互信息更有可比性因為它已經(jīng)減掉了背景偏差。4.2 返回結(jié)構(gòu)里有什么ak_covar返回的res是一個 MATLAB struct里面主要字段如下字段含義res.score標(biāo)準(zhǔn)化后的位點對得分矩陣對稱矩陣行列數(shù)等于輸入位點數(shù)res.zscore對應(yīng)的z分?jǐn)?shù)矩陣用于評估偏離隨機背景的程度res.pvalue置換檢驗得到的經(jīng)驗p值矩陣res.positions這些位點對應(yīng)原始MSA中的列號等價于之前info.positionsres.method本次運行使用的算法名第一次跑完建議先看size(res.score)確認(rèn)矩陣維度等于預(yù)處理后的列數(shù)而不是原始MSA列數(shù)。如果維度對不上多半前面某一步?jīng)]有把矩陣傳對。4.3 快速看最強的位點對得分矩陣是對稱矩陣第i行第j列和j行i列是同一個配對因此只看上三角就夠。下面這段代碼可以打印出得分最高的前10個位點對score res.score; % 只取上三角排除對角線 [I, J, V] find(triu(score, 1)); % 按得分降序排列 [Vs, idx] sort(V, descend); % 打印前10個候選 for k 1:min(10, numel(idx)) rowIdx I(idx(k)); colIdx J(idx(k)); fprintf(候選配對 bit %d - %d : score %.3f z %.2f p %.3f\n, ... rowIdx, colIdx, Vs(k), ... res.zscore(rowIdx, colIdx), ... res.pvalue(rowIdx, colIdx)); end這段代碼里的triu(score, 1)會保留矩陣對角線以上元素第二個參數(shù)1表示從對角線右側(cè)開始取直接去掉對角線。find返回的三個輸出分別是行列下標(biāo)和非零值。排序后打印時行的索引和列的索引要一起從I(idx(k))和J(idx(k))取避免排序后錯位?？吹降慕Y(jié)果里如果某個配對的p值特別低但z分?jǐn)?shù)卻不高要留意是否來自低頻殘基組合后面進階部分會提到怎么處理。4.4 參數(shù)組合怎么調(diào)實際分析很難一次到位的下面是我通常使用的參數(shù)組合參數(shù)推薦值說明method先mi后scaMI做快速初篩SCA做模塊驗證normalizezscore減掉隨機背景比原始分?jǐn)?shù)可解釋nperm1000追求統(tǒng)計穩(wěn)定性時使用nperm100調(diào)參階段先用快速模式有一點要特別注意切換method后不同算法的得分量綱完全不一樣SCA的權(quán)重甚至?xí)霈F(xiàn)負(fù)值不能拿MI的閾值去直接套SCA結(jié)果。調(diào)參階段我會先把nperm降到100甚至0這里nperm0表示只算原始分?jǐn)?shù)不做置換讓流程快速跑通確認(rèn)矩陣尺寸和輸出字段沒問題后再把置換次數(shù)提上去。5. 避坑與常見問題讓我返工的五個具體失誤5.1 沒有去冗余結(jié)果看起來“太完美”現(xiàn)象跑完IM后排名靠前的位點對幾乎都落在一個非常保守的結(jié)構(gòu)域上并且這些位點本身單獨看也很保守和前幾個位點的互信息高得像假數(shù)據(jù)一樣。原因MSA里如果包含大量相似度超過95%的近似序列它們在譜系上是重復(fù)抽樣會讓少數(shù)幾個常見殘基組合反復(fù)出現(xiàn)互信息被這些重復(fù)樣本刷高最終主導(dǎo)全局排名。協(xié)同進化分析真正需要的是多樣化的獨立序列而不是同一譜系的一堆副本。解決回到ak_filter_alignment把min_identity從0.9收緊到0.8必要時按構(gòu)建好的進化譜系進行加權(quán)抽樣。我一般會先打印去冗余前的序列數(shù)size(seqMat,1)和去冗余后size(alnClean,1)如果發(fā)現(xiàn)只去掉了不到10%的序列說明數(shù)據(jù)本身多樣性不足即使算出來結(jié)果也最好不要直接下強結(jié)論。5.2 把 gap 當(dāng)成一個普通字符現(xiàn)象得分矩陣?yán)镉幸徽麎K連續(xù)列在結(jié)果中排名特別高比如位點120到150。單獨看它們的原始互信息發(fā)現(xiàn)這些列在兩兩之間幾乎都是最高分但和功能沒有任何關(guān)系。原因這些列在MSA里是一大片gap區(qū)許多序列在這個區(qū)域整體缺失gap字符因此成對共現(xiàn)。在互信息計算眼里這等價于“有大量序列在同一個位置缺失”于是統(tǒng)計量表現(xiàn)出強關(guān)聯(lián)但這根本不是進化意義上的協(xié)同變異。解決用max_gap_frac先把高gap列過濾掉。如果某個區(qū)域沒有完全過濾干凈可以額外保留一個gap掩碼矩陣在調(diào)用ak_covar時把掩碼傳入至少不要讓-參與最后的字符頻次統(tǒng)計。我常用的做法是先把max_gap_frac定在0.3再人工看info.positions是否還有連續(xù)成片被保留的gap區(qū)域。5.3 過濾后的位點編號和結(jié)構(gòu)編號對不上現(xiàn)象結(jié)果顯示第42位和第88位是顯著共進化殘基對但把這兩個編號直接對應(yīng)到三維結(jié)構(gòu)時發(fā)現(xiàn)結(jié)構(gòu)里該區(qū)域根本不存在這兩個殘基或者對應(yīng)到了完全不同位置的殘基上。原因MSA本身不是從結(jié)構(gòu)編號開始的。原始比對里可能刪掉了一些空位而ak_filter_alignment又基于舊的MSA矩陣刪除了高gap列所有后續(xù)得到的編號都是“當(dāng)前矩陣的第幾列”而不是原始序列上的殘基編號更不是PDB文件里的編號。解決從一開始就保存info.positions每個分析結(jié)果里的位點編號都要通過它做一次映射把“當(dāng)前矩陣列號”轉(zhuǎn)成“原始MSA列號”再結(jié)合你讀入結(jié)構(gòu)文件時建立的序列到結(jié)構(gòu)編號映射做第二次轉(zhuǎn)換。這個映射我每次都會寫進腳本第一段而不是等結(jié)果出來后再手工拼。5.4 置換檢驗時間爆炸現(xiàn)象當(dāng)位點數(shù)到200以上、nperm設(shè)為1000時作業(yè)跑了幾個小時還沒結(jié)束。中間殺掉進程后檢查發(fā)現(xiàn)每個位點對都要重新掃一遍矩陣計算量直接膨脹。原因nperm是對每一位點對做置換重算位點數(shù)從100增加到200時兩兩配對數(shù)增加到4倍如果再疊加1000次置換總體耗時就是非線性上升。AKtoolbox默認(rèn)在置換時會重建整棵統(tǒng)計量矩陣保證p值的一致性但代價就是慢。解決調(diào)參階段先用nperm100甚至nperm0。如果最終報告里只關(guān)心排名前幾十的候選位點對可以只對候選對單獨做置換檢驗這樣計算量減少一個量級。另一個做法是使用塊置換而不是逐字符置換保持序列內(nèi)部相關(guān)性隨機性更貼近真實進化過程。5.5 小樣本硬套 SCA還拿譜去下結(jié)論現(xiàn)象只有20條序列跑SCA后得到的第一特征向量看起來有明確分層第一特征值解釋了很大方差但換用MI方法后這些“強信號位點”大多對不上。原因SCA在構(gòu)造協(xié)方差矩陣時需要足夠的獨立樣本去估計殘基之間的協(xié)變模式20條序列對L×L維矩陣來說自由度嚴(yán)重不足。少數(shù)幾條每序列自帶的特殊突變組合會被當(dāng)成穩(wěn)定特征放大看起來像模塊化信號實際上是抽樣噪聲。解決少于50條序列的MSA不要做SCA先用MI做初篩序列數(shù)足夠時再跑SCA并且對比兩種方法得到的強候選列表是否一致。如果兩個方法都指向同一組位點對結(jié)論才相對可信。小樣本下看到SCA譜再有“規(guī)律”也不要直接寫進報告除非用置換檢驗驗證過。6. 進階從得分矩陣到殘基配對網(wǎng)絡(luò)以及一個收尾習(xí)慣6.1 給得分矩陣做保守性屏蔽協(xié)同進化分?jǐn)?shù)高的位點對不一定都有生物學(xué)意義尤其是互信息方法它會把保守位點之間的微弱共變放得很大。一個實用的技巧是先計算每個位點的保守程度把過于保守的位點排除掉只保留有一定變異度、又表現(xiàn)出聯(lián)動信號的位點。% 計算每個位點的香農(nóng)熵-1表示完全保守0表示完全均勻 H ak_entropy(alnClean); conserv 1 - H; % 設(shè)置一個保守性閾值太保守的位點不進入最終網(wǎng)絡(luò) keep conserv 0.2; scoreClean res.score .* (keep keep);這里的ak_entropy返回每個MSA位點的歸一化香農(nóng)熵值conserv越接近1代表位點越保守。keep是一個邏輯向量keep keep會把它擴展成一個邏輯矩陣只有兩個位點都足夠可變時對應(yīng)位置上的得分才保留下來其他位置全部置零。這個操作能明顯減少保守位點主導(dǎo)的問題讓網(wǎng)絡(luò)更聚焦在有真實進化約束的殘基上。6.2 從矩陣到邊列表拿到清理后的得分矩陣下一步是把矩陣轉(zhuǎn)成網(wǎng)絡(luò)分析工具能讀的邊列表。我一般輸出一個三列TSV文件分別是位點A、位點B、得分。% 取上三角轉(zhuǎn)成網(wǎng)絡(luò)邊表 thr prctile(scoreClean(scoreClean eps), 99); [I, J, W] find(scoreClean thr); % 寫出邊表三列分別表示起始位點、目標(biāo)位點、權(quán)重 edges [I, J, W]; writematrix(edges, covar_edges.tsv, Delimiter, tab);閾值我通常設(shè)為得分矩陣非零部分的高分位點第99百分位數(shù)優(yōu)先保留最強信號。scoreClean thr會生成一個邏輯稀疏矩陣find提取出的三列就是邊表。你可以在網(wǎng)絡(luò)可視化軟件里用這個邊表畫殘基相互作用網(wǎng)絡(luò)也可以把它和三維結(jié)構(gòu)中的殘基距離做對比看那些顯著共進的殘基對是否在空間上確實接近。最后說一個讓我改變工作習(xí)慣的教訓(xùn)以前我拿到新MSA就直接跑協(xié)同進化結(jié)果總要先浪費一兩輪在檢查數(shù)據(jù)質(zhì)量和編號映射上。從那以后我每次跑AKtoolbox之前都強制自己走一遍“序列數(shù)、冗余度、gap比例、編號映射”四項檢查確認(rèn)無誤再分析。這個習(xí)慣幫我省掉了很多返工時間希望也能幫到你。本文還有配套的精品資源點擊獲取