转录组表达数据分析的一些可视化

通过前面的讲解，我们顺利的了解了GEO数据库以及如何下载其数据，得到我们想要的表达矩阵，也学会了两个常用的套路分析得到的表达矩阵，就是GSEA分析和差异分析。也通过超几何分布检验的方法成功的理解了我们的统计学显著的差异表达基因的生物学功能。包括 GO/KEGG数据库以及 Reactome和Msigdb数据库的理解。

历史目录： 解读GEO数据存放规律及下载，一文就够解读SRA数据库规律一文就够从GEO数据库下载得到表达矩阵一文就够 GSEA分析一文就够（单机版+R语言版）根据分组信息做差异分析- 这个一文不够的差异分析得到的结果注释一文就够

但是我们的整个芯片数据分析流程居然缺少一个最重要的环节，就是可视化。

我们这个系列可以说是讲解的非常清楚了，但是无论再清楚的文字版教程，也总会有疏漏的地方，或者说传达不清楚的地方，理论上有视频教程是最好的了。

表达矩阵的可视化

这个可视化非常丰富：http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html

帮助我们理解数据的分布，分组，表达量高低等等。

首先加载一些R包

library(CLL)
library(ggplot2)
library(reshape2)
library(gpairs)
library(corrplot)

加载内置的测试数据：

data(sCLLex)
sCLLex=sCLLex[,1:8] ## 样本太多，我就取前面8个
group_list=sCLLex$Disease
exprSet=exprs(sCLLex)
head(exprSet)

主要用到ggplot2这个包，需要把我们的宽矩阵用reshape2包变成长矩阵

library(reshape2)
exprSet_L=melt(exprSet)
colnames(exprSet_L)=c('probe','sample','value')
exprSet_L$group=rep(group_list,each=nrow(exprSet))
head(exprSet_L)

这里就展现一个最简单的 boxplot

p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)

img

表达矩阵的一些QC

比较传统的就是主成分分析和相关性分析咯，可视化如下；

样本PAC

img

差异分析结果的火山图：

对芯片差异分析结果进行注释基因ID，并绘制火山图，判断具有统计学显著的差异基因列表。

img

火山图代码如下：

my_volcano <- function(DEG,t_p,t_FC=0,t_marker){
  ## example: print(my_volcano(limma_results[,c(5,1)],0.01,0.6,0))
  library(ggplot2)
  DEG=na.omit(DEG)
  colnames(DEG)=c('p','logFC')
  DEG$gene <- rownames(DEG)
  DEG[DEG$p < 1e-10,'p']=1e-10

  if (t_FC == 0) {
    t_FC <- with(DEG, mean(abs(logFC)) + 2 * sd(abs(logFC)))
  }

  #DEG$change <- as.factor(DEG$p<t_p & abs(DEG$logFC) > t_FC)
  DEG$change= as.factor(ifelse(DEG$p<t_p & abs(DEG$logFC) > t_FC, 
                                    ifelse(DEG$logFC > t_FC,"UP", "DOWN"), "NOT"))
  this_tile <- paste0("Cutoff for logFC is ", round(t_FC, 3), 
                      "nThe number of up gene is ", 
                      nrow(DEG[DEG$change == "UP", ]), 
                      "nThe number of down gene is ", 
                      nrow(DEG[DEG$change == "DOWN", ]))

  p <- ggplot(data=DEG, aes(x=logFC, y =-log10(p),color =change)) +
    geom_point() +
    scale_color_manual(values =c('blue',"black","red"))+
    geom_hline(yintercept = -log10(t_p),lty=4,lwd=0.6,alpha=0.8)+
    geom_vline(xintercept = c(t_FC,-t_FC),lty=4,lwd=0.6,alpha=0.8)+
    theme_bw()+
    theme(panel.border = element_blank(),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),   
          axis.line = element_line(colour = "black"))+ggtitle(this_tile) +
    labs( x="log2 (fold change)",y="-log10 (q-value)")+
    theme(plot.title = element_text(hjust = 0.5))
  if (t_marker !=0 ) p = p+geom_text(DEG=subset(DEG, abs(logFC) > t_marker), aes(label=gene),col="green",alpha = 0.5)
  return(p)
}

写到这里，我们的表达矩阵数据分析系列教程就告一段落啦，大家应该明白，我们的教程是针对表达矩阵，不仅仅是芯片得到的表达矩阵，也包括RNA-seq等技术得到的表达矩阵，所以我才会讲解GEO+SRA数据库。希望大家能活学活用。

看了一下，这个系列，大家的留言并不积极，不知道是不是follow的人比较少。

我们简单做一个调查吧，看看大家是否需要我们生信技能树团队录制视频教程来指导大家入门R语言处理表达矩阵呢？