论文
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.fablast
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"))
怎么添加上下的刻度，另外再写推文介绍吧
还有一个问题是 水稻这一段序列看起来和下面这个是反向的，那么如果水稻序列取反向互补，那么原来的基因位置坐标应该如何转换，这个暂时想不明白
推文记录的是自己的学习笔记，内容可能会存在错误，请大家批判着看，欢迎大家指出其中的错误示例数据和代码可以给推文点赞，然后点击在看，最后留言获取
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
