频率分布直方图，或者频数分布直方图，是一种用来可视化数据的分布情况的绘图，在生物医学领域应用广泛，比如展示高通量测序结果的测序读数分布等。
Spatiotemporally-resolved mapping of RNA binding proteins via functional proximity labeling reveals a mitochondrial mRNA anchor promoting stress recovery
例如上面这个图展示的是两个定量蛋白组样品的频率分布直方图，有意思的是该图在x轴上下两个方向展示两组样品，在区分两组的同时又能很好的比较二者的差异。
今天我们就来学习一下频率分布直方图的画法。
用的测试数据是我从附表里面随便选的两列，不同基因的表达值。
library(ggplot2)
library(reshape2)
library(tidyverse)
data <- read.table("data.txt",header=T,sep="\t")
data_new <- melt(data,id="ID") #还是前面学过的语法，长矩阵转化成短矩阵。
colnames(data_new) <- c("ID","Sample","Value")
先来一个简单版本的，一个变量的情况。我们先用hist函数测试。
X128 <- data_new %>% filter(Sample=="X128N")
X130 <- data_new %>% filter(Sample=="X130C")
hist(X128$Value,
breaks = 14, #指定直方图的X轴区间，可以是向量分割自己指定
col = "red",
xlab = "Fold Change(Log2)",
ylab = "Frequency",
main = "test",
border = "black",
freq = FALSE,
density = 12,
angle = 45,
labels = T, #添加直方图bar上的label
ylim=c(0,0.8)
)
#添加密度线
lines(density(X128$Value),
col = "black",
lwd = 3)
#添加外框线
box()
下面我们还是测试最常用的ggplot。
ggplot(X128,aes(Value))+geom_histogram()
ggplot(X128,aes(Value))+
geom_histogram(stat = 'bin',bins = 20, #设定间距的个数
fill='darkgreen',
color='gray')+
theme_bw()
也可以通过设置 binwidth 参数的值，该参数值会覆盖 bins 参数的值，所以只要设置其中一个参数就可以了
ggplot(X128)+
geom_histogram(aes(Value),
stat = 'bin',bins = 20,
fill='darkgreen',
color='gray')+
geom_histogram(aes(Value,y = -..count..), #可以画出反方向的
stat = 'bin',bins = 20,
fill='blue',
color='gray')+
theme_bw()+
theme(axis.title = element_blank())
我们再来添加密度曲线上去。
ggplot(X128,aes(Value))+
geom_histogram(stat = 'bin',bins = 20,
fill='darkgreen',
color='gray')+
theme_bw()+
geom_freqpoly(bins = 20,binwidth = 0.5,size=1.5,color="red")+
theme(axis.title = element_blank())
ggplot(X128)+
geom_histogram(aes(Value),
stat = 'bin',bins = 20,
fill='darkgreen',
color='gray')+
geom_freqpoly(aes(Value),bins = 20,binwidth = 0.5,size=1.5,color="red")+
geom_histogram(aes(Value,y = -..count..),
stat = 'bin',bins = 20,
fill='blue',
color='gray')+
geom_freqpoly(aes(Value,y = -..count..),bins = 20,binwidth = 0.5,size=1.5,color="red")+
theme_bw()+
theme(axis.title = element_blank())
默认的geom_histogram是画的count，我们也可以通过density来画密度。
ggplot(X128)+
geom_histogram(aes(x=Value,y=..density..),
stat = 'bin',bins = 20,
fill='darkgreen',
color='gray')+
geom_histogram(aes(x=Value,y = -..density..),
stat = 'bin',bins = 20,
fill='blue',
color='gray')+
theme_bw()+
theme(axis.title = element_blank())
geom_freqpoly()的另一个方式是geom_density()，但底层密度计算是复杂的，从而导致有时结果很难解释，它们总是假设数据是连续的、无界的、平滑的。这两个函数是针对单个连续数值变量进行统计，但仍然可以比较不同的subgroup，举例：ggplot(diamonds, aes(price, fill = cut)) + geom_histogram(binwidth = 500)和ggplot(diamonds, aes(price, colour = cut)) + geom_freqpoly(binwidth = 500)，即histogram设置aes的fill参数，freqpoly设置aes的color参数。另一种可选方案当然是分面啦。
ggplot(X128)+
geom_histogram(aes(x=Value,y=..density..),
stat = 'bin',bins = 20,
fill='darkgreen',
color='gray')+
geom_density(aes(Value,y=..density..),bins = 20,binwidth = 0.5,size=1.5,color="red")+
geom_histogram(aes(x=Value,y = -..density..),
stat = 'bin',bins = 20,
fill='blue',
color='gray')+
geom_density(aes(Value,y = -..density..),bins = 20,binwidth = 0.5,size=1.5,color="red")+
theme_bw()+
theme(axis.title = element_blank())
我们还可以添加背景填充色，以及设置背景填充色。
dense=data.frame(density(X128$Value)[c("x","y")]) #获得密度分布数据
ggplot(X128)+
geom_histogram(aes(x=Value,y=..density..),
stat = 'bin',bins = 20,
fill='gray',
color='gray')+
geom_density(aes(Value,y=..density..),size=1.5,color="red")+
geom_area(data=subset(dense,x<2),aes(x,y,fill="Label 1"),alpha=0.4)+
geom_area(data=subset(dense,x>=2 & x<3),aes(x,y,fill="Label 2"),alpha=0.4)+
geom_area(data=subset(dense,x>=3 & x<5),aes(x,y,fill="Label 3"),alpha=0.4)+
geom_area(data=subset(dense,x>=5),aes(x,y,fill="Label 4"),alpha=0.4)+
scale_fill_manual("Test Tile",breaks=c("Label 1","Label 2","Label 3","Label 4"),
values=c("Label 1"="red","Label 2"="blue","Label 3"="purple","Label 4"="cyan")) #自定义颜色
接下来，我们来绘制多个变量的情况。
ggplot()+
geom_histogram(data=X128,aes(x=Value,y=..density..),
stat = 'bin',bins = 20,
fill='lightgreen',
color='gray')+
geom_density(data=X128,aes(Value,y=..density..),size=1.5,color="red")+
geom_histogram(data=X130,aes(x=Value,y = -..density..),
stat = 'bin',bins = 20,
fill='lightblue',
color='gray')+
geom_density(data=X130,aes(Value,y = -..density..),size=1.5,color="red")+
theme_bw()+
theme(axis.title = element_blank())
这样子，我们就绘制了一个镜像的直方图。
ggplot(data_new, aes(Value, after_stat(density), colour = Sample)) +
geom_freqpoly(bins = 40)
两组数据简单的密度曲线。
#多变量直方图
默认是堆积直方图的效果，和柱状图的调整是类似的，通过position来调整
ggplot(data_new, aes(Value, after_stat(density), fill= Sample)) +
geom_histogram(color="#e9ecef", alpha=0.6, position = 'identity') +
#geom_histogram(color="#e9ecef", alpha=0.6, position = 'stack') +
#geom_histogram(color="#e9ecef", alpha=0.6, position = 'dodge') +
#geom_histogram(color="#e9ecef", alpha=0.6, position = 'fill') +
scale_fill_manual(values=c("#377eb8", "#4daf4a"))
这是两个变量直方图放在一起的样子，还是不如镜像直方图直观。
也可以利用我们前面用过的分面技巧，分开绘制。
#分面直方图
ggplot(data_new, aes(Value, after_stat(density), fill= Sample)) +
geom_histogram(alpha = 0.6, bins = 40) +
geom_freqpoly(bins = 40)+
facet_wrap(~ Sample) +
theme(legend.position = "none")
下面，我们试着画一下，我们开始在paper中看到的图。
ggplot()+
geom_histogram(data=X128,aes(x=Value,y=..density..),
stat = 'bin',bins = 20,
fill='#2AC643',
color='white')+
geom_histogram(data=X130,aes(x=Value,y = -..density..),
stat = 'bin',bins = 20,
fill='gray60',
color='white')+
scale_y_continuous(label=abs)+
theme_classic(base_size = 15)+ #换个背景主题
scale_x_continuous(limits = c(1,7),
breaks = c(1,2,3,4,5,6,7),
expand = c(0,0))+
theme(panel.border = element_rect(size = 1,fill='transparent'),
legend.position = 'none', #去掉图例
axis.text = element_text(colour = 'black'))+
geom_vline(xintercept =median(X128$Value),linetype=2,cex=1)+ #添加辅助线
labs(x='X128/X130',y='Frequency')+ #自定义轴标题
annotate('text',x=median(X128$Value)+0.1,y=-0.6,
label = round(median(X128$Value),digits = 2),
size=4,color='black')+
annotate('text',x=2,y=0.7,label='Known nuclear RBPs',size=6,color='#2AC643')+
annotate('text',x=2,y=-0.7,label='Non-nuclear Non-RBPs',size=6,color='grey50')+ #添加文本标签
geom_segment(aes(y = 0.78, yend = 0.78,x=median(X128$Value), xend =median(X128$Value)+0.15),arrow = arrow(length = unit(0.2, "cm"),type="closed"),
size=0.5)+ #添加箭头
annotate('text',x=median(X128$Value)+0.3,y=0.78,
label ="retain",digits = 2,
size=4,color='black') #添加文本
网友评论