使用R語言完成序列比對及進化樹美化
爾云間? 一個專門做科研的團隊
???
系統(tǒng)發(fā)育樹可以描述一組生物之間的系譜關(guān)系，如人和黑猩猩的親緣關(guān)系最為接近，而與大猩猩則次之。
拿血紅蛋白舉例，人和黑猩猩是從大約700w年由同一祖先進化而來，它們的血紅蛋白只有一個氨基酸的差異，和大猩猩的血紅蛋白差異則較大。
香蕉。一種連動物都不是的生物，居然跟人類的基因有50%的相似性。
我們可以根據(jù)生物的基因序列構(gòu)建，也可以根據(jù)蛋白質(zhì)序列構(gòu)建。如下圖所示，一個有根的系統(tǒng)發(fā)育樹代表了一種進化歷史模型，該模型由樹節(jié)點之間的祖先-后代關(guān)系以及“sister”或“cousin”有機體在不同相關(guān)性水平上的聚類所描述。在傳染病研究中，系統(tǒng)發(fā)育樹通常是根據(jù)病原體的基因或基因組序列構(gòu)建的，以顯示哪個病原體樣本在遺傳上更接近另一個樣本，從而深入了解潛在的未觀察到的流行病學(xué)聯(lián)系和爆發(fā)的潛在來源。
下面小果帶大家使用R語言完成序列比對和進化樹的繪制：
序列下載：
打開NCBI ,選擇Nucleotide數(shù)據(jù)庫，輸入Coronavirus ，回車，選擇你需要的序列，選擇send to進行下載，
保存成fasta格式，即以>開頭的文件，
下面是R語言實操部分，小果只推薦跑10條左右的短序列，如果太多還是使用mega等軟件比較好：
###BiocManager::install("msa")
require(msa)
mySequenceFile <- readAAStringSet(“sequence.txt”)
#多序列比對
myFirstAlignment <- msa(mySequenceFile)
head(mySequenceFile)
library(ggplot2)
require(seqinr)
myAlignment <- msaConvert(myFirstAlignment, type="seqinr::alignment")
d <- dist.alignment(myAlignment, "identity")
#構(gòu)建NJ樹
require(ape)
tree <- nj(d)
#畫樹并輸出到PDF文件ggtree.pdf
require(ggtree)
- ##環(huán)狀圖
- p1<-ggtree(tree, layout='circular', ladderize=FALSE, size=0.8, branch.length="none",col="red")+
- ? geom_tiplab2(hjust=-0.3)+
- ? geom_tippoint(size=1.5,col="blue")+
- ? geom_nodepoint(color="black", alpha=1/4, size=2) +
- ? theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5)))
- # 圖例位置、文字大小
###長方形圖
p2<- ggtree(tree, layout='rectangular', size=0.8, col="deepskyblue3") +
?
? geom_tiplab(size=3, color="purple4", hjust=-0.05)+
? geom_tippoint(size=1.5, color="deepskyblue3")+
? geom_nodepoint(color="pink", alpha=1/4, size=5)+
? theme_tree2()
ggsave(p1, file="tree_circular.pdf", width=9, height=9)
ggsave(p2, file="tree_rectangular.pdf", width=9, height=9)
進化樹的美化大家可以使用itol或者自己閱讀ggtree包的參數(shù)了，好了今天小果的分享就到這里了，關(guān)注小果，下期會給大家?guī)斫o多生信入門實操知識。
推薦閱讀
