首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >转录组表达数据分析的一些可视化

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

作者头像
生信技能树
发布2018-03-29 16:08:17
1.6K0
发布2018-03-29 16:08:17
举报
文章被收录于专栏:生信技能树生信技能树

通过前面的讲解,我们顺利的了解了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语言处理表达矩阵呢?

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2018-03-21,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 表达矩阵的可视化
  • 表达矩阵的一些QC
  • 差异分析结果的火山图:
相关产品与服务
数据库
云数据库为企业提供了完善的关系型数据库、非关系型数据库、分析型数据库和数据库生态工具。您可以通过产品选择和组合搭建,轻松实现高可靠、高可用性、高性能等数据库需求。云数据库服务也可大幅减少您的运维工作量,更专注于业务发展,让企业一站式享受数据上云及分布式架构的技术红利!
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档