安装并加载必须的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$Disease
exprSet=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)

boxplot

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

histogram

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

boxplot

p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
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)