【生信技能树】GEO数据库挖掘 P6 5 了解你的表达矩阵
# 检查matrix表达量是否合适
# 载入内参基因的表达量,看是否与相应的基因名相统一
# GAPDH 以及 ACTB(b-actin的基因名)是常用的内参
# 转换exprset的rownames为基因名
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
# 检查matrix表达量是否合适
# 载入内参基因的表达量,看是否与相应的基因名相统一
# GAPDH 以及 ACTB(b-actin的基因名)是常用的内参
# 转换exprset的rownames为基因名
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
后续检查样本是否正常表达一些基因。
【方法1】用GAPDH和ACTB基因作为阳性参照,看是否表达量正常。(应保守的高表达)
【方法2】用melt将样本变成纵列后,重命名以及予以分组信息,使用ggplot2看基因表达量的分布。
# 检查1,常见基因表达量,确定管家基因是否弄错
exprSet['GAPDH',]
exprSet['ACTB',]
# 检查2 检验基因分布图
if(T){
library(reshape2)
exprSet_L=melt(exprSet)
colnames(exprSet_L)=c('probe','sample','value')
# group_list获取每个样本的分组信息
group_list=c('Control','Control','Control','Vemurafenib','Vemurafenib','Vemurafenib')
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)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
print(p)
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)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
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)
# 如果来自不同样本或基线不匹配,需使用sv包中的combine函数来矫正
# 检查样本是否有问题,是否有异常值或不一致的分布
}
【方法3】用hclust包看以下样本间的关系,是否相同处理聚类到了一起。
grouplist可以根据网址上的样本信息予以确定,也可以根据table确定。
# 检查3 利用hclust包进行聚类,简单看一下进化关系,是否每个样本都类似
if(T){
# pdata=pData(data_information)
group_list=as.character(pdata[,2])
group_list
dim(exprSet)
exprSet[1:5,1:5]
## hclust
colnames(exprSet)=paste(group_list,1:6,sep='')
# Define nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
cex = 0.7, col = "blue")
hc=hclust(dist(t(exprSet)))
par(mar=c(5,5,5,10))
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)
}
【方法4】PCA主成分分析看一下样本情况。
# 检查4 利用PCA来看样本分布是否一致
if(F){
# BiocManager::install('ggfortify')
library(ggfortify)
df=as.data.frame(t(exprSet))
df$group=group_list
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")
df=as.data.frame(t(exprSet))
dat.pca <- PCA(df, graph = FALSE)#现在dat最后一列是group_list,需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = group_list, # color by groups
# palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
}