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