以前做序列比对通常都是MEGA来做，然后有时候会用genedoc来展示，但是经常碰到一个问题就是没办法得到矢量图。所以一直在想用R能不能展示多序列比对的结果，尤其是核心关键domain的区域。
下面是paper中的一个结果，看着就很清晰。
library("phylotools")
library("ggplot2")
library("tidyr")
library("dplyr")
df <- read.fasta("data.fasta")
df
#用ggplot展示的话，我们就需要拆分成单个碱基的形式，有点类似长矩阵变成短矩阵的感觉。
new.df <- df %>% separate(seq.text,paste0("col",str_pad(1:28,2,side = "left",pad = "0")),'') %>% select(-col01) %>%mutate(seq.name=factor(seq.name,levels = rev(seq.name))) %>% pivot_longer(!seq.name)
下面就可以用ggplot画图了
p1<-ggplot(new.df,aes(x=name,y=seq.name))+
geom_tile(fill="white",alpha=0)+
geom_text(aes(label=value))+
theme_bw()+
theme(axis.text.x = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
panel.border = element_blank())
p1
下面添加两侧的数字：
df1<-data.frame(x=0,
y=11:1,
label=c(671,668,669,666,744,743,756,736,706,747,759)) #左侧
df2<-data.frame(x=28,
y=11:1,
label=c(695,692,693,690,770,665,777,761,727,768,780)) #右侧
p2 <- p1+
geom_text(data=df1,
aes(x=x,y=y,label=label),
inherit.aes = FALSE,hjust=1)+
geom_text(data=df2,
aes(x=x,y=y,label=label),
inherit.aes = FALSE,hjust=0)+
coord_equal(xlim = c(-1,29))
p2
下面就是添加颜色了：
我们把控制颜色的数据单独放在一个文件中，我只准备了部分。
df3 <- read.table("tmp.txt",sep="\t",header=T)
df3
p2 +
geom_tile(data=df3,aes(x=x,y=y,fill=group),
color="white",show.legend = FALSE)+
geom_text(data=df3,aes(x=x,y=y,label=label),
color="white")+
scale_fill_manual(values = c("#00adef","#ed1b24"))+
labs(x=NULL,y=NULL)
这样经过控制，就可以获得文中的效果了。
网友评论