「比较基因组学」简单个性化之画KaKs图

个性化需求

鉴定得到的wgd gene pair，TBtools计算得到TBtools.kaks.tab.xls，绘制Ks为横坐标，Ka为纵坐标的点图，对各数据进行分别以斜率k=1,k=0.5,k=0.05进行分区，不同分区使用不同颜色。
我已经快四年没有自己画过图了，平时比较忙，一般都是改改写写好的脚本，快忘光了。

运行命令

1.去掉TBtools软件计算结果文件TBtools.kaks.tab.xls中$3，$4中有-0.0、NaN、Infinity的行，其中还有重复的title，也一并去掉。

$awk '$3!~/-0.0/ &&  $3!~/NaN/ {print $0}' TBtools.kaks.tab.xls | awk '$4!~/-0.0/ &&  $4!~/NaN/  && $4!~/Infinity/ {print $0}' | awk '!visited[$0]++'  >TBtools.kaks.tab.final.xls

2.linux 命令行运行

$/share/nas2/genome/biosoft/R/3.6.1/bin/Rscript KaKs.R -i TBtools.kaks.tab.final.xls -o KaKs
Warning message:
Removed 9 rows containing missing values (geom_point). 
null device 
          1 
Warning message:
Removed 9 rows containing missing values (geom_point). 
null device 
          1 

3.脚本封装，KaKs.R 脚本内容如下：

#!/share/nas2/genome/biosoft/R/3.5.2/bin/Rscript
# _*_ coding: utf-8 _*_
library(getopt)
library(ggplot2)

 
command=matrix(c( 
  'help', 'h', 0,'loical', '帮助文档',
  'input', 'i', 2, 'character', '处理后的TBtools计算的KaKs文件',
  'output', 'o', 1, 'character', '结果文件'),byrow=T,ncol=5)
args=getopt(command)


if (!is.null(args$help) || is.null(args$input) || is.null(args$output) ) {
  cat(paste(getopt(command, usage = T), "\n"))
  q(status=1)
}
 

#print(args$input)
infile <- args$input


data <- read.table(infile,sep = "\t",header=T)
data$Total <- NULL
df1<-data[,c(3:4)]
color=c("red","green","blue")
df_abline <- data.frame(intercept=c(0,0,0),slope=c(1, .5,.05),linetype=factor(c(1,3,5)))
plot=ggplot(df1,aes(x=Ks,y=Ka))+geom_point()+geom_abline(data=df_abline, aes(intercept=intercept, slope=slope, color=color))+scale_colour_discrete(name="Ka/Ks",labels = c("k=0.05","k=0.5","k=1")) +xlim(0,5) +ylim(0,5)+theme(axis.title.x=element_text(size=20), axis.title.y=element_text(size=20), axis.text.x=element_text(size=10), axis.text.y=element_text(size=10),legend.text = element_text(size=10),plot.margin=unit(c(1, 1, 1, 1),'cm'),legend.position = "right"，legend.key = element_blank(),  panel.background = element_blank(), panel.border = element_rect(color="black",fill = "transparent"))
pdf(paste(args$output, 'pdf', sep = '.'), width = 8, height = 8)
plot
dev.off()
png(paste(args$output, 'png', sep = '.'), width = 4000, height = 4000, res = 600, units = 'px')
plot
dev.off()

4.画图结果如下：

KaKs.png