1.背景知识
今天的目的是为了从一个5分组数据中筛选出两个分组，去做后续分析。扩展了一下，5个分组的数据其实也可以直接一步到位作差异分析的。
数据编号：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.下载和整理数据
geo_download函数可以完成GEO芯片数据的下载、匹配和整理，一步到位。
#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函数下载下来的数据包含表达矩阵、临床信息表格、GPL编号信息。
3.提取表达矩阵和分组信息
这个表达矩阵需要取个log。分组信息需要在pd表格里面找，对于这个数据，直接在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"
去除里面的数字就是合格的分组信息向量了。
Group = str_remove_all(geo$pd$title,"\\d")
table(Group)
## Group
## aHCC CL eHCC IL NL
## 13 10 13 10 10
可以看到5个分组及每个分组的样本数量。
4.提取自己需要的两个分组对应的数据
需要神奇符号%in%，x %in% y表示对x的每个元素进行判断，判断是否在y中存在，存在返回TRUE，不存在返回FALSE，生成了一个逻辑值向量。
用这个逻辑值向量对表达矩阵和分组信息分别取子集就可以啦！
因为表达矩阵的每一列和pd的每一行对应的是相同的样本，所以按照相同的条件取子集就实现了对样本的筛选。
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函数，方便你我他，把差异分析结果表格、差异基因和常见的几张图一起输出啦。
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
一样的代码是不？对，这个包是我写的的，就很丝滑。
