繼續(xù)上一篇筆記，這一個視頻講的是對表達矩陣進行各種的探索，所以視頻的名稱就叫“了解你的表達矩陣”，視頻地址是：https://www.bilibili.com/video/BV1is411H7Hq?p=6
經(jīng)過過濾探針（GEO數(shù)據(jù)庫視頻學習筆記（ID轉(zhuǎn)換）），我們需要把ExprSet表達矩陣的行名（探針I(yè)D）換成基因名（從上一篇得到的ids表里提?。?/p>
#現(xiàn)在的表達矩陣是18821個基因，6個樣品
> dim(exprSet)
[1] 18821 6
#將ids表里的行匹配到exprSet里
> ids = ids[match(rownames(exprSet),ids$probe_id),]
> dim(ids) #現(xiàn)在ids里的基因數(shù)和表達矩陣一致了
[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
到此，表達矩陣就算是完成了。
（一）查看內(nèi)參表達情況
下面可以看一下內(nèi)參基因的表達情況，比如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算是很高的表達量了。如果你的管家基因的表達量低于平均水平，那么有可能你的探針I(yè)D和基因名轉(zhuǎn)換的時候搞錯了。
（二）查看表達矩陣分布圖
接下來看一下表達矩陣的一些分布圖：
> library(ggplot2)
> library(reshape2)
#把表達矩陣轉(zhuǎn)化成長list
> exprSet_L = melt(exprSet)
> colnames(exprSet_L) = c('group','sample','value')
#這里group的情況可以去GEO網(wǎng)頁上顯示的來命名，這里6個樣品，分別是對照和藥物處理，各3個重復，這里我只區(qū)分對照和處理，就不標明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)
這里要注意：如果你上面的圖畫出來，每一個樣品中間那條線基本在同一水平線上，這些樣品就可以在一起進行比較。如果相差很多，是不能進行比較的。說明不是一類的東西，或者說存在批次效應(yīng)。
你還可以怎么獲得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)
現(xiàn)在你可以根據(jù)這個hclust圖，看看是否這個圖的分類和GEO網(wǎng)站上真正的樣品分類是不是一樣的，是不是對照在一組，處理在一組，如果不是，那肯定是某一步除了問題。這里你會發(fā)現(xiàn)，樣品名還是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語言習題的答案
