1.背景知識
今天的目的是為了從一個5分組數(shù)據(jù)中篩選出兩個分組，去做后續(xù)分析。擴展了一下，5個分組的數(shù)據(jù)其實也可以直接一步到位作差異分析的。
數(shù)據(jù)編號：GSE54238
實驗設計： In the study, we profiled lncRNA/mRNA expression in 10 normal livers (NL), 10 chronic inflammatory livers (IL), 10 cirrhotic livers (CL), 13 early HCC (eHCC) and 13 advanced HCC (aHCC) samples.
2.下載和整理數(shù)據(jù)
geo_download函數(shù)可以完成GEO芯片數(shù)據(jù)的下載、匹配和整理，一步到位。
#install.packages("tinyarray")
library(tinyarray)
library(tidyverse)
geo = geo_download("GSE54238")
str(geo,max.level = 1)
## List of 3
## $ exp: num [1:23845, 1:56] 2299.2 25939.9 25.3 236.5 345.4 ...
## ..- attr(*, "dimnames")=List of 2
## $ pd :'data.frame': 56 obs. of 10 variables:
## $ gpl: chr "GPL16955"
如上，geo_download函數(shù)下載下來的數(shù)據(jù)包含表達矩陣、臨床信息表格、GPL編號信息。
3.提取表達矩陣和分組信息
這個表達矩陣需要取個log。分組信息需要在pd表格里面找，對于這個數(shù)據(jù)，直接在pd的第一列(title)里面。
exp = log2(geo$exp+1)
geo$pd$title
## [1] "NL1" "NL2" "NL3" "NL4" "NL5" "NL6"
## [7] "NL7" "NL8" "NL9" "NL10" "IL1" "IL2"
## [13] "IL3" "IL4" "IL5" "IL6" "IL7" "IL8"
## [19] "IL9" "IL10" "CL1" "CL2" "CL3" "CL4"
## [25] "CL5" "CL6" "CL7" "CL8" "CL9" "CL10"
## [31] "eHCC1" "eHCC2" "eHCC3" "eHCC4" "eHCC5" "eHCC6"
## [37] "eHCC7" "eHCC8" "eHCC9" "eHCC10" "eHCC11" "eHCC12"
## [43] "eHCC13" "aHCC1" "aHCC2" "aHCC3" "aHCC4" "aHCC5"
## [49] "aHCC6" "aHCC7" "aHCC8" "aHCC9" "aHCC10" "aHCC11"
## [55] "aHCC12" "aHCC13"
去除里面的數(shù)字就是合格的分組信息向量了。
Group = str_remove_all(geo$pd$title,"\\d")
table(Group)
## Group
## aHCC CL eHCC IL NL
## 13 10 13 10 10
可以看到5個分組及每個分組的樣本數(shù)量。
4.提取自己需要的兩個分組對應的數(shù)據(jù)
需要神奇符號%in%，x %in% y表示對x的每個元素進行判斷，判斷是否在y中存在，存在返回TRUE，不存在返回FALSE，生成了一個邏輯值向量。
用這個邏輯值向量對表達矩陣和分組信息分別取子集就可以啦！
因為表達矩陣的每一列和pd的每一行對應的是相同的樣本，所以按照相同的條件取子集就實現(xiàn)了對樣本的篩選。
k = Group %in% c("NL", "CL");table(k)
## k
## FALSE TRUE
## 36 20
exp1 = exp[,k]
Group1 = factor(Group[k],levels = c("NL", "CL"))
ncol(exp1)
## [1] 20
length(Group1)
## [1] 20
5.探針注釋
library(GEOquery)
a = getGEO(geo$gpl)
b = a@dataTable@table[c(1,6)]
colnames(b) = c("probe_id","symbol")
b = b[b$symbol!="",]
head(b)
## probe_id symbol
## 1 NM_000015 NAT2
## 2 NM_000016 ACADM
## 3 NM_000017 ACADS
## 4 NM_000018 ACADVL
## 5 NM_000020 ACVRL1
## 6 NM_000023 SGCA
6.二分組差異分析一步到位
get_deg_all函數(shù)，方便你我他，把差異分析結果表格、差異基因和常見的幾張圖一起輸出啦。
dcp1 = get_deg_all(exp1,Group1,b)
str(dcp1,max.level = 1)
## List of 3
## $ deg :'data.frame': 10035 obs. of 10 variables:
## $ cgs :List of 3
## $ plots:A patchwork composed of 3 patches
## - Autotagging is turned off
## - Guides are collected
##
## Layout:
## 3 patch areas, spanning 3 columns and 1 rows
##
## t l b r
## 1: 1 1 1 1
## 2: 1 2 1 2
## 3: 1 3 1 3
head(dcp1$deg)
## logFC AveExpr t P.Value adj.P.Val
## 1 -2.944973 9.188322 -9.353001 2.941332e-09 7.013607e-05
## 2 2.471296 9.458068 7.924233 5.449593e-08 6.454766e-04
## 3 -2.667746 8.150188 -7.161227 2.894538e-07 9.860037e-04
## 4 -1.660002 8.636170 -6.999874 4.161179e-07 9.936383e-04
## 5 -2.045595 10.231269 -6.658487 9.068426e-07 1.801972e-03
## 6 -1.937557 8.059010 -6.464551 1.420836e-06 1.903384e-03
## B probe_id symbol change ENTREZID
## 1 10.934721 NM_021098 CACNA1H down 8912
## 2 8.341558 NM_006074 TRIM22 up 10346
## 3 6.825299 NM_004213 SLC28A1 down 9154
## 4 6.493086 NM_001039199 TTPAL down 79183
## 5 5.777318 NM_001018011 ZBTB16 down 7704
## 6 5.363146 NM_173483 CYP4F22 down 126410
dcp1$plots
你會看到熱圖的聚類與分組不匹配是吧。這不是個錯誤哦。詳見：分組聚類的熱圖
7.多分組的差異分析也可以搞定
Group = factor(Group,levels = c("NL", "IL", "eHCC", "CL", "aHCC"))
dcp2 = get_deg_all(exp,Group,b)
str(dcp2,max.level = 1)
## List of 3
## $ deg :List of 4
## $ cgs :List of 4
## $ plots:A patchwork composed of 4 patches
## - Autotagging is turned off
## - Guides are collected
##
## Layout:
## 4 patch areas, spanning 2 columns and 2 rows
##
## t l b r
## 1: 1 1 1 1
## 2: 2 1 2 1
## 3: 1 2 1 2
## 4: 2 2 2 2
dcp2$plots
一樣的代碼是不？對，這個包是我寫的的，就很絲滑。
