多分组芯片数据的差异分析-升级版

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

一样的代码是不？对，这个包是我写的的，就很丝滑。