在使用Bioconductor 分析芯片数据（一）中我们使用的是芯片的原始数据（即*.cel文件）进行分析，但在大部分情况下我们是可以利用GEO数据库中已有的表达矩阵来进行数据分析的。
本次我们选用GSE17975数据，这是一个双通道芯片数据，研究Kawasaki disease的致病源头。但是奇怪的是作者没有上传芯片源文件（*.cel文件），因此我们直接使用 表达矩阵数据进行分析。

除了芯片外，我们还选择了一个Kawasaki disease研究相关转录组数据，利用这两种方法得到的差异基因的交集进行后续分析。

1、下载数据

这里依旧使用GEOquery包来下载数据：

> setwd('D:/TCGA/microarray_analysis/GSE17975/data/')
# 直接下载表达矩阵
> gmat <- getGEO('GSE17975',destdir = '.')
# 查看下里面是什么
> gmat
$GSE17975_series_matrix.txt.gz
ExpressionSet (storageMode: lockedEnvironment)
assayData: 22839 features, 3 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: GSM450153 GSM450154 GSM450155
  varLabels: title geo_accession ... disease state:ch2 (48 total)
  varMetadata: labelDescription
featureData
  featureNames: AGhsA010102 AGhsA010103 ... AGhsS010824 (22839 total)
  fvarLabels: ID Description ... SPOT_ID (13 total)
  fvarMetadata: Column Description labelDescription
experimentData: use 'experimentData(object)'
Annotation: GPL1291

> class(gmat)
[1] "list"
> class(gmat[[1]])
[1] "ExpressionSet"
attr(,"package")
[1] "Biobase"
# 得到表达矩阵
> exprset <- exprs(gmat[[1]]); dim(exprset)
[1] 22839     3  # 可以看到表达矩阵中共有22839个基因（features），3个样本
# 使用pData还可以查看更详细的实验相关信息
> pdata <- pData(exprset)； dim(pdata)
[1]  3 48   # 有3行48列实验数据信息

### 最后查看一下表达矩阵信息
> head(exprset)
            GSM450153 GSM450154 GSM450155
AGhsA010102    0.8335    0.0257    0.0841
AGhsA010103    0.0649   -0.3715   -0.2740
AGhsA010105   -0.3237   -0.7539   -0.2793
AGhsA010106    0.4844   -0.6897   -0.4521
AGhsA010108   -1.7418   -0.8783   -3.5804
AGhsA010109   -0.3401   -1.2311   -0.5332
> class(exprset)
[1] "matrix"
# 转换成data.frame数据结构方便后面使用
> exprset <- as.data.frame(exprset)

从上面可以发现，第一，虽然我们是直接下载的表达矩阵，但是其中也含有许多的实验数据信息以及芯片平台信息；第二通过直接查看表达矩阵他已经经过了标准化。

GEO Platform (GPL) 芯片平台
GEO Sample (GSM) 样本ID号
GEO Series (GSE) study的ID号
GEO Dataset (GDS) 数据集的ID号

2、筛选差异基因

limma进行差异表达分析需要准备好3个数据：表达矩阵，分组矩阵，差异比较矩阵。表达矩阵我们上面已经从GEO下载得到，只需要将芯片ID注释为基因ID即可（也可在差异分析后进行注释ID），但这里的芯片平台比较罕见，找了半个小时没找到对应的注释包，因此直接下载其对应的注释包写代码进行注释。

> gpl1291 <- getGEO(filename='GPL1291.soft')
> colnames(Table(gpl1291))
 [1] "ID"                    "Description"           "CDS_ID"               
 [4] "GB_ACC"                "Entrez Gene ID"        "Symbol"               
 [7] "GeneName"              "IPI ID"                "GO Molecular Function"
[10] "GO Biological Process" "GO Cellular Component" "ORF"                  
[13] "SPOT_ID"
# 从上面可以看到GPL1291.soft文件中第一列为芯片ID，第5列和第6列为我们常用的ID，这里使用第6列Gene Symbol进行注释
> symbol <- Table(gpl1291)[c('ID','Symbol')]; head(symbol)
           ID  Symbol
1 AGhsA010101   IFNA8
2 AGhsA010102    OGG1
3 AGhsA010103 KIR2DL2
4 AGhsA010104   FGFR2
5 AGhsA010105   GAGE1
6 AGhsA010106     VCX
### 最后我们只需要将两个数据框merge一下就可以了
 > exprset1 <- merge(x=exprset,y=symbol,by='ID',all.x=T);head(exprset1)
           ID GSM450153 GSM450154 GSM450155  Symbol
1 AGhsA010102    0.8335    0.0257    0.0841    OGG1
2 AGhsA010103    0.0649   -0.3715   -0.2740 KIR2DL2
3 AGhsA010105   -0.3237   -0.7539   -0.2793   GAGE1
4 AGhsA010106    0.4844   -0.6897   -0.4521     VCX
5 AGhsA010108   -1.7418   -0.8783   -3.5804  EEF1A1
6 AGhsA010109   -0.3401   -1.2311   -0.5332    DAZ1
### 查看一下注释有没有一对多或者未注释上的情况
> dim(exprset);dim(exprset1)
[1] 22839     4
[1] 22839     5
> exprset1.na <- exprset1[is.na(exprset1$Symbol)]; dim(exprset1.na)
[1] 22839     0
### 最后删除ID列，并把symbol列调到第一列
> exprset2 <- exprset1[,-1]
> exprset3 <- exprset2[c(4,1,2,3)]; head(exprset3)
   Symbol GSM450153 GSM450154 GSM450155
1    OGG1    0.8335    0.0257    0.0841
2 KIR2DL2    0.0649   -0.3715   -0.2740
3   GAGE1   -0.3237   -0.7539   -0.2793
4     VCX    0.4844   -0.6897   -0.4521
5  EEF1A1   -1.7418   -0.8783   -3.5804
6    DAZ1   -0.3401   -1.2311   -0.5332

从上面可以看到芯片ID注释的结果很好，全部注释上了而且没有一对多的情况。
但是真的没有问题吗？会不会不是NA而是别的错误没有被我们发现？让我们再来做一次检查。

### 查看一下有没有空白字符
> failed <- exprset3[exprset3$Symbol=='',]
> head(failed)
    Symbol GSM450153 GSM450154 GSM450155
16           -0.6712   -0.9269   -1.8417
50           -1.2345    0.1686   -0.5564
56           -0.5735   -1.0439   -1.5821
100           2.1925   -0.2059    0.4510
102          -2.1779   -0.6038   -0.3111
120          -0.4600    0.3663   -0.5886
> dim(failed)
[1] 7025    4   # (ﾟДﾟ≡ﾟДﾟ) 惊呆了有没有，竟然有这么多没有注释上
### 难道是merge出错了？不应该呀，查看一下官方给的注释文件
> official.failed <- symbol[symbol$Symbol=='',]
> dim(official.failed)
[1] 10985     2     # 再次惊呆，官方给的注释文件竟然有这么多空白。。。
### 最后，去除这些未注释的ID
> exprset4 <- exprset3[exprset3$Symbol!='',]; dim(exprset4)
[1] 15814     4

最后我们得到了15814个注释到的基因，相比之前少了7025个空的。所以，看起来简简单单的一步注释，也到处都是坑啊。
由于数据比较老，这里各列的数据其实都是表示logFC，原文作者也是直接根据logFC选择的差异基因进行后续分析的，因此这里我们也是直接使用3组平行试验均值的|logFC|>1筛选差异基因进行后续分析。

> exprset4$mean <- (exprset4$GSM450153+exprset4$GSM450154+exprset4$GSM450155)/3
# 设定阈值|mean| >1 来筛选差异基因
> microexprgenes.differ <- exprset4[abs(exprset4$mean)>1,]
> microexprgenes.differ <- microexprgenes.differ[c(1,2,3,4)]
> dim(microexprgenes.differ)
[1] 938   4

这里我们筛选得了938个差异基因。然后去一下重复的gene symbol：

> microexprgenes.differ.dup <- microexprgenes.differ[duplicated(microexprgenes.differ$Symbol),]
> dim(microexprgenes.differ.dup)
[1] 99  4
> head(microexprgenes.differ.dup)
       Symbol GSM450153 GSM450154 GSM450155
28     EEF1A1   -1.5649   -1.5187   -4.8573
677    FCGR1A    2.7853    2.5568    2.1375
1304    MATR3   -0.6394   -1.8890   -1.5522
2662 CCNB1IP1   -1.3921   -1.7959   -1.0923
3349     DHX9   -0.6461   -1.2042   -3.1203
4523    PSMB3    1.4895    1.1023    1.1603
> microexprgenes.differ <- microexprgenes.differ[!duplicated(microexprgenes.differ$Symbol),]
> dim(microexprgenes.differ)
[1] 839   4

最后我们得到了839个差异基因。

3、转录组差异基因筛选

这里直接使用的数据同样也是关于Kawasaki disease研究的，GEO数据集为GSE64486，我们这里直接采用已经处理好的表达矩阵，但是我们这里使用的表达矩阵数据为GSE64486_untreated_vs_control_foldgenes_broad.txt.gz，即未治疗的与对照组。
由于这里的数据是处理好的，我们可以直接设置logFC和P-value阈值来筛选差异基因。这里设置的阈值为P-value<0.05, |logFC|>1。

# 先读取数据
> exprgenes <- read.table('GSE64486_untreated_vs_control_foldgenes_broad.txt',header = T, stringsAsFactors = F,sep='\t');head(exprgenes)
               ID      P.Value      Q.Value    adj.P.Val Gene.Symbol RefSeq.ID Entrez.ID Fold.Change
1 ENSG00000211896 7.601798e-25 1.645219e-20 1.645219e-20       IGHG1        NA      3500   101.77408
2 ENSG00000138755 3.384964e-09 1.495083e-06 1.495083e-06       CXCL9        NA      4283    45.22891
3 ENSG00000211897 3.600406e-16 8.657976e-13 8.657976e-13       IGHG3        NA      3502    45.06493
4 ENSG00000211899 9.557181e-17 2.711637e-13 2.711637e-13        IGHM        NA      3507    41.89091
5 ENSG00000132465 1.444665e-14 2.233298e-11 2.233298e-11         IGJ        NA      3512    36.47332
6 ENSG00000196735 3.890287e-21 2.806518e-17 2.806518e-17    HLA-DQA1        NA      3117    35.58685
# 筛选差异基因
> exprgenes.differ <- exprgenes[exprgenes$P.Value<0.05 & abs(exprgenes$Fold.Change)>1,]
> dim(exprgenes);dim(exprgenes.differ)
[1] 43285     8
[1] 7099    8

可以看到，这里一共筛选得到了7099个，得到的基因比较多，可能阈值太宽了，暂且先这样接着做吧。
最后，还是检查一下作者提供的数据注释的gene symbol有没有空的或者重复的：

> failed.symbol <- exprgenes.differ[duplicated(exprgenes.differ$Gene.Symbol),]
> head(failed.symbol)
                 ID      P.Value      Q.Value    adj.P.Val Gene.Symbol RefSeq.ID Entrez.ID Fold.Change
120 ENSG00000252493 1.048966e-09 5.470419e-07 5.470419e-07                    NA        NA    7.509095
124 ENSG00000235300 3.582412e-09 1.566310e-06 1.566310e-06                    NA        NA    7.389821
132 ENSG00000251320 4.767764e-13 5.291607e-10 5.291607e-10                    NA        NA    7.135815
160 ENSG00000252862 7.716183e-16 1.757868e-12 1.757868e-12                    NA        NA    6.585053
167 ENSG00000251301 1.863953e-05 1.750134e-03 1.750134e-03                    NA        NA    6.419032
202 ENSG00000262151 4.473742e-05 3.500927e-03 3.500927e-03
> dim(failed.symbol)
[1] 2281    8
### 去除那些没有gene symbol 信息的
> exprgenes.differ <- exprgenes.differ[!duplicated(exprgenes.differ$Gene.Symbol),]
> dim(exprgenes.differ)
[1] 4818    8

可以看到，一共2281个没有注释到gene symbol，去除后一共得到了4818个有gene symbol信息的差异基因。

我们可以发现，无论是芯片数据还是转录组数据，ID注释的空白或者重复都是一个容易忽略的问题。有机会再进行较为深入的研究一下。

4、共有差异基因选取

前面我们分别得到了关于川崎病研究的芯片和转录组的差异基因，现在我们将这两种数据合并取交集，以期得到与川崎病相关更为密切的基因。这里操作也比较简单，直接使用R的merge函数即可。

> intersectgenes <- merge(x=microexprgenes.differ, y=exprgenes.differ, all=F,by.x='Symbol',by.y = 'Gene.Symbol')
> dim(intersectgenes)
[1] 128  11
> head(intersectgenes)
    Symbol GSM450153 GSM450154 GSM450155              ID      P.Value     Q.Value   adj.P.Val RefSeq.ID Entrez.ID Fold.Change
1     AIM1   -0.7155   -0.8391   -2.3808 ENSG00000112297 4.518321e-02 0.287315306 0.287315306        NA       202    2.555354
2      AK5   -1.9269   -0.6967   -0.5628 ENSG00000154027 7.907520e-05 0.005415775 0.005415775        NA     26289    1.350051
3    ANXA3    4.2110    0.8687   -0.1016 ENSG00000138772 3.122002e-02 0.229588586 0.229588586        NA       306   -2.328736
4 ARHGAP15   -0.0725   -1.5821   -2.0469 ENSG00000075884 6.321948e-05 0.004553170 0.004553170        NA     55843    5.064183
5    ASGR2    1.5563    1.4054    1.2066 ENSG00000161944 1.474010e-02 0.146976591 0.146976591        NA       433    2.061535
6      ATM   -1.4344   -0.4961   -1.9324 ENSG00000149311 1.403235e-02 0.142613364 0.142613364        NA       472    2.313447

可以看到，这里我们一共得到了128个共有基因，我们将使用这128个共有基因进行后续分析。讲这些基因写入文件保存方便后面使用。

> write.table(intersectgenes,'../intersectgenes_logFC_broad.txt',row.names = F,sep='\t')
> write.table(intersectgenes$Symbol,'../intersectgenes.txt',row.names = F,sep='\t')

做个韦恩图可视化看看：

vn <- venn.diagram(list(tissu=exprgenes.differ$Gene.Symbol, peripheral.blood=microexprgenes.differ$Symbol),filename='../intersectgenes.png',fill=c('red','green'),alpha=c(0.5,0.5),lty=2,cat.fontface=2,col='black',cat.col=c('red','green'))

真丑！！！！

总结

在这一部分，我们利用了两份GEO数据库中的内容，其中一份是芯片数据，另一个是转录组数据，对这两个实验数据线选定阈值筛选出差异基因，然后综合两个数据寻找共同的差异基因。
这部分操作比较简单，都是直接用的原作者给出的表达矩阵，然后进行merge取交集即可。但是在操作过程中仍然有两点是值得注意的：第一是各种基因ID之间的相互转换关系，这里可以利用各个数据库提供的ID对应文件进行转换即可；第二点，在本次实践中我们发现了大量的未成功注释的、空白的或重复的gene symbol，这可以利用R中的duplicated函数在每次进行数据分析时进行一次检查。