GEO数据库视频学习笔记（了解你的表达矩阵）

继续上一篇笔记，这一个视频讲的是对表达矩阵进行各种的探索，所以视频的名称就叫“了解你的表达矩阵”，视频地址是：https://www.bilibili.com/video/BV1is411H7Hq?p=6

经过过滤探针（GEO数据库视频学习笔记（ID转换）），我们需要把ExprSet表达矩阵的行名（探针ID）换成基因名（从上一篇得到的ids表里提取）：

#现在的表达矩阵是18821个基因，6个样品
> dim(exprSet)
[1] 18821     6
#将ids表里的行匹配到exprSet里
> ids = ids[match(rownames(exprSet),ids$probe_id),]
> dim(ids) #现在ids里的基因数和表达矩阵一致了
[1] 18821     2
> rownames(exprSet) = ids$symbol
> head(exprSet)
          GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
LINC01128    8.75126    8.61650    8.81149    8.32067    8.41445    8.45208 
SAMD11       8.39069    8.52617    8.43338    9.17284    9.10216    9.14120
KLHL17       8.20228    8.30886    8.18518    8.13322    8.06453    8.15884
PLEKHN1      8.41004    8.37679    8.27521    8.34524    8.35557    8.44409
ISG15        7.72204    7.74572    7.78022    7.72308    7.53797    7.73401
AGRN         9.19237    9.10929    9.03668    9.94821    9.96994    9.99839

到此，表达矩阵就算是完成了。

（一）查看内参表达情况

下面可以看一下内参基因的表达情况，比如GAPDH:

> exprSet[rownames(exprSet)=="GAPDH",]
GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620 
   14.3187    14.3622    14.3638    14.4085    14.3569    14.3229 

可以看到大概是14左右，如何知道“14”是算表达高还是低呢？你可以对每一列的样品的基因做一个box看一下：

> boxplot(exprSet)

可以看到大部分样品的表达值在8左右，所以14算是很高的表达量了。如果你的管家基因的表达量低于平均水平，那么有可能你的探针ID和基因名转换的时候搞错了。

（二）查看表达矩阵分布图

接下来看一下表达矩阵的一些分布图：

> library(ggplot2) 
> library(reshape2)
#把表达矩阵转化成长list
> exprSet_L = melt(exprSet)
> colnames(exprSet_L) = c('group','sample','value')
#这里group的情况可以去GEO网页上显示的来命名，这里6个样品，分别是对照和药物处理，各3个重复，这里我只区分对照和处理，就不标明1,2,3了
> group_list = c(rep('control',3),rep('Vemurafenib',3))
#把长列表的表达矩阵里group一列命名成具体的样品名
> exprSet_L$group = rep(group_list,each = nrow(exprSet))
> head(exprSet_L)
    group     sample   value
1 control GSM1052615 8.75126
2 control GSM1052615 8.39069
3 control GSM1052615 8.20228
4 control GSM1052615 8.41004
5 control GSM1052615 7.72204
6 control GSM1052615 9.19237
> p = ggplot(exprSet_L,aes(x=sample,y= value,fill=group))+geom_boxplot()
> print(p)

这里要注意：如果你上面的图画出来，每一个样品中间那条线基本在同一水平线上，这些样品就可以在一起进行比较。如果相差很多，是不能进行比较的。说明不是一类的东西，或者说存在批次效应。

你还可以怎么获得group_list?

如果你用的是GEOquery包下载的表达矩阵：

gse <- getGEO('gse42872',GSEMatrix = TRUE, AnnotGPL = FALSE)
b = gse[[1]] #b是表达矩阵
tmp=pData(b)
View(tmp)

里面title这一列就是你的分组信息

也可以用小提琴图表示：

> p = ggplot(exprSet_L,aes(x=sample,y= value,fill=group))+geom_violin()
> print(p)
#图就不放了

（三）查看离群值

你还可以看这些样品里是否有离群值：

#看是否基因表达分布都差不多
> p = ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample,nrow =4)
> print(p)

（四）hclust图

> hc = hclust(dist(t(exprSet)))
> plot(hc)

现在你可以根据这个hclust图，看看是否这个图的分类和GEO网站上真正的样品分类是不是一样的，是不是对照在一组，处理在一组，如果不是，那肯定是某一步除了问题。这里你会发现，样品名还是GSM***，看着很乱，也不好一眼看出分类是不是有问题，那么可以把表达矩阵的样品名改一下：

> colnames(exprSet)= paste(group_list,1:3,sep='')
> colnames(exprSet)
[1] "control1"     "control2"     "control3"     "Vemurafenib1"
[5] "Vemurafenib2" "Vemurafenib3"
> hc = hclust(dist(t(exprSet)))
> plot(hc)

这样一眼就能看出你的分类有没有问题了。

（五）PCA

> 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')

可以看到这个对照和处理的样品分的还是很开的

参考资料：
生信人的20个R语言习题的答案