以CLL数据包为例:
1、提取所需数据
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex #直接查看数据集(S4对象)中的信息
exprSet=exprs(sCLLex)
##sCLLex是依赖于CLL这个package的一个对象
pdata=pData(sCLLex)
group_list=as.character(pdata[,2]) #pdata[,2]提取后为因子,需转为字符串
exprSet[1:5,1:5] #当基因集较大时,查看部分数据可减少内存消耗
gpl = sCLLex@annotation #获取平台信息及注释包
2、安装注释包
#BiocManager::install('hgu95av2.db') #安装注释包
library(hgu95av2.db)
ls("package:hgu95av2.db") #查看注释包的内容
ids=toTable(hgu95av2SYMBOL) #toTable可将Bimap对象转为数据框结构
save(ids,exprSet,pdata,file = 'input.Rdata')
3、探索探针与基因symbol的对应关系
length(unique(ids$symbol)) #查看共有多少个symbol,一般有多个探针对应一个symbol
tail(sort(table(ids$symbol))) ##查看基因对应最多多少个探针
table(sort(table(ids$symbol))) #显示总的symbol对应的探针数量并table
plot(table(sort(table(ids$symbol))))
table(rownames(exprSet) %in% ids$probe_id) #注意两个数据集的先后顺序,%in%表示前者中的变量与后者一一比较,返回布尔值
4、对数据进行筛选,去重:
exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,] #取得与注释包交集后的表达矩阵
dim(exprSet)
ids=ids[match(rownames(exprSet),ids$probe_id),] #match相当于%in%,可以匹配返回逻辑值,也可以将前者数据的顺序匹配给后者,使表达矩阵中的probe_id顺序与ids中的一致
identical(ids$probe_id,rownames(exprSet)) #identical判断两个对象是否一致
dat=exprSet
ids$median=apply(dat,1,median) #ids新建median这一列,列名为median,同时对dat这个矩阵按行操作,取每一行的中位数,将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#先对ids$symbol排序,在此基础上按照ids$median中位数从大到小排列的顺序排序,将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项,'!'为否,即取出不重复的项,去除重复的gene ,保留每个基因最大表达量结果
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列,将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4] #保留每个基因ID第一次出现的信息
dim(dat)
5、选择对应的数据作图(基础作图):
exprSet=dat
#第一个样本的所有基因的表达量作图
boxplot(exprSet[,1]) #箱线图
p<-hist(exprSet[,1]) #直方图
barplot(exprSet[,1]) #条形图
plot(density(exprSet[,1])) #单独核密度图
polygon(density(exprSet[,1]),col="red", border="blue") #多边形绘制
6、ggplot2作图:
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)
##绘图
### ggplot2
library(ggplot2)
p=ggplot(exprSet_L,
aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p) ##箱线图,fill映射到group
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()+theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) #小提琴图,theme(axis.text.x = element_text() 调整x轴刻度文字的角度
print(p)
p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4) #直方图,用sample对值进行分面facet_wrap
print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4) #密度图,facet_wrap分面
print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()
print(p)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
p=p+stat_summary(fun="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()) #将图中的全部文本字体设为粗体,x轴标签设为30度倾斜,将坐标轴的标签文字关闭
print(p)