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

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

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-03-21

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏新智元

【并非愚人节】科学家创建可自我复制的神经网络,AI像生命体一样繁殖

1873
来自专栏范传康的专栏

使用ARKit开发AR游戏-基础篇:ARKit入门篇

11月到12月,使用新出的ARkit开发了一个AR游戏,梳理下开发过程的经验,整理成文。 计划是一个系列,入门篇主要是收集的资料整合。

4574
来自专栏新智元

让神经网络替你编程:如何用深度学习实现程序自动合成

【新智元导读】本文介绍了训练神经网络学习用复杂的函数式语言(FlashFill DSL)进行编程取得的成功,标志着神经程序合成方面一个令人兴奋的突破。 ● 作...

4198
来自专栏华章科技

数据分析图的十大错误,你占了几个?

优秀的数据可视化依赖优异的设计,并非仅仅选择正确的图表模板那么简单。全在于以一种更加有助于理解和引导的方式去表达信息,尽可能减轻用户获取信息的成本。当然并非所有...

861
来自专栏IT大咖说

别急!看完文章再来说你懂TensorFlow

1641
来自专栏PPV课数据科学社区

【必看工具】可视化图表表达的10个错误。

数据可视化是一个沟通复杂信息的强大武器。通过可视化信息,我们的大脑能够更好地抓取和保存有效信息,增加信息的印象。但如果数据可视化做的较弱,反而会带来负面效果。错...

2956
来自专栏月色的自留地

从锅炉工到AI专家(11)(END)

2217
来自专栏机器学习实践二三事

【Google 年度顶级论文】机器学习系统,隐藏多少技术债?

原文在此:google原文 1. 介绍 随着机器学习(ML)社群持续积累了几年对于活跃系统(live systems)的经验,一种让人不舒服的趋势广泛地浮出...

26610
来自专栏AI科技大本营的专栏

观点 | 哈哈,TensorFlow被吐槽了吧

作者 | Nico 参与 | shawn 今天,一篇吐槽TensorFlow的文章在网上刷屏,到底是怎么回事呢?来看这位作者的抱怨有没有道理。 每隔几个月,我都...

29711
来自专栏机器人网

一个简单的多机器人编队算法实现--PID

用PID进行领航跟随法机器人编队控制 课题2:多机器人编队控制 研究对象:两轮差动的移动机器人或车式移动机器人 研究内容:平坦地形,编队的保持和避障,以及避障和...

4817

扫码关注云+社区

领取腾讯云代金券