使用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)建的,以顯示哪個病原體樣本在遺傳上更接近另一個樣本,從而深入了解潛在的未觀察到的流行病學聯(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)形進化樹
?
##環(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多生信入門實操知識。
