好久没更新了，最近太忙了，甚至双十一买的东西现在还没发货我都给忘了。（题外话）
接下来要给大家分享的是拿到一段区域的区间后如何把这段区域的基因展现出来。 代码略为粗糙，如果哪位朋友有更合适的方法或者改进方法，欢迎留言。
- 1.首先把候选区域的正反义链的基因文件提取出来，如下：
echo "startpos endpos ano"|tr " " "\t">positive-strand.chr4.gwas.txt #为了便于R中merge函数合并
cat test.gtf|awk '{if($1==4&&$3=="gene"&&$5>=69597209&&$4<=73300862&&$7=="+"){print$0}}'|cut -f4,5,9>>positive-strand.chr4.gwas.txt #正义链
echo "startmin endmin ano"|tr " " "\t">minus-strand.chr4.gwas.txt
cat test.gtf|awk '{if($1==4&&$3=="gene"&&$5>=69597209&&$4<=73300862&&$7=="-"){print$0}}'|cut -f4,5,9>>minus-strand.chr4.gwas.txt #负义链
- 2.用R语言处理文件并画图
library(ggplot2)
data1=read.table(file="positive-strand.chr4.gwas.txt",header = T,sep = "\t")
data2=read.table(file="minus-strand.chr4.gwas.txt",header = T,sep = "\t")
data=merge(data1,data2,by="ano",all=T) #正反义链互补的部分都会变成NA
ggplot(data)+scale_x_continuous(breaks =seq(70000,73000,500),limits=c(69597.209,73300.862),expand=(c(0,0)))+scale_y_continuous(breaks =seq(0,5,2),limits=c(0,5))+theme_bw(base_size = 16)+theme(panel.grid = element_blank(),axis.line = element_line(colour = 'black'),panel.background = element_rect(fill = "transparent"),axis.text.y =element_blank(),axis.ticks.y=element_blank(),axis.text.x = element_text(size=9.5,face="bold", hjust = 0.5, vjust = 0.5))+labs(x ='Chromsome 4 Position (kb)', y = expression("")) +geom_rect(aes(xmin=data$startpos/1000, xmax=data$endpos/1000, ymin=1, ymax=2),fill='red',alpha = 0.8)+geom_rect(aes(xmin=data$startmin/1000, xmax=data$endmin/1000, ymin=3, ymax=4),fill='#93e3e5',alpha = 0.8)
ggsave("genestructure.tiff",width=12,height=2.5)
图像如下
image.png
