表达矩阵可视化大全

貌代码被折叠了,大家需要阅读原文才能复制粘贴我代码在Rstudio里面直接运行,几分钟就可以学会15个图的制作!

basic visualization for expression matrix

jmzeng1314@163.com

March 14, 2017

安装并加载必须的packages

如果你还没有安装,就运行下面的代码安装:

BiocInstaller::biocLite('CLL')install.packages('corrplot')install.packages('gpairs')install.packages('vioplot')

如果你安装好了,就直接加载它们即可

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

加载内置的测试数据:

data(sCLLex)sCLLex=sCLLex[,1:8] ## 样本太多,我就取前面8个

group_list=sCLLex$DiseaseexprSet=exprs(sCLLex)head(exprSet)
##           CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL CLL16.CEL
## 1000_at    5.743132  6.219412  5.523328  5.340477  5.229904  4.920686
## 1001_at    2.285143  2.291229  2.287986  2.295313  2.662170  2.278040
## 1002_f_at  3.309294  3.318466  3.354423  3.327130  3.365113  3.568353
## 1003_s_at  1.085264  1.117288  1.084010  1.103217  1.074243  1.073097
## 1004_at    7.544884  7.671801  7.474025  7.152482  6.902932  7.368660
## 1005_at    5.083793  7.610593  7.631311  6.518594  5.059087  4.855161
##           CLL17.CEL CLL18.CEL
## 1000_at    5.325348  4.826131
## 1001_at    2.350796  2.325163
## 1002_f_at  3.502440  3.394410
## 1003_s_at  1.091264  1.076470
## 1004_at    6.456285  6.824862
## 1005_at    5.176975  4.874563
group_list
## [1] progres. stable   progres. progres. progres. progres. stable   stable  
## Levels: progres. stable

接下来进行一系列绘图操作

主要用到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)
##       probe    sample    value    group
## 1   1000_at CLL11.CEL 5.743132 progres.
## 2   1001_at CLL11.CEL 2.285143 progres.
## 3 1002_f_at CLL11.CEL 3.309294 progres.
## 4 1003_s_at CLL11.CEL 1.085264 progres.
## 5   1004_at CLL11.CEL 7.544884 progres.
## 6   1005_at CLL11.CEL 5.083793 progres.

boxplot

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

vioplot

#library(vioplot)#?vioplot#vioplot(exprSet)#do.call(vioplot,c(unname(exprSet),col='red',drawRect=FALSE,names=list(names(exprSet))))p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()print(p)

histogram

p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)print(p)

density

p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density() print(p)

gpairs

library(gpairs)gpairs(exprSet
       #,upper.pars = list(scatter = 'stats') 
       #,lower.pars = list(scatter = 'corrgram')
       )

cluster

out.dist=dist(t(exprSet),method='euclidean')out.hclust=hclust(out.dist,method='complete')plot(out.hclust)

PCA

pc <- prcomp(t(exprSet),scale=TRUE)pcx=data.frame(pc$x)pcr=cbind(samples=rownames(pcx),group_list, pcx) p=ggplot(pcr, aes(PC1, PC2))+geom_point(size=5, aes(color=group_list)) +
  geom_text(aes(label=samples),hjust=-0.1, vjust=-0.3)print(p)

heatmap

choose_gene=names(sort(apply(exprSet, 1, mad),decreasing = T)[1:50])choose_matrix=exprSet[choose_gene,]choose_matrix=scale(choose_matrix)heatmap(choose_matrix)
library(gplots)
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
heatmap.2(choose_matrix)
library(pheatmap)pheatmap(choose_matrix)

DEG && volcano plot

library(limma)
## 
## Attaching package: 'limma'
## The following object is masked from 'package:BiocGenerics':
## 
##     plotMA
design=model.matrix(~factor(group_list))fit=lmFit(exprSet,design)fit=eBayes(fit)DEG=topTable(fit,coef=2,n=Inf)with(DEG, plot(logFC, -log10(P.Value), pch=20, main="Volcano plot"))      
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,                              ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
                       )this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',]))g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle( this_tile  ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red'))  ## corresponding to the levels(res$change)print(g)

ggplot画图是可以切换主题的

其实绘图有非常多的细节可以调整,还是略微有点繁琐的!

p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()print(p)
p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")p=p+theme_set(theme_set(theme_bw(base_size=20)))p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())print(p)

可以很明显看到,换了主题之后的图美观一些。

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

原文发表时间:2017-03-14

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏iOSDevLog

iTerm2 科学上网

683
来自专栏小小挖掘机

windows下使用word2vec训练维基百科中文语料全攻略!(一)

训练一个聊天机器人的很重要的一步是词向量训练,无论是生成式聊天机器人还是检索式聊天机器人,都需要将文字转化为词向量,时下最火的词向量训练模型是word2vec,...

2566
来自专栏杂七杂八

tensorflow-gpu版安装

需要环境 Anaconda CUDA cuDNN 注:tensorflow1.4用的是cuda8,cudnn6;tensorflow用的是cuda9,cudn...

2864
来自专栏漫漫深度学习路

pytorch学习笔记(十五):pytorch 源码编译碰到的坑总结

2017.11.17 最近打算学习一下 pytorch 源码,所以按照官网的教程从头编译了一下 pytorch 。在编译的过程中,碰到了两个坑,在这里记录一下。...

3187
来自专栏编程

厉害了,用Python一行代码实现人脸识别

摘要: 1行代码实现人脸识别,1. 首先你需要提供一个文件夹,里面是所有你希望系统认识的人的图片。其中每个人一张图片,图片以人的名字命名。2. 接下来,你需要准...

1728
来自专栏菩提树下的杨过

jQuery1.3以上版本"@"的问题

jQuery1.3.2已经发布好一段时间了,近日把原来的项目全部从jQuery1.2.6升级为1.3.2了.但是有一点要注意,1.3版以上的jQuery在根据选...

18110
来自专栏浅探ARKit

SceneKit动态加载.dae模型步骤详解

1.打开你的Xcode,在 /Contents/Developer/usr/bin/ 路径里找到 copySceneKitAssets 、 scntool 这2...

47413
来自专栏Golang语言社区

【Go 语言社区】把Go程序变小的办法

把Go程序变小的办法是: go build -ldflags “-s -w” (go install类似) -s去掉符号表(然后panic时候的stack tr...

3496
来自专栏祥子的故事

Tensorflow | win10中安装tensorflow-0.12.1 (0.12.1以后的版本安装均适用)

4097
来自专栏计算机视觉life

OpenCV学习入门(二):Image Watch神器

Image Watch是在visual studio 2012及以上版本上使用的一款OpenCV工具,能够在调试过程中实时显示内存中矩阵Mat(存放图像,数组等...

1815

扫描关注云+社区