散点图归类作图脚本:
图例:
代码:
library(ggplot2)
b <- read.table("WB.RNA_PEP_Quadrant.txt",header = T)
b <- as.data.frame(b)
cva <- cor(x = b$log2FC_RNA,y = b$log2FC_PEP,method = "pearson")
r <- cor.test(x = b$log2FC_RNA, y = b$log2FC_PEP)
pva <- r$p.value
title <-paste("Pearson's correlation:",cva,"\np-value:",pva)
##
pdf(file = "WB.RNA_PEP_Quadrant.pdf",width = 10,height = 6.18)
p <- ggplot(data = b,aes(b$log2FC_PEP, b$log2FC_RNA)) + theme_minimal()
p <- p + geom_point(aes(colour=factor(b$Quadrant)),show.legend = FALSE)
p <- p + geom_hline(aes(yintercept = -1),linetype="dashed") + geom_hline(aes(yintercept = 1),linetype="dashed")
p <- p + geom_vline(aes(xintercept = -0.66),linetype="dashed") + geom_vline(aes(xintercept = 0.66),linetype="dashed")
p <- p + xlab(label = "log2(ratio of protein)") + ylab(label = "log2(ratio of transcript)")
p <- p + scale_x_continuous(breaks = seq(-3,3,1)) + scale_y_continuous(breaks = seq(-10,10,1))
p <- p + theme(panel.border = element_rect(colour = "black",fill = NA),panel.background = element_rect(fill = NA))
p <- p + ggtitle(label = title) + theme(plot.title = element_text(hjust = 0.5))
p
dev.off()
数据源:
写成模块:
#!/usr/bin/perl -w
use strict;
my $name = $ARGV[0] || die $!;
open RCODE, ">$name.R" || die $!;
print RCODE<<EOF;
library(ggplot2)
b <- read.table("$name.RNA_PEP_Quadrant.txt",header = T)
b <- as.data.frame(b)
cva <- cor(x = b\$log2FC_RNA,y = b\$log2FC_PEP,method = "pearson")
r <- cor.test(x = b\$log2FC_RNA, y = b\$log2FC_PEP)
pva <- r\$p.value
title <-paste("Pearson's correlation:",cva,"\\np-value:",pva)
##
pdf(file = "$name.RNA_PEP_Quadrant.pdf",width = 10,height = 6.18)
p <- ggplot(data = b,aes(b\$log2FC_PEP, b\$log2FC_RNA)) + theme_minimal()
p <- p + geom_point(aes(colour=factor(b\$Quadrant)),show.legend = FALSE)
p <- p + geom_hline(aes(yintercept = -1),linetype="dashed") + geom_hline(aes(yintercept = 1),linetype="dashed")
p <- p + geom_vline(aes(xintercept = -0.66),linetype="dashed") + geom_vline(aes(xintercept = 0.66),linetype="dashed")
p <- p + xlab(label = "log2(ratio of protein)") + ylab(label = "log2(ratio of transcript)")
p <- p + scale_x_continuous(breaks = seq(-3,3,1)) + scale_y_continuous(breaks = seq(-10,10,1))
p <- p + theme(panel.border = element_rect(colour = "black",fill = NA),panel.background = element_rect(fill = NA))
p <- p + ggtitle(label = title) + theme(plot.title = element_text(hjust = 0.5))
p
dev.off()
EOF
close RCODE;
system("Rscript $name.R");
system("convert $name.RNA_PEP_Quadrant.pdf -density 300 $name.RNA_PEP_Quadrant.png");
__END__
欢迎讨论交流!