|
寫這篇文章的原因： 前段時間注意到一篇給CorMut提意見的文章，一直沒有時間回應，終于在假期用幾天時間把CorMut的代碼逐行研讀，現在做一個簡單的說明。所提問題主要為三個：第一，序列處理時刪除gap造成密碼子錯位。CorMut依照序列文件第一條的參考序列進行去除gap，不是去除所有的gap，去除的也是跟參考序列對應位置的gap，不可能造成密碼子錯位的發(fā)生。 第二，考慮到序列中可能存在混合堿基，所以有一個去除混合堿基的步驟，這是按照核苷酸為單位進行，遍歷到哪個核苷酸返回那個位置涉及的密碼子然后與相應參考序列位置的密碼子比較，看是否引起氨基酸突變，最終選擇引起突變的核苷酸替換混合堿基。數數也是沒有問題的。 第三，計算LOD值時會用到二項分布，原版本是直接用分布的公式進行計算，這樣在數值較大時的確會造成誤差（大樣本時choose函數會產生Inf等異常值），所以使用二項分布的dbinom會避免大樣本時異常值的產生，但是1000以下的小樣本對結果是沒有影響的，不可能產生NaN。這的確是當時沒有考慮到的問題。在新版本中已經進行了更新。 經過對代碼逐行檢查，并沒有發(fā)現其他原則性的問題。同時在這個版本對部分內容進行了優(yōu)化，對一些有分母的公式進行了一個平滑處理，避免了分母是0時造成異常值的情況；對畫圖功能增加了選擇網絡呈現形式的參數，用戶可以自由設置網絡以何種方式呈現。
就像其他軟件一樣，CorMut肯定時有bug的，過去我們也使用了大量數據進行測試，讓CorMut變的更加完美。CorMut始終是開放的，我們會積極的接受任務批評和建議，遇到問題的話，歡迎與我進行交流。
CorMut開發(fā)始末： 最早CorMut使用的主流模型是用在HIV耐藥位點識別上。當時并沒有現成的工具來解決這件事情，當時研究組正在開展HIV耐藥位點鑒定的工作，因此我們便首先根據這個算法編寫了初步的腳本，使用這個算法對自己的數據進行分析，并對耐藥位點實驗驗證，發(fā)現了幾個全新的耐藥位點。并對多套數據進行分析發(fā)現了與生物學相一致的結果。后來萌生了編寫一個軟件包的想法，然后就有了CorMut。原有的算法是ka/ks和條件ka/ks，單單這個肯定是不夠的，后來增加了互信息，杰卡德系數等算法，對CorMut進行了較大的擴充，并增加了繪制網絡圖的功能。2014年博士畢業(yè)之后，參加工作，研究方向發(fā)生了較大轉變，因為種種原因，CorMut也基本停止了更新。最近意識到CorMut不能丟掉，對CorMut的維護也是自我學習的一個過程。首先，在這里公布一個比較全面的使用指南，供參考。
CorMut使用指南： 相關突變的背景介紹： 在遺傳學中，Ka/Ks的比值是指每個非同義位點非同義取代的數目與每個同義位點同義取代數目的比值。Ka/Ks比值可以被用作度量作用于蛋白質編碼基因選擇壓力的指標。Ka/Ks等于1表示中性選擇，也就是說觀察到的非同義突變與同義突變的比值與期望的隨機突變模型的比值相匹配。因此，氨基酸改變既不被選擇也不被淘汰。Ka/Ks大于1表明氨基酸的改變是進化偏好的，也就是說這些突變能夠更強組織的適應性。這種不尋常的狀態(tài)可能反映了基因功能的改變或強迫機體去適應的環(huán)境條件改變。例如，高度變異的病毒，如HIV和HCV，對于新的抗病毒藥物的突變可能在接受這些藥物治療的感染者體內是受到陽性選擇的。相關突變（CorrelatedMutation）是進化生物學中的一個基本概念，在一個可編碼蛋白質的基因上，氨基酸突變會受到蛋白質功能的限制。基于蛋白質功能的限制，一個突變的發(fā)生可能引發(fā)另外一個補償突變的發(fā)生，因此，在特定條件下特定的突變模式就會形成，理解這些突變模式對于了解特定條件下基因的功能的改變至關重要。 CorMut組合了Ka/Ks比值和相關突變分析。軟件包中所使用的計算單個位點Ka/Ks的方法是由Chen等提出的。CorMut提供了計算單個位點或特異氨基酸Ka/Ks并發(fā)現他們之間相關突變的函數。CorMut提供了三種方法來發(fā)現相關突變，包括條件選擇壓力、互信息和Jaccard系數。計算主要有兩個步驟組成：第一，發(fā)現陽性選擇位點；第二，計算陽性選擇位點之間的突變關聯。 值得注意的是第一個步驟是可選擇的。同時，CorMut可以通過構建相關突變網絡方便的比較兩種條件下的相關突變關系。 CorMut原理介紹： 方法： 1. 單個密碼子和單個氨基酸取代Ka/Ks的計算 單個密碼子或單個特異的氨基酸取代的Ka/Ks是基于Chen的算法確定的。單個密碼子或單個特異的氨基酸取代(X2Y)的Ka/Ks由下列公式計算： 分子中，NY和NS分別表示某個氨基酸發(fā)生有義突變和無義突變的數目；然后使用分母中的隨機突變模型（也就是不存在選擇壓力情況下的模型）對該值進行校正，也就是式子的分母部分。在隨機突變模型中，ft和fv分別表示轉化和顛換的頻率。它們的計算公式為：ft=Nt/ntS和fv=Nv/nvS（其中S表示樣本總數；Nt和Nv表示所有樣本中RT區(qū)實際發(fā)生轉換和顛換的數量；nt和nv表示在所檢測的序列中可能發(fā)生轉換和顛換的數量（分別等于樣本的數量和兩倍的樣本數量）。分母中，nY,t和nSt分別表示序列中既發(fā)生顛換又發(fā)生無義突變的數目和既發(fā)生顛換又發(fā)生有義突變的數目，nY,v和nSv分別表示序列中既發(fā)生轉換又發(fā)生無義突變的數目和既發(fā)生轉換又發(fā)生有義突變的數目。同時，對突變使用對數比值比（logodds ratio，LOD）進行評估。LOD的計算公式如下：
N表示密碼子的突變總數；q計算公式如下： 如果LOD大于2，陽性選擇被認為是顯著的。 2. 條件Ka/Ks的計算 條件選擇比值可以定義為氨基酸Y突變時氨基酸Z突變的Ka/Ks與氨基酸Y沒有突變時氨基酸Z突變的Ka/Ks的比值，計算公式為： 公式中NZaYa表示同時存在氨基酸突變Y和Z的樣本總數； NZsYa表示Z同義突變， Y氨基酸突變的樣本總數；NZaYo和NZsYo分別表示在不發(fā)生氨基酸Y突變的情況下，氨基酸Z非同義突變和同義突變的樣本總數。 LOD計算公式如下： 3. 互信息 互信息(MutualInformation)是信息論里提到的一種的信息度量的方法，被用來度量兩個事件集合之間的相關性[37]。已有文獻報道互信息被用作度量殘基取代之間的相關性。也就是說，在一個具有N個密碼子的蛋白質的多序列多序列比對中N列中的每一列被認為是一個離散變量Xi(1<><>。當計算密碼子之間的互信息時，Xi會計算20種氨基酸的每一種的概率，而計算兩個密碼子特定氨基酸的互信息時，Xi只計算兩個部分的概率，也就特異的氨基酸取代或者其他取代?；バ畔⒚枋鰞蓚€變量Xi和Xj之間的依賴關系，通過關聯Xi和Xj兩列的相關指標來度量兩者依賴的強弱。
4. Jaccard 系數 Jaccard系數可以度量兩個變量之間的相似性，并且已經廣泛應用于度量相關突變[38-40]。對于一個突變對X和Y，Jaccard系數計算的公式為Nxy/(Nxy+Nx0+Ny0)，在式中，Nxy代表X和Y都發(fā)生突變的序列的條數，Nx0代表只有X突變而Y沒有突變的序列的數目，Ny0只有Y突變而X沒有突變的序列的數目。Jaccard系數被認為可以有效的避免X和Y都沒有發(fā)生突變序列條數較多引起的的稀有突變結果的夸大計算。
CorMut包含的主要函數及使用說明： ckaksCodon, ckaksAA, miCodon, miAA, jiAA設置了參數 'setPosition'，可以指定感興趣的位置進行計算，避免其他位點影響結果的展示，降低運算的復雜度。 CorMut實戰(zhàn)： 1. 多序列比對文件的處理 seqFormat函數能夠用普通堿基取代混合堿基，基于參考序列刪除序列比對中的空缺。由于一個混合堿基可能與多個普通堿基相對應，seqFormat會隨機選擇能夠引發(fā)氨基酸突變的普通堿基。這里，未治療和治療的HIVPR序列被用作為例子來演示函數的使用： library(CorMut) examplefile=system.file('extdata','PI_treatment.aln',package='CorMut') examplefile02=system.file('extdata','PI_treatment_naive.aln',package='CorMut') example=seqFormat(examplefile) example02=seqFormat(examplefile02) 2. 計算單個密碼子位置的kaks result=kaksCodon(example) fresult=filterSites(result) head(fresult)
3. 計算單個氨基酸突變的Ka/Ks result=kaksAA(example) fresult=filterSites(result) head(fresult) 4. 計算密碼子之間的條件Ka/Ks（條件選擇壓力） result=ckaksCodon(example) plot(result) fresult=filterSites(result) head(fresult) 5. 計算氨基酸突變之間的條件Ka/Ks（條件選擇壓力） result=ckaksAA(example) plot(result) fresult=filterSites(result) head(fresult)
6. 計算密碼子之間的信息熵 包含兩個矩陣的列表將會被返回，他們分別是互信息矩陣和p值矩陣。矩陣的第i行第j列元素代表第i個密碼子和第j密碼子之間互信息的值或p值。 result=miCodon(example) plot(result) fresult=filterSites(result) head(fresult) 7. 計算氨基酸突變之間的信息熵 包含兩個矩陣的列表將會被返回，他們分別是互信息矩陣和p值矩陣。矩陣的第i行第j列元素代表第i個氨基酸突變和第j個氨基酸突變之間互信息的值或p值。 result=miAA(example) plot(result) fresult=filterSites(result) head(fresult)
8. 計算氨基酸突變之間的Jaccard系數 包含兩個矩陣的列表將會被返回，他們分別是Jaccard系數矩陣和p值矩陣。矩陣的第i行第j列元素代表第i個氨基酸突變和第j氨基酸突變之間Jaccard系數的值或p值。 result=jiAA(example) plot(result) fresult=filterSites(result) head(fresult) 9. 通過構建相關突變網絡比較兩種條件下的相關突變 biCompare函數能夠比較兩種條件下的相關突變關系，比如HIV未治療和治療情況下相關突變關系的改變. biCompare 的結果能夠通過plot函數實現可視化。只有陽性選擇密碼子或氨基酸突變被展現在圖上，藍色圓圈代表第一種條件下獨特的陽性選擇密碼子或氨基酸，也就是說這些位點在第二種條件下是受到中性選擇或陰性選擇.而紅色圓圈代表第二種條件下出現的獨特的陽性選擇的密碼子或氨基酸?？梢酝ㄟ^控制Plot的參數plotUnchanged展現那些兩種條件下均為陽性選擇的位點，plotUnchanged為FALSE時可以實現這一點。 舉例: biexample=biCkaksAA(example02,example) #biexample=biCkaksCodon (example02,example) #biexample=biMICodon (example02,example) #biexample=biMIAA (example02,example) #biexample=biJIAA (example02,example) result=biCompare(biexample) plot(result)
|
|
|
