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 datadf <- read.delim('data1.txt',header = T)colnames(df)# [1] "type"      "gene_name" "ribo"      "mRNA"      "TE"# wide to longlongdf <- melt(df,id.vars = c('type','gene_name'),               value.name = 'log2change',               variable.name = 'expr')# checkhead(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')# summarytable(ribo$type)# CLIP only target   CLIP+IP target       Non-target# 3652             1261             2822# stasticcompare_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# plotpribo <- 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')# stasticcompare_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# plotprna <- 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')# stasticcompare_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# plotpte <- 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')

合并

加了分面和两条虚线美化一下:

# combinepribo + prna + pte

5结尾