跟着PNAS学作图：R语言ggplot2作图展示多序列比对结果

论文

The evidence remains clear: SARS-CoV-2 emerged via the wildlife trade

https://www.pnas.org/doi/10.1073/pnas.2214427119

在 饶毅科学 公众号 看到的推文 美国科学院院刊：逐条反驳新冠是实验室产物的错误说法

其中有一个图是

今天的推文我们来试着复现一下这个图

首先是一个多序列比对文件

读取数据

df <- phylotools::read.fasta("data/20221126/pnas.fasta") 
df

把序列拆分成一个碱基一列

df %>% 
  tidyr::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) -> new.df

这里有一个问题是序列是27个碱基，拆分的时候需要指定28列，然后第一列是空的，暂时没有想明白是为啥

首先是多序列比对

p1<-new.df %>% 
  ggplot(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))
  
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
p2

在自己感兴趣的地方添加背景色

部分示例数据

df3<-readxl::read_excel("data/20221126/data.xlsx")
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)

示例数据和代码可以给推文点赞，点击在看，最后留言获取

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！