GEO数据库挖掘—生信技能树B站视频

以下是B站生信技能树GEO数据库挖掘的课程笔记

主要内容及学习目的：

  • 介绍GEO数据库：了解数据存放位置；
  • 介绍GSE项目的3种下载方式；
  • 介绍ID转换：使用R语言技巧实现基因ID之间的转换，我们下载的数据通常使用的是不同的芯片探针，它们有不同的探针ID号我们需要把它转化成ENTREZID或SYMBOLID才能被大众认知；
  • 介绍表达矩阵的相关可视化及归一化：从GEO数据库下载的是作者处理好的矩阵，我们需要会判别它是否符合要求，并学会提取分组信息；
  • 比较各组基因的表达量得到差异表达基因list或感兴趣基因集；
  • 得到差异表达基因list后做富集分析；
  • 用GSEA软件做一些图；

通过阅读文章提炼GEO数据库挖掘的脉络：选择GSE ---> 得到表达矩阵 ---> control VS treatment 进行差异分析 ---> 得到差异表达基因list ---> 5大数据库的注释 ---> PPI等网络

GEO数据库挖掘分析思路

接下来我们按照上面的分析思路，一步一步进行讲解

1.了解GEO数据库，找到文章的GSE编号

参考文案：解读GEO数据存放规律及下载，一文就够

任何一篇GEO数据挖掘文章，都可以找到它的GSE编号，找到后我们把网址最后的GSE编号修改一下，直接去网页粘贴并转到就能看到该编号在GEO数据库的详细页面：
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42872

找到文章中的GSE编号

我们看下GEO数据库的主页
GEO数据库主页

我们只需要知道这三个概念就好：

  • GEO Platform (GPL)
  • GEO Sample (GSM)

  • GEO Series (GSE)


理解起来也很容易。一篇文章可以有一个或者多个GSE数据集，一个GSE数据集里面可以有一个或者多个GSM样本。而每个数据集都有着自己对应的芯片平台，就是GPL。


点开后向下滑动找到GPL的表格信息

把GPL表格信息向右滑动，找到gene_assignment那列，把//做为切割符，取出第二个字符就是真正的基因名，这时探针就和基因完美匹配啦~
得到真正的基因名

知道如何找到任何一篇文章的数据存放位置，接下来就要下载数据进行分析了。

2.下载数据

下载数据的3种方式：

一. 直接下载rawdata —— 不推荐使用

二. 从网页上直接下载表达矩阵 ---> 读入R里

表达矩阵

表达矩阵下载到本地后要读入到R里：

a = read.table(file="./GSE42872_series_matrix.txt.gz",
               header = T,sep = "\t",quote = "",fill = T,
               comment.char = "!")

在读入下载好的表达矩阵时，为什么要加那么多参数才能下载成功？
我们首先需要在电脑上解压并打开文本文件，根据文件的样子选择参数

GSE42872_series_matrix.txt.gz解压后

我们看a的rowname是行号，没有意义的，需要转成探针ID号即a的第一列：rownames(a)= a[,1]
a的第一列就是ID号

> head(a)
  X.ID_REF. X.GSM1052615. X.GSM1052616. X.GSM1052617. X.GSM1052618. X.GSM1052619. X.GSM1052620.
1   7892501       7.24559       6.80686       7.73301       6.18961       7.05335       7.20371
2   7892502       6.82711       6.70157       7.02471       6.20493       6.76554       6.24252
3   7892503       4.39977       4.50781       4.88250       4.36295       4.18137       4.73492
4   7892504       9.48025       9.67952       9.63074       9.69200       9.91324       9.65897
5   7892505       4.54734       4.45247       5.11753       4.87307       5.15505       3.99340
6   7892506       6.80701       6.90597       6.72472       6.77028       6.77058       6.77685
> rownames(a)
   [1] "1"    "2"    "3"    "4"    "5"    "6"    "7"    "8"    "9"    "10"   "11"   "12"   "13"  
> rownames(a)= a[,1]
> head(a)
        X.ID_REF. X.GSM1052615. X.GSM1052616. X.GSM1052617. X.GSM1052618. X.GSM1052619. X.GSM1052620.
7892501   7892501       7.24559       6.80686       7.73301       6.18961       7.05335       7.20371
7892502   7892502       6.82711       6.70157       7.02471       6.20493       6.76554       6.24252
7892503   7892503       4.39977       4.50781       4.88250       4.36295       4.18137       4.73492
7892504   7892504       9.48025       9.67952       9.63074       9.69200       9.91324       9.65897
7892505   7892505       4.54734       4.45247       5.11753       4.87307       5.15505       3.99340
7892506   7892506       6.80701       6.90597       6.72472       6.77028       6.77058       6.77685

此时，a的列名就是探针ID号了。但是现在还是不符合预期，我们还要把RefSeq ID那一列去掉，也就是去掉此时的第一列：a = a[,-1]

> a = a[,-1]
> head(a)
        X.GSM1052615. X.GSM1052616. X.GSM1052617. X.GSM1052618. X.GSM1052619. X.GSM1052620.
7892501       7.24559       6.80686       7.73301       6.18961       7.05335       7.20371
7892502       6.82711       6.70157       7.02471       6.20493       6.76554       6.24252
7892503       4.39977       4.50781       4.88250       4.36295       4.18137       4.73492
7892504       9.48025       9.67952       9.63074       9.69200       9.91324       9.65897
7892505       4.54734       4.45247       5.11753       4.87307       5.15505       3.99340
7892506       6.80701       6.90597       6.72472       6.77028       6.77058       6.77685

这就是由样本和探针组成的表达矩阵

三. 在R里使用GSE号和GEOquery包从GEO数据库上直接下载——最推荐使用下载方式

library(GEOquery)
eSet <- getGEO("GSE42872", 
               destdir = '.',  #下载在当前目录
               getGPL = F) #平台信息不要

使用以上代码就可以将GSE42872数据下载到R里当前工作目录并赋值给eSet，下载完成后要注意检查数据文件的完整性——看我们下载的数据大小是否大于等于官网上给的大小。如果我们下载的数据内存大于官网上的那没事儿，如果小于官网上的那下载的数据就不完整。

2.1得到表达矩阵

我们用方法2下载的a和用方法3下载的eSet都是GSE42872数据，但它们是不一样的：

> class(a)
[1] "data.frame"
> class(eSet)
[1] "list"

我们可以看到a是一个数据框，eSet是一个列表这里我们称它为对象。
得到eSet对象里包含着各种各样的信息：表达矩阵、芯片如何设计的、样本如何分组 等等~
eSet是一个大列表，我们需要从中提取出表达矩阵，才能进行后续的操作。
为什么？因为一个GSE号里面对应多种芯片平台数据，我们使用GSE号下载数据就会把所有芯片平台的数据整合到一个list里面，每个list里的元素存放一个平台的表达矩阵。我们的数据只有一个平台所以eSet列表里就只有一个元素：


使用列表取子集的方法提取eSet列表里的第一个元素：eSet[[1]]；并使用exprs函数把它转化成矩阵：exp <- exprs(eSet[[1]])

> eSet[[1]]
ExpressionSet (storageMode: lockedEnvironment)
assayData: 33297 features, 6 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: GSM1052615 GSM1052616 ... GSM1052620 (6 total)
  varLabels: title geo_accession ... cell type:ch1 (34 total)
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
  pubMedIds: 24469106 
Annotation: GPL6244 
> exp <- exprs(eSet[[1]])
> exp[1:4,1:4]
        GSM1052615 GSM1052616 GSM1052617 GSM1052618
7892501    7.24559    6.80686    7.73301    6.18961
7892502    6.82711    6.70157    7.02471    6.20493
7892503    4.39977    4.50781    4.88250    4.36295
7892504    9.48025    9.67952    9.63074    9.69200

这时我们再看由方法2得到的表达矩阵a，和由方法3得到的表达矩阵exp是一模一样的：

> head(a)
        X.GSM1052615. X.GSM1052616. X.GSM1052617. X.GSM1052618. X.GSM1052619. X.GSM1052620.
7892501       7.24559       6.80686       7.73301       6.18961       7.05335       7.20371
7892502       6.82711       6.70157       7.02471       6.20493       6.76554       6.24252
7892503       4.39977       4.50781       4.88250       4.36295       4.18137       4.73492
7892504       9.48025       9.67952       9.63074       9.69200       9.91324       9.65897
7892505       4.54734       4.45247       5.11753       4.87307       5.15505       3.99340
7892506       6.80701       6.90597       6.72472       6.77028       6.77058       6.77685
> head(exp)
        GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
7892501    7.24559    6.80686    7.73301    6.18961    7.05335    7.20371
7892502    6.82711    6.70157    7.02471    6.20493    6.76554    6.24252
7892503    4.39977    4.50781    4.88250    4.36295    4.18137    4.73492
7892504    9.48025    9.67952    9.63074    9.69200    9.91324    9.65897
7892505    4.54734    4.45247    5.11753    4.87307    5.15505    3.99340
7892506    6.80701    6.90597    6.72472    6.77028    6.77058    6.77685

基因ID之间的转换

entrez ID <---probe_id--->symbol ID
分两步走：

  • 过滤probe_id，得到每个基因所对应的唯一的probe_id
  • 得到probe_id与symbol ID这件的转换关系

使用R语言技巧实现基因ID之间的转换，我们下载的数据通常使用的是不同的芯片探针，它们有不同的探针ID（probe_id）我们需要把它转化成entrez ID或symbol ID才能被大众认知；

所以接下来，我们学习怎样将探针ID（probe_id）转换成entrez ID和symbol ID
ID转换的第一步必须要加载特定的R包，下载哪个包，需要根据GPL来定

> eSet[[1]]
ExpressionSet (storageMode: lockedEnvironment)
assayData: 33297 features, 6 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: GSM1052615 GSM1052616 ... GSM1052620 (6 total)
  varLabels: title geo_accession ... cell type:ch1 (34 total)
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
  pubMedIds: 24469106 
Annotation: GPL6244

我们看到GPL号是：GPL6244

> eSet[[1]]@annotation
[1] "GPL6244"

用R获取芯片探针与基因的对应关系三部曲-bioconductor里搜索GPL6244所对应的R包，发现是hugene10sttranscriptcluster

http://www.bio-info-trainee.com/1399.html

所以我们加载此包：
library(hugene10sttranscriptcluster.db) 注意加上后缀 .db

> library(hugene10sttranscriptcluster.db)
> ls("package:hugene10sttranscriptcluster.db")
 [1] "hugene10sttranscriptcluster"              "hugene10sttranscriptcluster.db"          
 [3] "hugene10sttranscriptcluster_dbconn"       "hugene10sttranscriptcluster_dbfile"      
 [5] "hugene10sttranscriptcluster_dbInfo"       "hugene10sttranscriptcluster_dbschema"    
 [7] "hugene10sttranscriptclusterACCNUM"        "hugene10sttranscriptclusterALIAS2PROBE"  
 [9] "hugene10sttranscriptclusterCHR"           "hugene10sttranscriptclusterCHRLENGTHS"   
[11] "hugene10sttranscriptclusterCHRLOC"        "hugene10sttranscriptclusterCHRLOCEND"    
[13] "hugene10sttranscriptclusterENSEMBL"       "hugene10sttranscriptclusterENSEMBL2PROBE"
[15] "hugene10sttranscriptclusterENTREZID"      "hugene10sttranscriptclusterENZYME"       
[17] "hugene10sttranscriptclusterENZYME2PROBE"  "hugene10sttranscriptclusterGENENAME"     
[19] "hugene10sttranscriptclusterGO"            "hugene10sttranscriptclusterGO2ALLPROBES" 
[21] "hugene10sttranscriptclusterGO2PROBE"      "hugene10sttranscriptclusterMAP"          
[23] "hugene10sttranscriptclusterMAPCOUNTS"     "hugene10sttranscriptclusterOMIM"         
[25] "hugene10sttranscriptclusterORGANISM"      "hugene10sttranscriptclusterORGPKG"       
[27] "hugene10sttranscriptclusterPATH"          "hugene10sttranscriptclusterPATH2PROBE"   
[29] "hugene10sttranscriptclusterPFAM"          "hugene10sttranscriptclusterPMID"         
[31] "hugene10sttranscriptclusterPMID2PROBE"    "hugene10sttranscriptclusterPROSITE"      
[33] "hugene10sttranscriptclusterREFSEQ"        "hugene10sttranscriptclusterSYMBOL"       
[35] "hugene10sttranscriptclusterUNIGENE"       "hugene10sttranscriptclusterUNIPROT"

通过命令：
ls("package:hugene10sttranscriptcluster.db")
我们可以看到这个包里面有很多数据集，想要得到probe_id和symbol的对应关系要用hugene10sttranscriptclusterSYMBOL数据集，用toTable函数提取数据集里面的信息：

> ids=toTable(hugene10sttranscriptclusterSYMBOL)
> head(ids)
  probe_id    symbol
1  7896759 LINC01128
2  7896761    SAMD11
3  7896779    KLHL17
4  7896798   PLEKHN1
5  7896817     ISG15
6  7896822      AGRN

现在我们查看下一共多少个基因？一万八千多个基因

> length(unique(ids$symbol))
[1] 18834

unique函数是用来：Extract Unique Elements 去除重复的symbol只提取不同的元素；length函数统计去重之后还有多少个基因。

再查看每个基因对应多少个探针：

> tail(sort(table(ids$symbol)))

  RPL41  UBTFL1  CDK11B  UBE2D3    IGKC LRRFIP1 
      6       6       8       8      10      10

可以看到有的基因设计了10个探针或8个探针....
table() 函数可以生成频数统计表，这里就是统计每个基因symbol出现的次数然后将其表格化；sort()函数将symbol出现的频率从小到大进行排序；tail()取最后6个即出现频率最大的6个。

> table(sort(table(ids$symbol)))

    1     2     3     4     5     6     8    10 
18072   599   132    16     5     6     2     2

table一下我们可以看到，18072个基因设计了1个探针；599个基因设计了2个探针；132个基因设计了3个探针.....也就是说大部分的基因只设计了1个探针。

其实一般基因都会设计很多探针的，我们下载的表达矩阵是作者处理之后的，把许多不好的探针都过滤掉了，我们处理作者的数据要默认人家做的是对的，否则就要下载原始数据自己处理。

> table(rownames(exp) %in% ids$probe_id)

FALSE  TRUE 
13470 19827

发现有13470个探针没有对应的基因名；19827个探针有对应的基因名。
x %in% y表示 x 的元素在y中吗？然后返回逻辑值。rownames(exp)即表达矩阵exp的行名是文章数据中用到的所有探针ID（probe_id）；ids$probe_id是具有对应基因的所有探针。所以返回的TRUE就是文章数据中有对应基因的探针数。

现在我们对探针进行过滤，把没有对应基因名的探针过滤掉：

> exp = exp[rownames(exp) %in% ids$probe_id,]

过滤的本质就是矩阵取子集，如：matrix[2,]意思就是取矩阵matrix的第2行和所有的列。同理，我们这里exp[rownames(exp) %in% ids$probe_id,]就是取具有对应基因的所有探针（行），和所有的列。

对比一下过滤之前和过滤之后的探针数量：

> table(rownames(exp) %in% ids$probe_id)

FALSE  TRUE 
13470 19827 
> dim(exp)
[1] 33297     6
> exp = exp[rownames(exp) %in% ids$probe_id,]
> dim(exp)
[1] 19827     6

可以发现过滤之前有33297个探针，过滤之后就剩下19827个探针了。

然后，我们使用match函数把ids里的探针顺序改一下，使ids里探针顺序和我们表达矩阵的顺序完全一样：

ids=ids[match(rownames(exp),ids$probe_id),]

match()函数返回的是一个位置向量，该向量记录着第一个参数中每个元素在第二个参数中的位置。所以，此时ids里的探针顺序与表达矩阵exp的探针顺序一一对应：

> head(ids)
  probe_id    symbol
1  7896759 LINC01128
2  7896761    SAMD11
3  7896779    KLHL17
4  7896798   PLEKHN1
5  7896817     ISG15
6  7896822      AGRN
> head(exp)
        GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
7896759    8.75126    8.61650    8.81149    8.32067    8.41445    8.45208
7896761    8.39069    8.52617    8.43338    9.17284    9.10216    9.14120
7896779    8.20228    8.30886    8.18518    8.13322    8.06453    8.15884
7896798    8.41004    8.37679    8.27521    8.34524    8.35557    8.44409
7896817    7.72204    7.74572    7.78022    7.72308    7.53797    7.73401
7896822    9.19237    9.10929    9.03668    9.94821    9.96994    9.99839

既然已经完全对应上，我们就可以通过probe_id将表达矩阵exp进行分组，将同一个symbol所对应的多个探针分成不同的组，并对每组探针进行统计：计算每组中每行探针表达量的平均值（也就是每个探针在6个样本中表达量的均值rowMeans(x)），再取平均值最大的那个探针作为该symbol所对应的唯一探针，该组中的其它探针过滤掉，这样每个symbol就对应一个探针了，看下代码是如何操作的：

tmp = by(exp,
         ids$symbol,
         function(x) rownames(x)[which.max(rowMeans(x))])
probes = as.character(tmp)
dim(exp)
exp = exp[rownames(exp) %in% probes,] # 过滤有多个探针的基因
dim(exp)

是不是没有理解上面代码在干些什么？没关系，我们详细解释一下：
by()函数在这里发挥的功能就是将表达矩阵exp中的探针分组，同一个symbol所对应的多个探针分到一组，并对每组探针进行统计得到symbol所对应的唯一探针。所以tmp里放着by()函数的统计结果即每个symbol所对应的唯一探针IDprobe_id，用probes = as.character(tmp)将结果变身为纯字符型向量：

> head(tmp)
INDICES
     A1CF       A2M     A2ML1   A3GALT2    A4GALT     A4GNT 
"7933640" "7960947" "7953775" "7914643" "8076497" "8090955" 
> head(probes)
[1] "7933640" "7960947" "7953775" "7914643" "8076497" "8090955"
>

学习by()函数如何完成以上操作的。《R语言实战》这本书上是这样描述的

使用by()分组计算描述性统计量，它可以一次返回若干个统计量。格式为：
by(data, INDICES, FUN)
其中data是一个数据框或矩阵；INDICES是一个因子或因子组成的列表，定义了分组；FUN是任意函数。

简单一句话理解就是：by()函数就是根据因子将整个data分成几个小的data.frame，然后进行运算处理。
同理，我们这里：
by(exp, ids$symbol, function(x) rownames(x)[which.max(rowMeans(x))])
第二个参数ids$symbol定义了分组，将第一参数—exp表达矩阵分成了若干个小矩阵，每个小矩阵里存放着同一个symbol所对应的所有探针。第三个参数是我们自己定义的函数：计算每个小矩阵中每行探针表达量的平均值（也就是每个探针在6个样本中表达量的均值rowMeans(x)），再取平均值最大的那个探针作为该symbol所对应的唯一探针which.max(rowMeans(x))。
by()函数就可以返回每个分组里的统计结果，即每个symbol所对应的唯一探针IDprobe_id。
这时，探针ID和基因symbol就一一对应了，将表达矩阵探针ID即exp表达矩阵的行名（rownames(exp)）换为基因symbol：

rownames(exp)=ids[match(rownames(exp),ids$probe_id),2]
> head(exp)
          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

此时，我们已经将探针ID转化成基因symbol了。
由上面的介绍我们可以看到，在转换ID中最重要的是根据GPL平台号找到所对应的R注释包，可是如果找不到GPL平台对应的R注释包怎么办呢？
答：我们不用GEO号进行下载，而是下载平台信息(GPL)，从平台信息中选择我们想要的列：探针名、基因名....
GPL里面的信息量特别大，下载特别考验网速。

gpl <- getGEO('GPL6480', destdir = ".")
colnames(Table(gpl))
head(Table(gpl)[,c(1,6,7)]) #看gpl对象中哪一列是我们想要的取出来，发现1/6/7列是我们想要的
write.csv(Table(gpl)[,c(1,6,7)],"GPL6480.csv") #把我们想要的部分即探针名对应的基因名....存起来

获取分组信息—group_list

分组信息就是告诉我们哪些组是control；哪些组是tumor。
使用pData函数获取分组信息—group_list：

pd <- pData(eSet[[1]])

pData()函数可以得到每个样本的描述信息，一般来说数据框的第一列(title列)描述了哪些是control；哪些是treatment。

pd <- pData(eSet[[1]])

根据第一列所描述的信息我们自己创建分组信息group_list：
方法一：使用stringr函数

library(stringr)
group_list = ifelse(str_detect(pd$title,"Control")==TRUE,"contorl","treat")
group_list

stringr包用于字符串的处理，str_detect是该包里的函数，用来确定一个字符向量能否匹配一种模式。它返回一个与输入向量具有同样长度的逻辑向量：

> str_detect(pd$title,"Control")
[1]  TRUE  TRUE  TRUE FALSE FALSE FALSE

这里的输入向量是数据框pd的第一列pd$title内容，即由6个元素组成的字符型向量。str_detect()函数会自动判断Control，是否存在于pd$title向量的每一个元素中，存在返回TRUE，否则返回FALSE。
str_detect函数处理后我们再使用 ifelse生成符合要求的分组信息group_list

> group_list
[1] "contorl" "contorl" "contorl" "treat"   "treat"   "treat"

方法二：自己造一个
我们已经知道了前三个是control后三个是treatment，那就自己生成一个符合要求的分组信息：

group_list=c(rep("control",times=3),rep("treat",times=3))
group_list

3. 检查表达矩阵

得到表达矩阵就是描述了某个基因在某个样本的表达量。有了这个表达矩阵我们可以做后面的分析，第一步就是确定我们得到的表达矩阵是否正确：

  • 查看管家基因的表达量
  • 检测分组之间是否有差异：PCA图、热图和hclust图等等

3.1 检验常见基因的表达量

查看典型管家基因（如：GAPDH、ACTB）的表达量，如果表达量高于正常值，说明我们数据没问题。
此时表达矩阵exp的行名已经由探针ID转换成基因名了，所以我们使用exp['GAPDH',]来提取该基因在所有样品中的表达量。

> exp['GAPDH',]
GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620 
   14.3187    14.3622    14.3638    14.4085    14.3569    14.3229 
> exp['ACTB',]
GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620 
   13.8811    13.9002    13.8822    13.7773    13.6732    13.5363

我们可以看到我们数据中两个管家基因的表达量都偏高，符合预期。为什么知道它偏高呢？画一个整体样本所有基因的表达量的boxplot：boxplot(exp)


发现大部分基因的表达量都在8-9，而GAPDH、ACTB在13-14左右，所以是偏高的。
假如，我们发现管家基因表达量特别低，那我们就要思考是不是在提取表达矩阵的时候哪里出了问题：比如ID转换的时候转换错了等等....

3.2 看表达矩阵的分布图—画图看各个样本的表达量

使用ggplot2画各个样本表达量的boxplot图

# 准备画图所需数据exp_L
library(reshape2)
head(exp)
exp_L = melt(exp)
head(exp_L)
colnames(exp_L)=c('symbol','sample','value')
head(exp_L)

# 获得分组信息
library(stringr)
group_list = ifelse(str_detect(pd$title,"Control")==TRUE,"contorl","treat")
group_list
exp_L$group = rep(group_list,each=nrow(exp))
head(exp_L)

# ggplot2画图 
library(ggplot2)
p = ggplot(exp_L,
         aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)

##boxplot图精修版
p=ggplot(exp_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
p=p+theme_set(theme_set(theme_bw(base_size=20)))
p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
print(p)

作图函数看起来复杂我们拆开：

准备作图所需要的数据exp_L ---> 获得分组信息并加到exp_L中 ---> ggplot2作图

我们先理解一下 exp_L 数据

> head(exp_L)
     symbol     sample   value   group
1 LINC01128 GSM1052615 8.75126 contorl
2    SAMD11 GSM1052615 8.39069 contorl
3    KLHL17 GSM1052615 8.20228 contorl
4   PLEKHN1 GSM1052615 8.41004 contorl
5     ISG15 GSM1052615 7.72204 contorl
6      AGRN GSM1052615 9.19237 contorl

> table(exp_L[,2])
GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620 
     18834      18834      18834      18834      18834      18834 

> dim(exp_L)
[1] 113004      4

> 18834*6
[1] 113004

由以上代码我们可以看到exp_L矩阵是这样分布的：每个基因（18834个）在第一个样本GSM1052615中的value值，每个基因（18834个）在第二个样本中的value值....以此类推一共有6个样本。

难点攻克：如何得到这样的exp_L矩阵呢？？？使用reshape2包

reshape2包是一套重构和整合数据集的绝妙的万能工具。大致用法就是，需要首先将数据融合（melt），以使每一行都是唯一的标识符-变量组合。然后将数据重塑（cast）为你想要的任何形状。在重铸过程中，你可以使用任何函数对数据进行整合。
我们这里只用到这个包里的数据融合（melt）功能。
数据集的融合（melt）是将它重构为这样一种格式：每个测量变量（每个基因在每个样本中的表达量）独占一行，行中带有要唯一确定这个测量所需的标识符变量（基因symbol和样本sample）。注意，必须指定要唯一确定每个测量所需的变量（也就是说基因symbol和样本sample必须对应唯一的表达量），而表示测量变量名的变量将由程序为你自动创建（即表达量独占一行后程序会自动创建表达量所对应的symbol和sample）。
说成人话就是，以前exp矩阵是一个基因在6个样本中的表达量占一行，melt后就会将表达量独占一行。一个表达量的值需要有两个定语才能唯一指定，即这个表达量是哪个样本（sample）中的哪个基因（symbol）的。

各个样本表达量的boxplot

从图中可以看到两个分组control和treat基本在一条线上，这样的数据说明可以进行后续比较，如果不在一条线上说明有批次效应batch infect，需要用limma包内置函数normalizeBetweenArrays人工校正一下(Normalization)：

library(limma) 
exp = normalizeBetweenArrays(exp)

关于画样本表达量的分布图，除了上面介绍的boxplot，ggplot2还可以画别的，看情况使用就好，不同的图有不同的展现方式但都在展现同一个问题那就是各个样本的表达量，看自己喜欢用就好：

p=ggplot(exp_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)

p=ggplot(exp_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
print(p)

p=ggplot(exp_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print(p)

p=ggplot(exp_L,aes(value,col=group))+geom_density() 
print(p)

3.3 检查样本分组信息

检查样本分组信息，一般看PCA图，hclust图

hclust

# 更改表达矩阵列名
head(exp)
colnames(exp) = paste(group_list,1:6,sep='')
head(exp)
# 定义nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19), 
                cex = 0.7, col = "blue")
# 聚类
hc=hclust(dist(t(exp)))
par(mar=c(5,5,5,10)) 
# 绘图
plot(as.dendrogram(hc), nodePar = nodePar,  horiz = TRUE)

画出图后我们发现，control和treatment很好的分开了，组内也很好的聚类到了一起说明数据过关。

PCA

library(ggfortify)
# 互换行和列，再dim一下
df=as.data.frame(t(exp))
# 不要view df，列太多，软件会卡住；
dim(df)
dim(exp)

exp[1:6,1:6]
df[1:6,1:6]

df$group=group_list 
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')
save(exp,group_list,file = "step2output.Rdata")

同样发现该分开的分开了，该聚在一起的聚在一起了，数据很好，符合预期。


PCA图

经过上面一系列的表达矩阵可视化，我们检查了表达矩阵发现是正确的，接下来就要做差异分析啦

4. 差异分析及可视化

芯片数据做差异分析最常用的就是limma包
使用这个包需要三个数据：

  • 表达矩阵(exp)
  • 分组矩阵(design)
  • 差异比较矩阵(contrast.matrix)

下面我们开始准备这三个输入数据：
表达矩阵(exp)我们早就得到了，不用再制作了；
我们也得到了存放分组信息的向量group_list，用它来制作我们的分组矩阵

4.1 limma包做差异分析输入数据的准备

输入数据—分组矩阵

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = "step2output.Rdata")
dim(exp)
library(limma)
# 做分组矩阵 
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exp)
design  #得到的分组矩阵
分组矩阵design：1代表是；0代表不是

输入数据—差异比较矩阵

> contrast.matrix<-makeContrasts(paste0(c("treat","contorl"),collapse = "-"),levels = design)
> contrast.matrix
         Contrasts
Levels    treat-contorl
  contorl            -1
  treat               1

contrast.matrix 这个矩阵声明，我们要把treat组和contor组进行差异分析比较：-1和1的意思是contorl是用来被比的，treat是来比的即：treat/contorl

到此，做差异分析所需要的三个矩阵就做好了：表达矩阵(exp)、分组矩阵(design)、差异比较矩阵(contrast.matrix)
我们已经制作好了必要的输入数据，下面开始讲如何使用limma包来进行差异分析！

4.2 limma包做差异分析

只有三个步骤：

  • lmFit
  • eBayes
  • topTable
##step1
fit <- lmFit(exp,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2)  ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
save(nrDEG,file = "DEGoutput.Rdata")

此时我们就得到差异分析矩阵（nrDEG），重点看logFC和P值：

差异分析结果

差异分析就是对每个基因都进行检验，检验基因的logFG是多大、平均表达量是多少、p.value是否显著等...

4.3 差异表达基因的可视化

用limma包得到差异分析表达矩阵后作图检查差异基因是否真的很差异

画热图

选差异最显著的前25个基因画热图，查看差异是否真的很显著

##热图
rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = "DEGoutput.Rdata")
load(file = "DEGinput.Rdata")

library(pheatmap)
choose_gene = head(rownames(nrDEG),25)
choose_matrix = exp[choose_gene,]
choose_matrix = t(scale(t(choose_matrix)))
pheatmap(choose_matrix)

火山图

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = "DEGoutput.Rdata")
colnames(nrDEG)
plot(nrDEG$logFC,-log10(nrDEG$P.Value))

DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                              ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
)
this_tile
head(DEG)
g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle( this_tile  ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red'))  ## corresponding to the levels(res$change)
print(g)

5.富集分析—KEGG、 GO

富集分析就是用常用的数据库来注释基因list。差异分析通过自定义的阈值挑选了有统计学显著的基因列表后我们其实是需要对它们进行注释才能了解其功能，最常见的就是GO/KEGG数据库注释，当然也可以使用Reactome和Msigdb数据库来进行注释。而最常见的注释方法就是超几何分布检验。超几何分布检验，运用到通路的富集概念就是“总共有多少基因（这个地方值得注意，主流认为只考虑那些在KEGG等数据库注释的背景基因），你的通路有多少基因，你的通路被抽中了多少基因（在差异基因里面属于你的通路的基因）。” 目的就是知道，哪些通路中的哪些基因的表达因为药物或者某些操作的作用发生了较大的变化，导致通路有较大改变。

5.1 KEGG pathway analysis

clusterProfiler包作KEGG富集分析

#clusterProfiler作kegg富集分析：
  library(clusterProfiler)
  gene_up= deg[deg$change == 'up','ENTREZID'] 
  gene_down=deg[deg$change == 'down','ENTREZID'] 
  gene_diff=c(gene_up,gene_down)
  gene_all = deg[,'ENTREZID']
  kk.up <- enrichKEGG(gene         = gene_up,
                      organism     = 'hsa',
                      universe     = gene_all,
                      pvalueCutoff = 0.9,
                      qvalueCutoff =0.9)
  head(kk.up)[,1:6]
  dim(kk.up)
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa',
                        universe     = gene_all,
                        pvalueCutoff = 0.9,
                        qvalueCutoff =0.9)
  head(kk.down)[,1:6]
  dim(kk.down)
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa',
                        pvalueCutoff = 0.05)
  head(kk.diff)[,1:6]
  
  class(kk.diff)
  #提取出数据框
  kegg_diff_dt <- kk.diff@result
  
  #根据pvalue来选,用于可视化
  down_kegg <- kk.down@result %>%
    filter(pvalue<0.05) %>%
    mutate(group=-1)
  
  up_kegg <- kk.up@result %>%
    filter(pvalue<0.05) %>%
    mutate(group=1)
  
  #可视化走起
  kegg_plot <- function(up_kegg,down_kegg){
    dat=rbind(up_kegg,down_kegg)
    colnames(dat)
    dat$pvalue = -log10(dat$pvalue)
    dat$pvalue=dat$pvalue*dat$group 
    
    dat=dat[order(dat$pvalue,decreasing = F),]
    
    g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + 
      geom_bar(stat="identity") + 
      scale_fill_gradient(low="blue",high="red",guide = FALSE) + 
      scale_x_discrete(name ="Pathway names") +
      scale_y_continuous(name ="log10P-value") +
      coord_flip() + theme_bw()+theme(plot.title = element_text(hjust = 0.5))+
      ggtitle("Pathway Enrichment") 
  }
  
  g_kegg <- kegg_plot(up_kegg,down_kegg)
  g_kegg
  
  ggsave(g_kegg,filename = 'kegg_up_down.png')
clusterProfiler作KEGG富集分析

GSEA作KEGG富集分析

GSEA是另一个常用的富集分析方法，目的是看看基因全局表达量的变化是否有某些特定的基因集合的倾向性。

  data(geneList, package="DOSE")
  head(geneList)
  length(geneList)
  names(geneList)
  boxplot(geneList)
  boxplot(deg$logFC)
  
  geneList=deg$logFC
  names(geneList)=deg$ENTREZID
  geneList=sort(geneList,decreasing = T)
  
  kk_gse <- gseKEGG(geneList     = geneList,
                    organism     = 'hsa',
                    nPerm        = 1000,
                    minGSSize    = 120,
                    pvalueCutoff = 0.9,
                    verbose      = FALSE)
  head(kk_gse)[,1:6]
  gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))
  
  down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
  up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
  
  gse_kegg=kegg_plot(up_kegg,down_kegg)
  print(gse_kegg)
  ggsave(gse_kegg,filename ='kegg_up_down_gsea.png')
GSEA作KEGG富集分析

5.2 GO database analysis

做GO数据集超几何分布检验分析，重点在结果的可视化及生物学意义的理解。

GO富集分析生物学意义

GO富集分析原理：有一个term注释了100个差异表达基因参与了哪个过程，注释完之后（模式生物都有现成的注释包，不用我们自己注释），计算相对于背景它是否显著集中在某条通路、某一个细胞学定位、某一种生物学功能。
对GO database analysis一般从三个层面进行：

  • Cellular component，CC 细胞成分
  • Biological process， BP 生物学过程
  • Molecular function，MF 分子功能

这三个层面具体是指：

  • Cellular component解释的是基因存在在哪里，在细胞质还是在细胞核？如果存在细胞质那在哪个细胞器上？如果是在线粒体中那是存在线粒体膜上还是在线粒体的基质当中？这些信息都叫Cellular component。
  • Biological process是在说明该基因参与了哪些生物学过程，比如，它参与了rRNA的加工或参与了DNA的复制，这些信息都叫Biological process
  • Molecular function在讲该基因在分子层面的功能是什么？它是催化什么反应的？
    立足于这三个方面，我们将得到基因的注释信息。

GO富集分析的R代码

#go富集分析--耗费时间灰常长，很正常
library(clusterProfiler)
#输入数据
gene_up= deg[deg$change == 'up','ENTREZID'] 
gene_down=deg[deg$change == 'down','ENTREZID'] 
gene_diff=c(gene_up,gene_down)
head(deg)

  #**GO分析三大块**
#细胞组分
    ego_CC <- enrichGO(gene = gene_diff,
                       OrgDb= org.Hs.eg.db,
                       ont = "CC",
                       pAdjustMethod = "BH",
                       minGSSize = 1,
                       pvalueCutoff = 0.01,
                       qvalueCutoff = 0.01,
                       readable = TRUE)
    #生物过程
    ego_BP <- enrichGO(gene = gene_diff,
                       OrgDb= org.Hs.eg.db,
                       ont = "BP",
                       pAdjustMethod = "BH",
                       minGSSize = 1,
                       pvalueCutoff = 0.01,
                       qvalueCutoff = 0.01,
                       readable = TRUE)
    #分子功能：
    ego_MF <- enrichGO(gene = gene_diff,
                       OrgDb= org.Hs.eg.db,
                       ont = "MF",
                       pAdjustMethod = "BH",
                       minGSSize = 1,
                       pvalueCutoff = 0.01,
                       qvalueCutoff = 0.01,
                       readable = TRUE)
    save(ego_CC,ego_BP,ego_MF,file = "ego_GPL6244.Rdata")
    rm(list = ls())  ## 魔幻操作，一键清空~   
  load(file = "ego_GPL6244.Rdata")
  #作图
  #第一种，条带图，按p从小到大排的
  barplot(ego_CC, showCategory=20,title="EnrichmentGO_CC")
  barplot(ego_BP, showCategory=20,title="EnrichmentGO_CC")
  #如果运行了没出图，就dev.new()
  #第二种，点图，按富集数从大到小的
  dotplot(ego_CC,title="EnrichmentGO_BP_dot")
  
  #保存
  pdf(file = "dotplot_GPL6244.pdf")
  dotplot(ego_CC,title="EnrichmentGO_BP_dot")
  dev.off()

纯代码版：


#GEO B站视频 纯代码篇

#下载加载包
cran_packages <- c('tidyr',
                   'tibble',
                   'dplyr',
                   'stringr',
                   'ggplot2',
                   'ggpubr',
                   'factoextra',
                   'FactoMineR',
                   'WGCNA') 
Biocductor_packages <- c('GEOquery',
                         'hgu133plus2.db',
                         "KEGG.db",
                         "limma",
                         "impute",
                         "GSEABase",
                         "GSVA",
                         "clusterProfiler",
                         "genefu",
                         "org.Hs.eg.db",
                         "preprocessCore",
                         "hugene10sttranscriptcluster.db")

for (pkg in c(Biocductor_packages,cran_packages)){
  require(pkg,character.only=T) 
}

# 下载数据
rm(list = ls())
library(GEOquery)
eSet <- getGEO("GSE42872", 
               destdir = '.',
               getGPL = F)


# 从eSet中提取表达矩阵exp
exp <- exprs(eSet[[1]])

head(exp)

# ID转换

##探针ID（probe_id）转换成symbol ID

eSet[[1]]@annotation
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")
ids=toTable(hugene10sttranscriptclusterSYMBOL)
head(ids)
length(unique(ids$symbol))
tail(sort(table(ids$symbol)))

table(sort(table(ids$symbol)))

table(rownames(exp) %in% ids$probe_id)
dim(exp)
exp = exp[rownames(exp) %in% ids$probe_id,]
dim(exp)


ids=ids[match(rownames(exp),ids$probe_id),]
head(ids)
head(exp)

tmp = by(exp,
         ids$symbol,
         function(x) rownames(x)[which.max(rowMeans(x))])

probes = as.character(tmp)
head(tmp)
head(probes)

dim(exp)
exp = exp[rownames(exp) %in% probes,]
dim(exp)
rownames(exp)=ids[match(rownames(exp),ids$probe_id),2]
head(exp)

pd <- pData(eSet[[1]]) # pData函数得到每个样本的描述信息
head(pd)

save(pd,exp,file = "step1output.Rdata")
save(exp,file = "DEGinput.Rdata")

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = "step1output.Rdata")

#####

#检查表达矩阵
##画典型基因表达量的boxplot
exp['GAPDH',] 
exp['ACTB',]
boxplot(exp)
boxplot(exp['GAPDH',])
boxplot(exp['ACTB',])

#各个样本表达量的boxplot
##准备画图所需数据exp_L
library(reshape2)
head(exp)
exp_L = melt(exp)
head(exp_L)
colnames(exp_L)=c('symbol','sample','value')
head(exp_L)
# 获得分组信息
library(stringr)

group_list = ifelse(str_detect(pd$title,"Control")==TRUE,"contorl","treat")

group_list

exp_L$group = rep(group_list,each=nrow(exp))
head(exp_L)
table(exp_L[,2])
dim(exp_L)
### ggplot2画图 
library(ggplot2)
p = ggplot(exp_L,
         aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)

# 对表达矩阵进行聚类绘图，并添加样本的临床表型数据信息(更改样本名)

## hclust
# 更改表达矩阵列名
head(exp)
colnames(exp) = paste(group_list,1:6,sep='')
head(exp)
# 定义nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19), 
                cex = 0.7, col = "blue")
# 聚类
hc=hclust(dist(t(exp)))
par(mar=c(5,5,5,10)) 
# 绘图
plot(as.dendrogram(hc), nodePar = nodePar,  horiz = TRUE)

## PCA

library(ggfortify)
# 互换行和列，dim一下
head(exp)
df=as.data.frame(t(exp))
# 不要view df，列太多，软件会崩掉；
dim(df)
dim(exp)

exp[1:6,1:6]
df[1:6,1:6]

df$group=group_list 
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')

save(exp,group_list,file = "step2output.Rdata")


###################################################
############用limma对芯片数据做差异分析############
###################################################

#差异分析——limma
rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = "step2output.Rdata")
dim(exp)

library(limma)
# 做分组矩阵 
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exp)
design  #分组矩阵

# 做比较矩阵

# contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
# contrast.matrix ##这个矩阵声明，我们要把treat组和contorl组进行差异分析比较
# -1和1的意思是contorl是用来被比的，treat是来比的
contrast.matrix<-makeContrasts(paste0(c("treat","contorl"),collapse = "-"),levels = design)
contrast.matrix
#到此，做差异分析所需要的三个矩阵就做好了：表达矩阵、分组矩阵、差异比较矩阵
#我们已经制作好了必要的输入数据，下面开始讲如何使用limma这个包来进行差异分析了！

##step1
fit <- lmFit(exp,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2)  ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)

save(exp,group_list,nrDEG,file = "DEGoutput.Rdata")

#用limma包得到差异分析表达矩阵后作图检查差异基因是否真的很差异
##热图
rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = "DEGoutput.Rdata")
load(file = "DEGinput.Rdata")

library(pheatmap)
choose_gene = head(rownames(nrDEG),25)
choose_matrix = exp[choose_gene,]
choose_matrix = t(scale(t(choose_matrix)))
pheatmap(choose_matrix)


##火山图
rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = "DEGoutput.Rdata")
colnames(nrDEG)
plot(nrDEG$logFC,-log10(nrDEG$P.Value))

DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                              ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
)
this_tile
head(DEG)
g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle( this_tile  ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red'))  ## corresponding to the levels(res$change)
print(g)

#富集分析

#富集分析准备工作：

##首先对差异表达矩阵nrDEG，进行加工
###1.把行名变成SYMBOL列
rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = "DEGoutput.Rdata")
library(dplyr)
deg = nrDEG
deg <- mutate(deg,symbol = rownames(deg))
head(deg)

###2.加change列：上调或下调，火山图要用

logFC_t = 1 #不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
change=ifelse(deg$P.Value>0.01,'stable', 
              ifelse( deg$logFC >logFC_t,'up', 
                      ifelse( deg$logFC < -logFC_t,'down','stable') )
)
deg <- mutate(deg,change)
head(deg)
table(deg$change)

###3.加ENTREZID列，后面富集分析要用
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(unique(deg$symbol), fromType = "SYMBOL",  #ID转换核心函数bitr
            toType = c( "ENTREZID"),
            OrgDb = org.Hs.eg.db)
head(s2e)
head(deg)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))

head(deg)

save(exp,group_list,deg,file = "enrich_input.Rdata")

#####################
######富集分析#######
#####################

rm(list = ls()) 
options(stringsAsFactors = F)
load(file = 'enrich_input.Rdata')

## 1.KEGG pathway analysis
#上调、下调、差异、所有基因

#clusterProfiler作kegg富集分析：
  library(clusterProfiler)
  gene_up= deg[deg$change == 'up','ENTREZID'] 
  gene_down=deg[deg$change == 'down','ENTREZID'] 
  gene_diff=c(gene_up,gene_down)
  gene_all = deg[,'ENTREZID']
  kk.up <- enrichKEGG(gene         = gene_up,
                      organism     = 'hsa',
                      universe     = gene_all,
                      pvalueCutoff = 0.9,
                      qvalueCutoff =0.9)
  head(kk.up)[,1:6]
  dim(kk.up)
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa',
                        universe     = gene_all,
                        pvalueCutoff = 0.9,
                        qvalueCutoff =0.9)
  head(kk.down)[,1:6]
  dim(kk.down)
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa',
                        pvalueCutoff = 0.05)
  head(kk.diff)[,1:6]
  
  class(kk.diff)
  #提取出数据框
  kegg_diff_dt <- kk.diff@result
  
  #根据pvalue来选,用于可视化
  down_kegg <- kk.down@result %>%
    filter(pvalue<0.05) %>%
    mutate(group=-1)
  
  up_kegg <- kk.up@result %>%
    filter(pvalue<0.05) %>%
    mutate(group=1)
  
  #可视化
  kegg_plot <- function(up_kegg,down_kegg){
    dat=rbind(up_kegg,down_kegg)
    colnames(dat)
    dat$pvalue = -log10(dat$pvalue)
    dat$pvalue=dat$pvalue*dat$group 
    
    dat=dat[order(dat$pvalue,decreasing = F),]
    
    g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + 
      geom_bar(stat="identity") + 
      scale_fill_gradient(low="blue",high="red",guide = FALSE) + 
      scale_x_discrete(name ="Pathway names") +
      scale_y_continuous(name ="log10P-value") +
      coord_flip() + theme_bw()+theme(plot.title = element_text(hjust = 0.5))+
      ggtitle("Pathway Enrichment") 
  }
  
  g_kegg <- kegg_plot(up_kegg,down_kegg)
  g_kegg
  
  ggsave(g_kegg,filename = 'kegg_up_down.png')


#gsea作kegg富集分析：

  data(geneList, package="DOSE")
  head(geneList)
  length(geneList)
  names(geneList)
  boxplot(geneList)
  boxplot(deg$logFC)
  
  geneList=deg$logFC
  names(geneList)=deg$ENTREZID
  geneList=sort(geneList,decreasing = T)
  
  kk_gse <- gseKEGG(geneList     = geneList,
                    organism     = 'hsa',
                    nPerm        = 1000,
                    minGSSize    = 120,
                    pvalueCutoff = 0.9,
                    verbose      = FALSE)
  head(kk_gse)[,1:6]
  gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))
  
  down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
  up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
  
  gse_kegg=kegg_plot(up_kegg,down_kegg)
  print(gse_kegg)
  ggsave(gse_kegg,filename ='kegg_up_down_gsea.png')


### 2.GO database analysis 

#go富集分析
library(clusterProfiler)
#输入数据
gene_up= deg[deg$change == 'up','ENTREZID'] 
gene_down=deg[deg$change == 'down','ENTREZID'] 
gene_diff=c(gene_up,gene_down)
head(deg)

  #**GO分析三大块**
#细胞组分
    ego_CC <- enrichGO(gene = gene_diff,
                       OrgDb= org.Hs.eg.db,
                       ont = "CC",
                       pAdjustMethod = "BH",
                       minGSSize = 1,
                       pvalueCutoff = 0.01,
                       qvalueCutoff = 0.01,
                       readable = TRUE)
    #生物过程
    ego_BP <- enrichGO(gene = gene_diff,
                       OrgDb= org.Hs.eg.db,
                       ont = "BP",
                       pAdjustMethod = "BH",
                       minGSSize = 1,
                       pvalueCutoff = 0.01,
                       qvalueCutoff = 0.01,
                       readable = TRUE)
    #分子功能：
    ego_MF <- enrichGO(gene = gene_diff,
                       OrgDb= org.Hs.eg.db,
                       ont = "MF",
                       pAdjustMethod = "BH",
                       minGSSize = 1,
                       pvalueCutoff = 0.01,
                       qvalueCutoff = 0.01,
                       readable = TRUE)
    save(ego_CC,ego_BP,ego_MF,file = "ego_GPL6244.Rdata")
    rm(list = ls()) 
  load(file = "ego_GPL6244.Rdata")

  #第一种，条带图，按p从小到大排的
  barplot(ego_CC, showCategory=20,title="EnrichmentGO_CC")
  barplot(ego_BP, showCategory=20,title="EnrichmentGO_CC")
  #如果运行了没出图，就dev.new()
  #第二种，点图，按富集数从大到小的
  dotplot(ego_CC,title="EnrichmentGO_BP_dot")
  
  #保存
  pdf(file = "dotplot_GPL6244.pdf")
  dotplot(ego_CC,title="EnrichmentGO_BP_dot")
  dev.off()

特别感谢小洁老师激发了我学习GEO数据库挖掘的兴趣；有些图片还有富集分析的代码就来自小洁老师的课件哦

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 206,311评论 6赞 481
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 88,339评论 2赞 382
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 152,671评论 0赞 342
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 55,252评论 1赞 279
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 64,253评论 5赞 371
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,031评论 1赞 285
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,340评论 3赞 399
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,973评论 0赞 259
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 43,466评论 1赞 300
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,937评论 2赞 323
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,039评论 1赞 333
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,701评论 4赞 323
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,254评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,259评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,485评论 1赞 262
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 45,497评论 2赞 354
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,786评论 2赞 345