表达矩阵可视化大全

2024-07-31 00:40:35

貌代码被折叠了，大家需要阅读原文才能复制粘贴我代码在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.874563group_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': ## ## lowessheatmap.2(choose_matrix)

library(pheatmap)pheatmap(choose_matrix)

DEG && volcano plot

library(limma)## ## Attaching package: 'limma'## The following object is masked from 'package:BiocGenerics': ## ## plotMAdesign=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画图是可以切换主题的

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