以CLL數(shù)據(jù)包為例:
1爹梁、提取所需數(shù)據(jù)
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex #直接查看數(shù)據(jù)集(S4對象)中的信息
exprSet=exprs(sCLLex)
##sCLLex是依賴于CLL這個package的一個對象
pdata=pData(sCLLex)
group_list=as.character(pdata[,2]) #pdata[,2]提取后為因子,需轉(zhuǎn)為字符串
exprSet[1:5,1:5] #當基因集較大時念链,查看部分數(shù)據(jù)可減少內(nèi)存消耗
gpl = sCLLex@annotation #獲取平臺信息及注釋包
2、安裝注釋包
#BiocManager::install('hgu95av2.db') #安裝注釋包
library(hgu95av2.db)
ls("package:hgu95av2.db") #查看注釋包的內(nèi)容
ids=toTable(hgu95av2SYMBOL) #toTable可將Bimap對象轉(zhuǎn)為數(shù)據(jù)框結(jié)構(gòu)
save(ids,exprSet,pdata,file = 'input.Rdata')
3掂墓、探索探針與基因symbol的對應(yīng)關(guān)系
length(unique(ids$symbol)) #查看共有多少個symbol,一般有多個探針對應(yīng)一個symbol
tail(sort(table(ids$symbol))) ##查看基因?qū)?yīng)最多多少個探針
table(sort(table(ids$symbol))) #顯示總的symbol對應(yīng)的探針數(shù)量并table
plot(table(sort(table(ids$symbol))))
table(rownames(exprSet) %in% ids$probe_id) #注意兩個數(shù)據(jù)集的先后順序跨嘉,%in%表示前者中的變量與后者一一比較,返回布爾值
4祠乃、對數(shù)據(jù)進行篩選,去重:
exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,] #取得與注釋包交集后的表達矩陣
dim(exprSet)
ids=ids[match(rownames(exprSet),ids$probe_id),] #match相當于%in%亮瓷,可以匹配返回邏輯值,也可以將前者數(shù)據(jù)的順序匹配給后者,使表達矩陣中的probe_id順序與ids中的一致
identical(ids$probe_id,rownames(exprSet)) #identical判斷兩個對象是否一致
dat=exprSet
ids$median=apply(dat,1,median) #ids新建median這一列嘱支,列名為median挣饥,同時對dat這個矩陣按行操作,取每一行的中位數(shù)亮靴,將結(jié)果給到median這一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#先對ids$symbol排序于置,在此基礎(chǔ)上按照ids$median中位數(shù)從大到小排列的順序排序,將對應(yīng)的行賦值為一個新的ids
ids=ids[!duplicated(ids$symbol),]#將symbol這一列取取出重復(fù)項八毯,'!'為否,即取出不重復(fù)的項话速,去除重復(fù)的gene ,保留每個基因最大表達量結(jié)果
dat=dat[ids$probe_id,] #新的ids取出probe_id這一列泊交,將dat按照取出的這一列中的每一行組成一個新的dat
rownames(dat)=ids$symbol#把ids的symbol這一列中的每一行給dat作為dat的行名
dat[1:4,1:4] #保留每個基因ID第一次出現(xiàn)的信息
dim(dat)
5、選擇對應(yīng)的數(shù)據(jù)作圖(基礎(chǔ)作圖):
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() 調(diào)整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") #在箱線圖框內(nèi)加均值的點
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()) #將圖中的全部文本字體設(shè)為粗體,x軸標簽設(shè)為30度傾斜铅乡,將坐標軸的標簽文字關(guān)閉
print(p)