朝气蓬勃 2019-03-18
R语言作图点滴积累,今天要记录ggplot2中散点图的做法。散点图算是数据展示中非常基础的一种方法了吧,一般用于展示两个变量之间的关系(比如线性相关)。按照惯例,我每次讲一个新的图都会用一个我实际使用的例子,先来说说今天做这个散点图的例子吧:
问题描述:
我有两个变量,一个变量是蛋白相互作用网络中节点的度(degree),另一个变量是这些蛋白参与形成的复合物的个数;现在我想看看,一个蛋白参与形成的复合物个数(Number of complexes)是否跟它的度(degree)相关。
上面描述的这个问题,就是一个很典型的可以通过散点图来展示的例子。先来个数据快照(显示前10个数据点):
漂亮的散点图要一步一步来
1. 首先来个简单的散点图
#Import data > dat <- read.table("ProteinDegree_complex.txt",header = TRUE) #plot a simple scatter plot > library(ggplot2) > p <- ggplot(dat,aes(x=degree,y=complex)) + geom_point(shape=19)) + xlab("Degree") + ylab("Number of complexes") > p
note:
2. 如果想要拟合一条直线呢
#Import data > dat <- read.table("ProteinDegree_complex.txt",header = TRUE) #plot a simple scatter plot > library(ggplot2) > p <- ggplot(dat,aes(x=degree,y=complex)) + geom_point(shape=19)) + xlab("Degree") + ylab("Number of complexes") + geom_smooth(method = lm) > p
ggplot2 提供一个函数自动添加拟合的曲线(包括直线),当然该函数底层肯定是做了拟合分析的,比如线性回归分析等。
下面两个图分别是使用了"lm"和"loess":
拟合直线
拟合曲线
给散点图加了直线,可是这个直线拟合得怎么样,以及拟合的直线的参数ggplot2并没有提供,为了图的信息更完整,我们应当考虑给这个拟合的直线加上公式,以及拟合的R2值。这样之后,我们从图上可以得到些什么信息呢?
首先,我们可以很容易知道degree和complex数目是呈正相关关系的,通过拟合直线和公式可以知道两个变量的线性关系强弱;然后,然后就是一堆不知道是啥的黑点...
总感觉还是缺少些什么,仔细看看这个散点图,你会不会想知道图中degree很高且参与很多复合物的这几个蛋白是什么?这样的蛋白一定是生物细胞中十分重要的蛋白。所以呢,我们可以把最靠近右上角的前10个点给高亮出来,甚至给这些点表示label(基因名)。说干就干,看代码:
#Import data > dat <- read.table("WD40_complex_degree.out",header = TRUE) #edit the formula for the fitted line > formula <- sprintf("italic(y) == %.2f %+.2f * italic(x)", round(coef(dat.lm)[1],2),round(coef(dat.lm)[2],2)) r2 <- sprintf("italic(R^2) == %.2f",summary(dat.lm)$r.squared) labels <- data.frame(formula=formula,r2=r2,stringsAsFactors = FALSE) #plot the simple scatterplot > p <- ggplot(dat,aes(x=degree,y=complex,colour=degree>=63)) + geom_point(shape=19) + xlab("Degree of WD40 proteins") + ylab("Number of complexes") #linear regression analysis > dat.lm <- lm(complex ~ degree, data = dat) #add a line and labels for the formula > p <- p + geom_abline(intercept = coef(dat.lm)[1],slope = coef(dat.lm)[2]) + geom_text(data=labels,mapping=aes(x = 15,y=175,label=formula),parse = TRUE,inherit.aes = FALSE, size = 6) + geom_text(data=labels,mapping=aes(x = 15,y=165,label=r2),parse = TRUE,inherit.aes = FALSE, size = 6) + #add labels(gene name) for top 10 degree-ranked proteins annotate(geom = "text",x=annoText$degree-1,y=annoText$complex-2,label=annoText$WD40id, size=4.0) > p + theme(legend.position = "none") + theme(axis.title = element_text(size = 16), axis.text = element_text(size = 12,colour="black"))
代码一下子长了好多(囧),我们可以与前面的比较下,
最后的图是这个样子的:
最后的样子
pytyhon学习资料
python学习资料