其他
Cell 教我学画图之累积分布曲线
那就,端午节快乐吧
1前言
有些人 20 出头就已经相夫教子了，有些人 30 岁了却还在追求梦想，有的人为了车贷，房贷货一辈子，有些人买一辆摩托车就走遍大好河山，你到底要成为怎样的人，你要怎么活，为谁而活。---from online
2引言
模仿绘制一篇 cell 文章里的一个主图。文章讲的是 YTHDF1 识别 mRNA 上的 m6A 并促进 其 mRNA 的翻译。
文章标题:
N6-methyladenosine Modulates Messenger RNA Translation Efficiency
主图 2:
绘制 A,B,C 前三个图。
3绘图
读取并预处理数据
library(tidyverse)
library(ggplot2)
library(ggsci)
library(reshape2)
library(ggpubr)
library(patchwork)
# load data
df <- read.delim('data1.txt',header = T)
colnames(df)
# [1] "type" "gene_name" "ribo" "mRNA" "TE"
# wide to long
longdf <- melt(df,id.vars = c('type','gene_name'),
value.name = 'log2change',
variable.name = 'expr')
# check
head(longdf)
# type gene_name expr log2change
# 1 CLIP+IP target A4GALT ribo -1.306
# 2 Non-target AAAS ribo -0.354
# 3 Non-target AADAT ribo 0.318
# 4 Non-target AAED1 ribo -0.734
# 5 Non-target AAGAB ribo 0.932
# 6 CLIP only target AAK1 ribo 1.234
A 图
ribo <- longdf %>% filter(expr == 'ribo')
# summary
table(ribo$type)
# CLIP only target CLIP+IP target Non-target
# 3652 1261 2822
# stastic
compare_means(log2change~type,data = ribo,method = "wilcox.test",
ref.group = 'Non-target')
# A tibble: 2 x 8
# .y. group1 group2 p p.adj p.format p.signif method
# <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
# 1 log2change Non-target CLIP+IP target 1.13e-23 2.3e-23 <2e-16 **** Wilcoxon
# 2 log2change Non-target CLIP only target 5.46e-18 5.5e-18 <2e-16 **** Wilcoxon
# plot
pribo <- ggplot(ribo,
aes(x = log2change)) +
geom_hline(yintercept = 0.5,lty = 'dashed',size = 1,color = 'grey60') +
geom_vline(xintercept = 0,lty = 'dashed',size = 1,color = 'grey60') +
stat_ecdf(aes(color = type),geom = 'line',size = 1) +
theme_classic(base_size = 16) +
scale_color_futurama(name = '',label = c('CLIP only target' = 'CLIP only target (3652)',
'CLIP+IP target' = 'CLIP+IP target (1261)',
'Non-target' = 'Non-target (2822)')) +
scale_x_continuous(limits = c(-1.5,1.5),breaks = c(-1.5,0,1.5)) +
theme(aspect.ratio = 0.8,
legend.background = element_blank(),
legend.position = c(0.25,0.9),
strip.background = element_rect(colour = NA,fill = 'grey90')) +
facet_wrap(~expr,ncol = 3,scales = 'free') +
xlab('log2(siYTHDF1/siControl)') +
ylab('Cumulative Fraction') +
annotate(geom = 'text',x = 1,y = 0.15,
label = bquote('P = 5.46 x 10'^-18),size = 6,color = 'orange') +
annotate(geom = 'text',x = 1,y = 0.05,
label = bquote('P = 1.13 x 10'^-23),size = 6,color = 'red')
4B 图
#######################################################################
rna <- longdf %>% filter(expr == 'mRNA')
# stastic
compare_means(log2change~type,data = rna,method = "wilcox.test",
ref.group = 'Non-target')
# A tibble: 2 x 8
# .y. group1 group2 p p.adj p.format p.signif method
# <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
# 1 log2change Non-target CLIP+IP target 1.06e-22 2.1e-22 <2e-16 **** Wilcoxon
# 2 log2change Non-target CLIP only target 7.96e- 6 8 e- 6 8e-06 **** Wilcoxon
# plot
prna <- ggplot(rna,
aes(x = log2change)) +
geom_hline(yintercept = 0.5,lty = 'dashed',size = 1,color = 'grey60') +
geom_vline(xintercept = 0,lty = 'dashed',size = 1,color = 'grey60') +
stat_ecdf(aes(color = type),geom = 'line',size = 1,show.legend = F) +
theme_classic(base_size = 16) +
scale_color_futurama(name = '',label = c('CLIP only target' = 'CLIP only target (3652)',
'CLIP+IP target' = 'CLIP+IP target (1261)',
'Non-target' = 'Non-target (2822)')) +
scale_x_continuous(limits = c(-1.5,1.5),breaks = c(-1.5,0,1.5)) +
theme(aspect.ratio = 0.8,
legend.background = element_blank(),
legend.position = c(0.25,0.9),
strip.background = element_rect(colour = NA,fill = 'grey90')) +
facet_wrap(~expr,ncol = 3,scales = 'free') +
xlab('log2(siYTHDF1/siControl)') +
ylab('Cumulative Fraction') +
annotate(geom = 'text',x = 1,y = 0.15,
label = bquote('P = 8.0 x 10'^-6),size = 6,color = 'orange') +
annotate(geom = 'text',x = 1,y = 0.05,
label = bquote('P = 1.1 x 10'^-22),size = 6,color = 'red')
C 图
#######################################################################
te <- longdf %>% filter(expr == 'TE')
# stastic
compare_means(log2change~type,data = te,method = "wilcox.test",
ref.group = 'Non-target')
# A tibble: 2 x 8
# .y. group1 group2 p p.adj p.format p.signif method
# <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
# 1 log2change Non-target CLIP+IP target 1.11e-16 2.2e-16 < 2e-16 **** Wilcoxon
# 2 log2change Non-target CLIP only target 1.12e-15 1.1e-15 1.1e-15 **** Wilcoxon
# plot
pte <- ggplot(te,
aes(x = log2change)) +
geom_hline(yintercept = 0.5,lty = 'dashed',size = 1,color = 'grey60') +
geom_vline(xintercept = 0,lty = 'dashed',size = 1,color = 'grey60') +
stat_ecdf(aes(color = type),geom = 'line',size = 1,show.legend = F) +
theme_classic(base_size = 16) +
scale_color_futurama(name = '',label = c('CLIP only target' = 'CLIP only target (3652)',
'CLIP+IP target' = 'CLIP+IP target (1261)',
'Non-target' = 'Non-target (2822)')) +
scale_x_continuous(limits = c(-1.5,1.5),breaks = c(-1.5,0,1.5)) +
theme(aspect.ratio = 0.8,
legend.background = element_blank(),
legend.position = c(0.25,0.9),
strip.background = element_rect(colour = NA,fill = 'grey90')) +
facet_wrap(~expr,ncol = 3,scales = 'free') +
xlab('log2(siYTHDF1/siControl)') +
ylab('Cumulative Fraction') +
annotate(geom = 'text',x = 1,y = 0.15,
label = bquote('P = 1.1 x 10'^-15),size = 6,color = 'orange') +
annotate(geom = 'text',x = 1,y = 0.05,
label = bquote('P = 1.1 x 10'^-16),size = 6,color = 'red')
合并
加了分面和两条虚线美化一下:
# combine
pribo + prna + pte
5结尾
将此推文 分享到三个微信群里
(老俊俊生信群除外), 截图发我微信, 我把数据和代码分享给你哦!
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群满200人后需收取20元入群费用), 数据代码已上传至QQ群,欢迎加入下载。
群二维码：
老俊俊微信：
知识星球：
所以今天你学习了吗？
今天的分享就到这里了，敬请期待下一篇！
最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！
如果觉得对您帮助很大，赏杯快乐水喝喝吧！
往期回顾
◀Molecular Cell 文章 ribosome pausing 结果复现 (终)
◀Molecular Cell 文章 ribosome pausing 结果复现 (四)
◀Molecular Cell 文章 ribosome pausing 结果复现 (三)
◀Molecular Cell 文章 ribosome pausing 结果复现 (二) (PCR 去重)
◀Molecular Cell 文章 ribosome pausing 结果复现 (一)
◀...