個性化需求
鑒定得到的wgd gene pair，TBtools計算得到TBtools.kaks.tab.xls，繪制Ks為橫坐標，Ka為縱坐標的點圖擎椰，對各數據進行分別以斜率k=1,k=0.5,k=0.05進行分區(qū)，不同分區(qū)使用不同顏色创肥。
我已經快四年沒有自己畫過圖了达舒，平時比較忙，一般都是改改寫寫好的腳本叹侄，快忘光了休弃。
運行命令
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
