基因芯片的簡單分析

本文主要介紹一下基因芯片的簡單分析過程。

文檔參考

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 等的分析。

以上就是基因芯片分析的簡單流程。是很是簡化的，還有不少細節和相關理論基礎有待加深。