基因芯片的简单分析

本文主要介绍一下基因芯片的简单分析过程。

文档参考

Using Bioconductor for Microarray Analysis

affy

limma # 能够常常看看这个，比较全面

Processing Affymetrix Gene Expression Arrays

Tutorial: Analysing Microarray Data In Bioconductor

R包安装

主要须要两个包， affy, limma

> source("http://bioconductor.org/biocLite.R")
> biocLite("affy")
> biocLite("limma")

案例数据

GEO accession: GSE29349
GSE29349
是以水稻为研究对象， 分为野生型与突变体， 每组有三个重复。下载原始数据 GSE29349_RAW.tar 便可，解压，文件夹GSE29349_RAW里面包含了原始 .CEL 文件。

数据读入

> library(affy)
> myData <- ReadAffy(celfile.path="/media/文档/microarry/GSE29349_RAW/")

如今 myData 是一个 AffyBatch object , CEL 格式的信息都被读入到 myData 中。

数据查看

> myData
> str(myData)
> pData(myData)
> phenoData(myData)
> annotation(myData)
> probeNames(myData)

RMA normalization, 以及数据提取

> eset=rma(myData)
# 利用 rma 进行 Background correcting, Normalizing, Calculating Expression
> class(eset)
# 能够看到目前 eset 是一个 ExpressionSet object
# 后面使用 limma 包进行分析时就会用到 eset

固然，也能够同时将通过 rma 处理过的数据提取处理，保存在一个文本中。

> rma <- exprs(eset)
> class(rma)
# 目前 rma 是一个 matrix
> rma
> rma=format(rma, digits=6)
# 设置小数点位数
> write.table(rma, file="rma.txt", quote=FALSE, sep="\t")
# 数据写入
# 打开该文本，格式以下:
GSM725469_61_HB_Rice_.CEL.gz    GSM725470_62_HB_Rice_.CEL.gz    GSM725471_63_HB_Rice_.CEL.gz    GSM725472_9522-1_HB_Rice_.CEL.gz    GSM725473_9522-2_HB_Rice_.CEL.gz    GSM725474_9522-3_HB_Rice_.CEL.gz
AFFX-BioB-3_at   8.37703     8.70610     8.15687     7.74724     8.68766     8.71460
AFFX-BioB-5_at   8.48102     8.74516     8.29362     7.90427     8.75194     8.71383
AFFX-BioB-M_at   8.68125     8.92532     8.41924     8.15696     8.94650     8.96838
AFFX-BioC-3_at  10.12669    10.39479     9.79106     9.47079    10.33159    10.39901
AFFX-BioC-5_at   9.87941    10.20950     9.66993     9.35733    10.10600    10.22251
AFFX-BioDn-3_at 12.50246    12.79332    12.25677    11.86041    12.70061    12.75385
AFFX-BioDn-5_at 11.52168    11.77054    11.33690    10.96891    11.73023    11.80830
AFFX-CreX-3_at  13.88116    13.96916    13.73498    13.42818    13.98805    14.00421
...

使用 limma 包查找差别基因

主要使用上一步获得的 eset

> library(limma)
> pData(eset)
                                 sample
GSM725469_61_HB_Rice_.CEL.gz          1
GSM725470_62_HB_Rice_.CEL.gz          2
GSM725471_63_HB_Rice_.CEL.gz          3
GSM725472_9522-1_HB_Rice_.CEL.gz      4
GSM725473_9522-2_HB_Rice_.CEL.gz      5
GSM725474_9522-3_HB_Rice_.CEL.gz      6
# 能够看到前三个为 mutant 组， 后三个为 WT 组
> group <- c("mutant","mutant","mutant","WTvsmutant","WTvsmutant","WTvsmutant")
> design <- model.matrix(~factor(group))
> colnames(design) <- c("mutant", "WTvsmutant")
> design
# 之因此写"WTvsmutant"而不是"WT",由于这个方法不是按两种分别对待的，即第二组是比较的WTvsmutant, 所以后面的结果 logFC 的值若是是正的，则表明 WT 相对于 mutant 该基因表达上调。

> fit <- lmFit(eset, design)
> fit <- eBayes(fit)
> topTable(fit, coef=2,  n=100, adjust="BH")
# 这样就列出前100 个差别表达的基因列表。

                       logFC AveExpr     t P.Value adj.P.Val   B
Os.50472.2.S1_x_at      4.41     6.2  28.1 7.7e-10   4.4e-05 9.9
OsAffx.17912.1.S1_at    4.12     6.8  22.8 4.6e-09   1.3e-04 9.2
Os.54142.1.S1_at        3.88     5.9  21.2 8.5e-09   1.6e-04 8.9
Os.47752.1.S1_at        2.76     6.1  19.0 2.2e-08   2.4e-04 8.4
OsAffx.28938.1.S1_at    3.80     6.4  18.6 2.6e-08   2.4e-04 8.3
OsAffx.4263.1.S1_at     2.41     5.1  18.4 2.9e-08   2.4e-04 8.2
OsAffx.28280.1.S1_at    1.38     6.2  18.3 2.9e-08   2.4e-04 8.2
Os.7327.4.S1_at         2.16     5.4  17.9 3.6e-08   2.6e-04 8.1
Os.320.2.S1_a_at        2.95     5.8  17.4 4.7e-08   3.0e-04 8.0
...

# 上表中 logFC 表明 log2 以后的fold-change
# AveExpr 表明标准化(log2)以后的基因的表达平均值
# t 表明 t-test
# P.Value 表明p-value, 不过对于基因组水平来讲， p-value 的值不太可靠
# adj.P.Val 表明 q-value , 也能够说是 FDR value, 通常可取 FDR < 0.05 进行过滤
# B  表明 BH 的值

对结果的筛选

若是想提取出全部的差别表达基因，能够在上面 topTable 里把 n 设置为总数
能够把结果写入一个文本中
进行筛选差别基因，参考一些文献，能够按 FDR value < 0.05 而且 fold change >= 1.0 的标准去筛选。
在R中应该能够进行相关的筛选，但对Ｒ不是很熟悉，能够用 python 脚本去实现:

#!/usr/bin/python
# fiter the data

result = open("result.txt")
after_filter = open("after_filter.txt","a")
n = 0

for line in result:   # line is string
    list_line = line.split()  # transform string to a list
    if abs(float(list_line[1])) >= 1 and abs(float(list_line[5])) < 0.05:
        after_filter.write(line)
        n = n + 1

print "the total number of the gene is " + str(n)

result.close()
after_filter.close()

上面这个脚本能够输出一个结果文本，包括了全部符合两个条件的基因，这些基因理论上能够认为是差别基因，能够直接用于后面的各类分析， 好比 cluster , GO, KEGG 等的分析。

以上就是基因芯片分析的简单流程。是很是简化的，还有不少细节和相关理论基础有待加深。