这是ggplot2可视化专题的第一个实例操作

ggplot2的基本思路见前文总论：基于ggplot2的RNA-seq转录组可视化：总述和分文目录

【ggplot2绘图具体策略第二篇】：单/多个基因在组间同时展示的多种选择：分组小提琴图、分组/分面柱状图、单基因蜂群点图拼图绘制

【ggplot2绘图具体策略第三篇】：配对样本基因表达趋势：优化前后的散点连线图+拼图绘制

在完成差异表达分析后，我们通常想使用图片对分组间的所有测到的基因表达情况进行总体的展示。此时，组间标准化counts、或FPKM或TPM的相关性散点图即可以完成这一目的，在展示所有基因在两组间的表达水平的同时，可以选择性突出我们重点关注的hits。

此外在展示差异基因的层面，除了热图能够通过色阶反映基因在样本间的表达离散程度外，瀑布图更能直观地反映差异基因的foldchange以及整体的变化趋势。

此前在文献阅读过程中发现几张figures还挺喜欢的，分别是测序后FoldChange和p值的总体展示瀑布图、测序结果间相关性高密度散点图，和带errorbar注释的测序结果间相关性散点图。

the figures that we're going to imitate

致谢以下研究人员，引用如有不妥请联系我删除。

A图来源：Ye, L.; Park, J.J.; Dong, M.B.; Yang, Q.; Chow, R.D.; Peng, L.; Du, Y.; Guo, J.; Dai, X.; Wang, G., et al. In vivo CRISPR screening in CD8 T cells with AAV-Sleeping Beauty hybrid vectors identifies membrane targets for improving immunotherapy for glioblastoma. Nat Biotechnol 2019, 37, 1302-1313, doi:10.1038/s41587-019-0246-4.

B图来源：Zhao, Y.; Tyrishkin, K.; Sjaarda, C.; Khanal, P.; Stafford, J.; Rauh, M.; Liu, X.; Babak, T.; Yang, X. A one-step tRNA-CRISPR system for genome-wide genetic interaction mapping in mammalian cells. Sci Rep 2019, 9, 14499, doi:10.1038/s41598-019-51090-3.

C图来源：Chow, R.D.; Wang, G.; Ye, L.; Codina, A.; Kim, H.R.; Shen, L.; Dong, M.B.; Errami, Y.; Chen, S. In vivo profiling of metastatic double knockouts through CRISPR-Cpf1 screens. Nat Methods 2019, 16, 405-408, doi:10.1038/s41592-019-0371-5.

下面就尝试来用前文介绍的TCGA白人LUSC肺鳞癌mRNA-seq公共数据复制一下图片大致的形式。

数据获取

基于第一篇文章从TCGA数据库下载并整合清洗高通量肿瘤表达谱-临床性状数据，我们下载并清洗TCGA数据库中white人种的LUSC肺鳞癌mRNA-seq转录组counts数据和FPKM数据。

随后根据第二篇文章TCGA数据整合后进行DESeq2差异表达分析和基于R的多种可视化对counts数据进行了基于DESeq2的差异分析。

现在假设我们已经获得：
（样本是1到344为cancer，345到386为normal）

（1）dds_DE: the object generated by DESeq2 containing all info about the differential expression analysis.
（2）resSigAll: the object containing all up-regulated and down-regulated significant genes extracted by results and subset function from dds_DE.
（3）resSigUp: the subset object of resSigAll containing all up-regulated genes.
（4）resSigDown: the subset object of resSigAll containing all down-regulated genes.
（5）expr_vst: vst-transformed normalized counts matrix generated by DESeq2.

需要R包：

ggplot2 (作图)、DESeq2 (从差异分析对象中提取一些数据)、ggrepel (不会发生重叠的基于ggplot2的文字注释包)、IDPmisc（绘制像素化高密度散点图），和dplyr (进行文字处理)。

install.packages('ggplot2')
install.packages('DESeq2')
install.packages('ggrepel')
install.packages('dplyr')
install.packages('IDPmisc')
library('ggplot2')
library('DESeq2')
library('ggrepel')
library('dplyr')
library('IDPmisc')

基因表达量rank瀑布图（A）

我们首先需要想好ggplot2需要的数据是什么。观察图A，我们可知绘图需要：
（1）所有在DESeq2中分析的qualified genes的log2FC（x轴）；
（2）基因名和log2FC ranks（y轴）；
（3）基因FC的校正后p值（点的size）。

因此使用DESeq2包的results函数获得（几乎）所有基因的基因名、log2FC、-log10(p.adj)、表达rank，并按rank倒序排列。由于需要用颜色突出top5和bottom5差异基因，因此dataframe还需加一列'hits'将基因分为top5/medium/bottom5。准备将该数据框传递给画图主体函数。

由于只对top/bottom5差异基因进行文字注释，因此还需要准备一个只包含它们的数据框传递给注释函数。

pre_ranked_all_genes <- as.data.frame(results(dds_DE, alpha=0.9999,contrast=c("sample_type","cancer","normal")))
pre_ranked_all_genes <- pre_ranked_all_genes[order(pre_ranked_all_genes$log2FoldChange, decreasing = F),]
trend_all <- sapply(pre_ranked_all_genes$log2FoldChange, function(x){if(x>0) 'up' else 'down'})
pre_ranked_all_genes <- data.frame('gene_names'=rownames(pre_ranked_all_genes), pre_ranked_all_genes, 
                                   'trend'=trend_all, 'rank'=1:nrow(pre_ranked_all_genes), 
                                   'hits'='medium', stringsAsFactors = F)
pre_ranked_all_genes[1:5, 'hits'] <- 'bottom5'
pre_ranked_all_genes[(nrow(pre_ranked_all_genes)-4):nrow(pre_ranked_all_genes), 'hits'] <- 'top5'
pre_ranked_all_genes$padj <- -log10(pre_ranked_all_genes$padj)
set.seed(12)
to_be_pointed_out_all <- pre_ranked_all_genes[c(1:5, (nrow(pre_ranked_all_genes)-4):nrow(pre_ranked_all_genes)),]

随后可进行绘图，思路和使用参数解释如下（可对比代码理解）：

◈ 含所有基因的数据框传递给ggplot函数 (data=)，x轴为rank，y轴为log2FoldChange，color以hits列分类 (aes())，作为绘图对象中的全局变量。

◈ 调用散点图函数 (geom_point()) 并在继承变量的基础上补充参数，即本函数控制的散点图中点的size由padj列控制 (aes())。

◈ 添加y轴截距为2，0和-2的水平虚线 (geom_hline(xintercept=, size=, linetype=))，添加x轴位于上/下调基因间的竖直虚线 (geom_vline(yintercept=, ...))。

◈ 选择不继承ggplot函数的输入变量 (inherit.aes=F)，将仅含top/bottom5基因信息的小dataframe传入ggrepel包的文字注释函数 (geom_text.repel(data=,mapping=))，指定与全局相同的x/y轴、注释文字和颜色信息 (aes(x=, y=, label=, color=))，并限定注释大小 (size)，排列方式 (direction='x'/'y')，和注释位置的x轴范围限制 (xlim=c(m,n))。

◈ 手动设置ggplot函数传入color分类变量的对应颜色 (scale_color_manual(values=c('type1'='color1', ...)))；设置size对应图例的名称 (scale_size_continuous(name=))。

◈ 最后设置x轴小标题(xlab())，在默认白色标题(theme_bw())的基础上修改分割线为白色，并居中图例标题(theme())。

fall_plot <- ggplot(pre_ranked_all_genes, aes(x=rank, y=log2FoldChange, color=hits))+
  geom_point(aes(size=padj))+
  geom_hline(yintercept = c(2,-2), linetype=2, size=0.25)+
  geom_hline(yintercept = c(0), linetype=1, size=0.5)+
  geom_vline(xintercept = sum(pre_ranked_all_genes$trend == 'down')+0.5, linetype=2, size=0.25)+
  ggrepel::geom_text_repel(inherit.aes = F, data = to_be_pointed_out_all, aes(x=rank, y=log2FoldChange, label=gene_names, color=hits), 
                           size=3, direction = 'y', xlim = c(1000, 15000))+
  scale_color_manual(values=c('bottom5'='blue', 'top5'='red', 'medium'='black'))+
  scale_size_continuous(name='-log10(FDR)')+
  scale_y_continuous(breaks=c(-5, -2, 0, 2, 5, 10))+
  xlab('rank of differentially expressed genes')+
  theme_bw()+
  theme(panel.grid = element_line(color='white'), legend.title.align = 0.5)
fall_plot

感觉还原度还挺不错的！

the plot we get

另一个版本的只含差异基因的瀑布图：

这个版本我们只选取差异基因作图，随机选取5个基因进行注释，颜色代表foldchange，而不再将size与校正p值关联。

整体流程区别不大。作图的代码可以尝试自行理解。差别主要体现在：
ggplot函数中aes()映射的color变量由上面的三分类变量（top5/medium/bottom5）变成了连续变量（log2FoldChange），因此颜色的设置发生变化。

这里使用了scale_color_gradient2()设置颜色，其特点为可以设置3种颜色 (low=, high=, mid=, midpoint=) 根据指定值由低到高的渐变。

#draw the foldchange plot with specific annotations.
pre_ranked_sig_genes <- as.data.frame(resSigAll)[,c(2,6)] %>% data.frame('gene_names'=resSigAll@rownames, stringsAsFactors = F)
pre_ranked_sig_genes <- pre_ranked_sig_genes[order(pre_ranked_sig_genes$log2FoldChange, decreasing = T),]
trend <- sapply(pre_ranked_sig_genes$log2FoldChange, function(x){if(x>0) 'up' else 'down'})
pre_ranked_sig_genes <- data.frame(pre_ranked_sig_genes, 'trend'=trend,
                                   'rank'=1:nrow(pre_ranked_sig_genes), stringsAsFactors = F)
#randomly select 5 genes to be annotated in the plot.
set.seed(12)
to_be_pointed_out <- pre_ranked_sig_genes[sample(1:nrow(pre_ranked_sig_genes), 5),]

fall_plot <- ggplot(pre_ranked_sig_genes, aes(x=rank, y=log2FoldChange, color=log2FoldChange))+
  geom_point(size=1)+
  geom_hline(yintercept = c(2,-2), linetype=2, size=0.25)+
  geom_hline(yintercept = c(0), linetype=1, size=0.5)+
  geom_vline(xintercept = 1636.5, linetype=2, size=0.25)+
  scale_color_gradient2(low="navy", high="firebrick3", mid="white", midpoint=0)+
  geom_point(inherit.aes = F, data=to_be_pointed_out, aes(x=rank, y=log2FoldChange), size = 3, color='black')+
  geom_point(inherit.aes = F, data=to_be_pointed_out, aes(x=rank, y=log2FoldChange), size = 2, color='yellow')+
  ggrepel::geom_text_repel(inherit.aes = F, data = to_be_pointed_out, aes(x=rank, y=log2FoldChange, label=gene_names), size=3)+
  xlab('rank of differentially expressed genes')+
  theme_bw()+
  theme(panel.grid = element_line(color='white'), legend.title.align = 0.5)
fall_plot

得到的图形为:

another version of FoldChange plot

高密度表达相关性散点图（B）

这个figure的意义是为了说明测序结果的可重复性和可比较性，即主要的散点集中在散点图的对角线上，而色域能很好地反映这一聚集特点。

绘图需要：
（1）在cancer样本中的所有基因表达值（这里用vst counts）均数（y轴）；
（2）在normal样本中的所有基因表达值（这里用vst counts）均数（x轴）；
（3）散点的聚集密度信息（色域）。

其实已经有其他的R包或ggplot2中的函数可以完成类似的功能。但感觉没有达到article中的效果。article中的figure是直接对散点进行上色，而不是像素化的效果：

（1）IDPmisc包：iplot函数将画图区域按pixs命令的分辨率划分为若干方格，按落在每个方格中的点密度进行绘图。

pixs参数设置像素边长，colramp参数设置渐变色彩。

expr_for_visualization <- data.frame('cancer'=rowMeans(expr_vst[,1:344]), 'normal'=rowMeans(expr_vst[,345:386]), stringsAsFactors = F)
library(IDPmisc)
IDPmisc::iplot(x = expr_for_visualization$normal, y=expr_for_visualization$cancer, pixs = 0.4, 
               colramp=colorRampPalette(colors = c('black','blue','yellow','red')),
               xlab='normal', ylab='cancer', main='correlation scattorplot with density')

（2）ggplot2包stat_bin2d()函数：将画图区域分为30*30的方格进行密度计算和可视化。

expr_for_visualization <- data.frame('cancer'=rowMeans(expr_vst[,1:344]), 'normal'=rowMeans(expr_vst[,345:386]), stringsAsFactors = F)
ggplot(expr_for_visualization, aes(x=normal, y=cancer))+stat_bin2d()+theme_bw()+
  theme(panel.grid = element_line(colour = 'white'))

得到的结果如下，都是像素化效果。左边还好，右边丑哭：

density plot drawn by IDPmisc (left)-density plot drawn by ggplot2-stat_bin2d() function (right)

所以我自己用ggplot2进行了自定义作图：

图片的绘制思想：绘制表达均值的散点图，以散点密度为连续变量定义散点的color。

具体做法1：直接将绘图区域划分网格并计算网格内密度

设定将绘图区域分为70*70的网格，使用cut函数进行x/y轴的分组，使用table函数进行网格密度频数统计得到freq，最终merge函数左合并至坐标dataframe，让每个散点都获得自己所在网格的密度数。

#Draw the high-density scatterplot.
intended_pixel_num = 70
expr_for_visualization <- data.frame('cancer'=rowMeans(expr_vst[,1:344]), 'normal'=rowMeans(expr_vst[,345:386]), stringsAsFactors = F)
cancer_notches <- seq(from = min(expr_for_visualization$cancer-0.1),to = max(expr_for_visualization$cancer+0.1), length.out = (intended_pixel_num + 1))
normal_notches <- seq(from = min(expr_for_visualization$normal-0.1),to = max(expr_for_visualization$normal+0.1), length.out = (intended_pixel_num + 1))
cancer_groups <- as.character(cut(expr_for_visualization$cancer, breaks = cancer_notches))
normal_groups <- as.character(cut(expr_for_visualization$normal, breaks = normal_notches))
expr_for_visualization <- data.frame('genes'=rownames(expr_for_visualization), expr_for_visualization,
                                     'groups'=paste(cancer_groups, normal_groups, sep = '_'), stringsAsFactors = F)
frequency <- as.data.frame(table(expr_for_visualization$groups))
frequency$Var1 <- as.character(frequency$Var1)
colnames(frequency) <- c('groups','freq')
expr_for_visualization <- merge(expr_for_visualization, frequency, by='groups', all.x = T)

然后进行ggplot2的绘图：

◈ 含所有基因在cancer和normal样本间vst counts均值的数据框传递给ggplot函数，x轴为normal，y轴为cancer，color由freq列连续变量控制，作为绘图对象中的全局变量。

◈ 用geom_point继承输入变量作散点图，散点size固定。

◈ scale_color_gradientn()函数设置连续变量的多种颜色渐变。

◈ ggtitle()设置figure标题，统一使用自己习惯的作图主题，并将标题居中。

high_density_scatterplot <- ggplot(expr_for_visualization, aes(x=normal, y=cancer, colour=freq))+
  geom_point(size=0.1)+
  scale_color_gradientn(colours=c('black','blue','yellow','red'))+
  ggtitle('correlation scattorplot with density')+
  theme_bw()+theme(panel.grid = element_line(colour = 'white'), 
                   plot.title = element_text(hjust = 0.5))
high_density_scatterplot

获得的图片和第二种方法同时展示在下面。方便比较。

具体做法2：将绘图区域划分网格后，计算网格内和网格周围部分区域密度

可以看到，虽然是对散点着色，但在高密度区由于散点密集，分割像素化的痕迹还是很明显。所以这次扩大每个网格计算密度的范围，使得每个网格取密度时相互重叠，或许可以获得更好的过渡效果。

设定将绘图区域分为100*100的网格，使用两轮for循环（用for循环真是要了R的老命...），逐个网格进行分组标记和网格内散点密度数标记。

随后的ggplot2绘图命令几乎是完全相同的。不再过细解释。

intended_pixel_num <- 100
#data integration and formating.
expr_for_visualization <- data.frame('cancer'=rowMeans(expr_vst[,1:344]), 'normal'=rowMeans(expr_vst[,345:386]),
                                     'series'=rep('0', nrow(expr_vst)), 'counts'= rep(0, nrow(expr_vst)), stringsAsFactors = F)
cancer_notches <- seq(from = min(expr_for_visualization$cancer-0.1),to = max(expr_for_visualization$cancer+0.1), length.out = (intended_pixel_num + 1))
normal_notches <- seq(from = min(expr_for_visualization$normal-0.1),to = max(expr_for_visualization$normal+0.1), length.out = (intended_pixel_num + 1))
n <- 0
for(rows in 1:(length(cancer_notches)-1)) {
  for (columns in 1:(length(normal_notches)-1)) {
    n = n+1
    series_num <- as.character(n)
    cancer_interval <- cancer_notches[3]-cancer_notches[2]
    normal_interval <- normal_notches[3]-normal_notches[2]
    indices <- dplyr::intersect(which(expr_for_visualization$cancer>=(cancer_notches[rows]) & expr_for_visualization$cancer<(cancer_notches[rows+1])),
                     which(expr_for_visualization$normal>=(normal_notches[columns]) & expr_for_visualization$normal<(normal_notches[columns+1])))
    counts <- dplyr::intersect(which(expr_for_visualization$cancer>=(cancer_notches[rows]-cancer_interval*0.5) & expr_for_visualization$cancer<(cancer_notches[rows+1]+cancer_interval*0.5)),
                                which(expr_for_visualization$normal>=(normal_notches[columns]-normal_interval*0.5) & expr_for_visualization$normal<(normal_notches[columns+1]+normal_interval*0.5)))
    if(length(indices) >0){
      expr_for_visualization[indices,'series'] <- series_num
      expr_for_visualization[indices,'counts'] <- length(counts)
    }
  }
}

#draw the plot.
high_density_scatterplot <- ggplot(expr_for_visualization, aes(x=normal, y=cancer, colour=counts))+
  theme_bw()+theme(panel.grid = element_line(colour = 'white'), plot.title = element_text(hjust = 0.5))+
  geom_point(size=0.02)+ggtitle('correlation scattorplot with density')+
  scale_color_gradientn(colours=c('black','blue','yellow','red'))
high_density_scatterplot

得到的结果说明，第二种采样取密度的方法得到的结果很好地减轻了像素化的程度。
【不过我还是不知道article的原做法...它的色域渐变非常流畅！如果有朋友有相关经验，欢迎交流！】

the high-density scatterplot drawn by method1 (left)-the high-density scatterplot drawn by method2 (right)

注释特定基因及errorbar的表达相关性散点图（C）

图片的特点就在于对特定的outliers做了表达量的errorbar注释，增加了拟合直线和置信区间，以及左上角进行了少许图文注释。

绘图需要：
（1）在cancer样本中的所有基因表达值（这里用vst counts）均数（y轴）；
（2）在normal样本中的所有基因表达值（这里用vst counts）均数（x轴）；
（3）上下调差异基因列表（进行color标识）；
（4）top/bottom3基因在cancer和normal组中vst counts的均数、标准差sd（errorbar注释）。

数据整理：获得表达均值即在图中的坐标后，根据差异分析的结果将散点分为up/no/down三类。然后提取top/bottom3基因名，计算sd值，并与mean值整合成新的小dataframe用于特定注释。

#draw the scatterplot with specific annotation.
vst_coordinate <- data.frame('cancer'=rowMeans(expr_vst[,1:344]), 'normal'=rowMeans(expr_vst[,345:386]), 'trend'='no', stringsAsFactors = F)
vst_coordinate[resSigUp@rownames,'trend'] <- 'up'
vst_coordinate[resSigDown@rownames,'trend'] <- 'down'
#get the top 3 up/down-regulated genes and calculate the sd values.
ressigall_df <- as.data.frame(resSigAll)
ressigall_df <- ressigall_df[order(ressigall_df$log2FoldChange, decreasing = T),]
top_genes <- rownames(ressigall_df)[c(1:3, (nrow(ressigall_df)-2):nrow(ressigall_df))]
#get the small dataframe containing sd and mean vst values of top genes.
sd_2d_calculation <- function(gene_name){
  sd_cancer <- sd(expr_vst[gene_name,1:344])
  sd_normal <- sd(expr_vst[gene_name,345:386])
  data.frame('sd_cancer'=sd_cancer,'sd_normal'=sd_normal)
}
sd_df <- do.call(rbind,lapply(top_genes, sd_2d_calculation))
sd_df <- cbind(vst_coordinate[top_genes,], sd_df)

进行ggplot2作图：
◈ 含所有基因在cancer和normal样本间vst counts均值的数据框传递给ggplot函数，x轴为normal，y轴为cancer，作为绘图对象中的全局变量。

◈ 用stat_smooth()继承输入变量进行线性回归拟合 (method='lm')，添加拟合虚线 (linetype=2) 并可视化95%置信区间 (level=0.95)。

◈ geom_point()函数绘制主体散点，并补充参数，以trend列作为散点颜色分类依据。

◈ scale_color_manual()函数手动设置散点类别对应颜色，上调基因用firebrick3，下调基因用navy，不显著基因用black。

◈ geom_errorbar()函数和geom_errorbarh()函数均不继承全局参数，而用小dataframe作为输入数据，分别添加top genes竖直和水平方向的errorbar (误差棒的方向需设置最大值和最小值，如errorbar函数为ymax, ymin；errorbarh函数为xmax，xmin)，并设置bar顶端横线的长度 (width/height) 和bar的粗细 (size)。

◈ 另外调用一个geom_point()函数，输入top genes dataframe进行top genes的黄色点标注。

◈ 使用ggrepel包的geom_label_repel()函数，输入top genes dataframe进行top hits的文字注释。（由于散点较多，所以使用了带文字框的label函数，而不是仅进行文字注释的text函数）

◈ 使用annotate()函数进行文字的小范围自由注释，需输入好注释的类型（如'text'）、color、size、label (注释内容)，以及在x/y轴分别的位置。

◈ 设置figure标题，统一使用自己习惯的作图主题，并将标题居中。

p <- ggplot(vst_coordinate, aes(x=normal, y=cancer))
#add the fitting line. level is the level of confidence interval.
p <- p+stat_smooth(method='lm', level = 0.95, linetype=2, color='black')
#add the points.
p <- p+geom_point(aes(color=trend), size=0.5)+
  scale_color_manual(values = c('no'='black','up'='firebrick3','down'='navy'))
#add error bars to top hits. widths and heights are parameters of the tips at the top of bars.
p <- p+geom_errorbar(inherit.aes = F, data = sd_df, aes(x=normal, ymin=cancer-sd_cancer, ymax=cancer+sd_cancer), width=0.3)+
  geom_errorbarh(inherit.aes = F, data = sd_df, aes(y=cancer, xmin=normal-sd_normal, xmax=normal+sd_normal), size=0.5, height=0.4)
#emphasize the top hits.
p <- p+geom_point(inherit.aes = F, data = sd_df, aes(x=normal, y=cancer), size=2.5, color='black')+
  geom_point(inherit.aes = F, data = sd_df, aes(x=normal, y=cancer), size=1.5, color='yellow')
#annotate the hits with gene names.
p <- p+geom_label_repel(data = sd_df, aes(label = rownames(sd_df)), 
                   color='black', fontface='bold', size = 3, 
                   segment.color = 'black', xlim = c(7, 20))
#add some text annotations on the upper left of the panel.
p <- p+annotate('text', label=paste0('Upregulated genes (n=', length(resSigUp@rownames),')'), color='firebrick3', x=9.6, y=19.5)+
  annotate('text', label=paste0('Downregulated genes (n=', length(resSigDown@rownames),')'), color='navy', x=10, y=18.7)+
  annotate('point', color='black', size=3, x=5.8, y=17.9)+
  annotate('point', color='yellow', size=2, x=5.8, y=17.9)+
  annotate('text', label='Top3 up/downregulated DE hits', color='black', x=10.5, y=17.9)
#Set the theme of the plot.
p <- p+theme_bw()+
  theme(panel.grid = element_line(color = 'white'), legend.position = 'none')+
  xlab('vst counts of normal samples')+ylab('vst counts of cancer samples')
p

获得的图片结果为：

scattorplot with specific annotations.jpg

因为是转录组mRNA-seq数据，不是原文中的较小范围的CRISPR-screening，所以图片没有那么美观。但是样式模仿到了就好！

不知道大家是否对ggplot2的绘图范式有了一定的认识，希望本文能有所帮助~喜欢就点个赞吧！