跟着Science学作图：R语言ggplot2作图展示基因组局部区域的共线性

论文

Convergent selection of a WD40 protein that enhances grain yield in maize and rice

science.abg7985.pdf

https://www.science.org/doi/10.1126/science.abg7985

今天的推文我们来试着复现一下论文中的Fig3b

jcvi这个工具可以做这个基因组局部的共线性，jcvi的链接

https://github.com/tanghaibao/jcvi

如果有数据用ggplot2来做可能可定制性会高一些

准备数据

每个区间的bed文件

水稻

Chr4 28500000 28600000

玉米

2 17650000 18050000

然后用bed文件和对应的gff文件取交集提取区间内的基因

bedtools intersect -a ../../maize/PhytozomeV13/Zmays/RefGen_V4/annotation/Zmays_493_RefGen_V4.gene.gff3 -b maize.bed | awk '$3=="gene" {print}' > maize.gene.gff

bedtools intersect -a ../../rice/Phytozome/PhytozomeV12/Osativa/annotation/Osativa_323_v7.0.gene.gff3 -b rice.bed | awk '$3=="gene" {print}' > rice.gene.gff

这个区间内的共线性关系如果有现成的就好了，我这里的处理方式是把两段区间内的序列提取出来，然后做blast，然后用blast的结果作为共线性的关系（我这里仅仅是为了获得作图数据，不太确定这种方式作为共线性是否合理）

samtools faidx ../../maize/PhytozomeV13/Zmays/RefGen_V4/assembly/Zmays_493_APGv4.fa 2:17650000-18050000 > maize.fa

samtools faidx ../../rice/Phytozome/PhytozomeV12/Osativa/assembly/Osativa_323_v7.0.fa Chr4:28500000-28600000 > rice.fa

blast

makeblastdb -in rice.fa -dbtype nucl -out rice
blastn -query maize.fa -db rice -outfmt 6 > rice.maize.blastn

作图

library(tidyverse)
library(gggenes)
library(ggforce)

rice<-read_delim("D:/R_4_1_0_working_directory/env001/data/20230503/rice.gene.gff",
                 delim = "t",
                 col_names = FALSE) %>% 
  mutate(X9=str_replace(str_sub(X9,1,17),"ID=",""),
         X7=case_when(
           X7 == "-" ~ -1,
           TRUE ~ 1
         ),
         X10=case_when(
           X9 == "LOC_Os04g47990" ~ "A",
           TRUE ~ "B"
         ),
         X1=2) %>% 
  select(X4,X5,X7,X9,X10,X1)

rice

maize<-read_delim("D:/R_4_1_0_working_directory/env001/data/20230503/maize.gene.gff",
                 delim = "t",
                 col_names = FALSE) %>%  
  mutate(X9=str_replace(str_sub(X9,1,17),"ID=",""),
         X7=case_when(
           X7 == "-" ~ -1,
           TRUE ~ 1
         ),
         X10=case_when(
           X9 == "Zm00001d002641" ~ "A",
           TRUE ~ "B"
         ),
         X1=1) %>% 
  select(X4,X5,X7,X9,X10,X1)

maize

rice.maize.blastn<-read_delim("D:/R_4_1_0_working_directory/env001/data/20230503/rice.maize.blastn",
                              col_names = FALSE,
                              delim = "t") %>% 
  mutate(X7=X7+17650000-1+10700000,
         X8=X8+17650000-1+10700000,
         X9=X9+28500000-1,
         X10=X10+28500000-1) %>% 
  select(X7,X8,X9,X10)

rice.maize.blastn

myabc<-function(df){
  x1<-c()
  
  for(i in 1:nrow(df)){
    
    x1<-append(x1,c(t(df)[,i][c(3,1,2,4)] %>% as.vector()))
  }
  return(x1)
}

x1<-myabc(rice.maize.blastn)
x1
y1<-rep(c(2,1,1,2),nrow(rice.maize.blastn))
y1
group1<-rep(as.character(1:400)[1:nrow(rice.maize.blastn)],each=4)
group1

ggplot()+
  geom_diagonal_wide(data = data.frame(x=x1,y=y1,group=group1),
                     aes(x=x,y=y,group=group),color="grey",fill="grey")+
  annotate(geom="segment",x=28500000,xend=28600000,y=2,yend=2)+
  geom_gene_arrow(data=rice,
                  aes(xmin=X4,xmax=X5,y=X1,
                      forward=X7,fill=X10),
                  arrowhead_height = unit(3, "mm"), 
                  arrowhead_width = unit(1, "mm"),
                  show.legend = FALSE)+
  annotate(geom="segment",x=28500000-150000,xend=28600000+160000,y=1,yend=1)+
  geom_gene_arrow(data=maize,
                  aes(xmin=X4+10700000,xmax=X5+10700000,y=X1,
                      forward=X7,fill=X10),
                  arrowhead_height = unit(3, "mm"), 
                  arrowhead_width = unit(1, "mm"),
                  show.legend = FALSE)+
  theme_bw()+
  theme(panel.grid = element_blank(),
        panel.border = element_blank(),
        axis.title = element_blank(),
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank())+
  scale_fill_manual(values = c("#f0505a","#65b691"))

怎么添加上下的刻度，另外再写推文介绍吧

还有一个问题是 水稻这一段序列看起来和下面这个是反向的，那么如果水稻序列取反向互补，那么原来的基因位置坐标应该如何转换，这个暂时想不明白

推文记录的是自己的学习笔记，内容可能会存在错误，请大家批判着看，欢迎大家指出其中的错误

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