所用數(shù)據(jù)為一個(gè)屬內(nèi)不同種不同群體的葉綠體基因組序列顾稀,數(shù)量為80條博肋。
發(fā)現(xiàn)用全長序列建樹的時(shí)候汪厨,不適合選用太多外類群赃春,否則ML法中會(huì)導(dǎo)致屬內(nèi)分枝的枝長特別短。原因應(yīng)該是基因間隔區(qū)和內(nèi)含子區(qū)域序列位點(diǎn)的差異較大劫乱。
枝長含義
NJ:表示遺傳距離织中;
MP:性狀狀態(tài)變換的替換數(shù);
ML/BI:該分枝上的相對進(jìn)化數(shù)量(遺傳變異量)衷戈;每個(gè)位點(diǎn)上的替換數(shù)(一般以每位點(diǎn)多少次核苷酸替換或氨基酸取代來表示)狭吼。
遺傳距離
大多數(shù)情況以序列來說遺傳距離就是兩個(gè)OTU(個(gè)體、群體殖妇、物種或基因家族)之間序列的差異值刁笙。
序列比對
多序列比對用mafft得到的結(jié)果較為準(zhǔn)確,muscle比對的速度較快。
多序列比對的絕大多數(shù)算法都是基于漸進(jìn)比對的概念疲吸。簡單來說就是先從兩個(gè)序列的比對開始座每,逐漸添加新序列,直到所有的序列都加入為止摘悴。但是不同的添加順序會(huì)產(chǎn)生不同的比對結(jié)果峭梳。所以由最相似的兩個(gè)序列開始比對,由近到遠(yuǎn)逐步完成最為可靠蹂喻。
mafft --thread 15 --auto 80-AcoeOut.fasta > 80-AcoeOut_aln.fasta
##比對時(shí)如果不清楚什么參數(shù)合適葱椭,加個(gè)參數(shù) --auto,軟件可以自動(dòng)幫你處理
挑選保守位點(diǎn)進(jìn)行下一步建樹
序列比對完后口四,用于建樹的序列位點(diǎn)必須保證具有良好的同源性孵运。所以需要?jiǎng)h除序列分歧很大的區(qū)域和gap區(qū)域。
我用的軟件為Gblocks蔓彩,主要目的是把有g(shù)ap的位點(diǎn)全部去除治笨,參數(shù)為-b5=n
,其余的選項(xiàng)有-b5=h
粪小,h表示half 指去除在大于50%的序列中出現(xiàn)gap的位點(diǎn)。
Gblocks 80-AcoeOut_aln.fasta -t=d -b5=n
最大簡約法(軟件PAUP)
最大簡約法的樹長指所有性狀在一棵樹上的進(jìn)化改變總數(shù)抡句。
計(jì)算得到的結(jié)果可能會(huì)有許多樹長相等的簡約樹探膊,此時(shí)需要計(jì)算它們的一致樹。分為strict consensus和semistrict consensus等待榔,strict表示100%逞壁,在所有簡約樹中都出現(xiàn)的分枝,才會(huì)出現(xiàn)在一致樹中锐锣,否則為梳子腌闯。這個(gè)閾值可以調(diào)。
一般文章中所用的系統(tǒng)樹的拓?fù)浣Y(jié)構(gòu)都為ML或BI樹雕憔,所以要把MP的bootstrap值標(biāo)到ML/BI法的底樹上姿骏。
1.cstatus ##查看數(shù)據(jù)特征,如簡約信息位點(diǎn)數(shù)量等
2.tstatus
3.define outgroup
4.analysis-branch and bound ##數(shù)據(jù)量大用heuristic search
5.describetrees
6.contree ##計(jì)算一致樹
7.bootstrap ##自舉值
8.print bootstrap consensus
進(jìn)化模型
DNA序列進(jìn)化就是序列位點(diǎn)上的核苷酸隨時(shí)間的變化斤彼,主要包括堿基替換分瘦、缺失和插入。
兩條比對好的DNA序列的同源位點(diǎn)之間很容易看出堿基的相同或不同琉苇,但是在漫長的進(jìn)化過程中實(shí)際發(fā)生了什么我們并不知道嘲玫。最常見的當(dāng)然是單次替換,但是當(dāng)進(jìn)化時(shí)間較長時(shí)并扇,已經(jīng)發(fā)生過替換的某些位點(diǎn)可能會(huì)再次發(fā)生替換去团,即多重替換。
DNA序列的進(jìn)化模型將DNA的進(jìn)化作為一系列隨機(jī)突變來描述,并明確定義了4種堿基之間相互的替換速率土陪。
DNA進(jìn)化模型的參數(shù)主要有4類:
- 堿基組成頻率
- 替換速率矩陣:指定了4種堿基之間在單位時(shí)間內(nèi)或給定分枝長度下相互替換的概率昼汗。
- 不變位點(diǎn)比例
- 位點(diǎn)之間速率的變異:不同位點(diǎn)之間替換速率的差異。
ML法和BI法都需要選擇合適的進(jìn)化模型旺坠。模型選擇軟件具有的模型越多乔遮,檢測結(jié)果越準(zhǔn)確,但建樹軟件不一定支持該模型取刃。
判斷模型與數(shù)據(jù)擬合好壞的標(biāo)準(zhǔn)主要有AIC和BIC等蹋肮。
最大似然法(軟件IQ-TREE)
似然值是當(dāng)模型(樹和進(jìn)化參數(shù))為真時(shí)能夠得到實(shí)際觀測數(shù)據(jù)的概率。似然值是觀測數(shù)據(jù)(即序列)的條件概率璧疗,其條件為計(jì)算似然值時(shí)依據(jù)的模型坯辩,而不是模型為真時(shí)的概率。
ML法建樹的過程是先選擇一個(gè)適合數(shù)據(jù)集的進(jìn)化模型崩侠,然后對指定拓?fù)浣Y(jié)構(gòu)的一棵樹優(yōu)化分枝長度漆魔,以使得該拓?fù)浣Y(jié)構(gòu)的似然值最大化。通過計(jì)算不同拓?fù)浣Y(jié)構(gòu)樹的似然值却音,將具有最大似然值的樹看成是指定模型下的能夠產(chǎn)生觀測數(shù)據(jù)的最佳估計(jì)改抡。
ML法采用的搜索方法主要是啟發(fā)式搜索,步驟如下:
- 通過NJ樹或逐步添加序列的方法構(gòu)建初始樹系瓢;
- 以初始樹為基礎(chǔ)通過各種分枝交換方法(TBR阿纤、SPR等)計(jì)算似然值,將最大似然值的樹保存夷陋,并作為下一輪重排的初始樹欠拾;
- 重復(fù)進(jìn)行分枝交換,直到不能增加似然值為止骗绕。重排的最后獲得的最大似然值樹即為ML樹藐窄。
建ML樹的軟件用RAxML的較多,但近來IQ-TREE的引用量一路上升酬土。綜合使用下來荆忍,個(gè)人感覺IQ-TREE的速度真快。
使用過程是下載了PhyloSuite的組件撤缴,從選模型到構(gòu)樹一站式操作還挺方便的东揣。注意下載好后首先要配置用于不同分析的插件。
貝葉斯推論法(軟件MrBayes)
BI法與ML法不同的是腹泌,前者根據(jù)提供的數(shù)據(jù)和選擇的替代模型尋找可能性最大的樹嘶卧,而ML法則是尋找合適的樹以使得數(shù)據(jù)的可能性最大。
推斷系統(tǒng)發(fā)育樹的步驟為:
- 選擇一些樹作為起始點(diǎn)凉袱;
- 判定這些樹的似然值芥吟;
- 修改樹的拓?fù)浣Y(jié)構(gòu)和分支長度侦铜;
- 計(jì)算出新樹的似然值;
- 新樹的似然值比舊樹大钟鸵,則接受新樹钉稍。
如此就構(gòu)成了一代,一次又一次的重復(fù)迭代棺耍,直到新樹的似然值不再有明顯變化贡未,即樹的似然值不再有顯著區(qū)別,參數(shù)已收斂為止蒙袍。如果沒有收斂俊卤,適當(dāng)?shù)脑黾哟鷺淅^續(xù)跑。
如何判斷參數(shù)是否已收斂
軟件運(yùn)行完畢后害幅,看結(jié)果文件的分離頻率平均標(biāo)準(zhǔn)差值(Average standard deviation of split frequencies)
該值<0.01時(shí)消恍,說明兩次運(yùn)行的結(jié)果差異很少,參數(shù)已收斂以现;
該值>0.05時(shí)狠怨,需要繼續(xù)運(yùn)行。
同樣是在PhyloSuite中運(yùn)行
##只輸出偶數(shù)行
sed -n '2~2p' test.txt
##只輸出奇數(shù)行
sed -n '1~2p' test.txt
## 刪除空行
sed '/^$/d' test.txt
## 計(jì)算行數(shù)
sed -n '$=' test.txt