基因芯片的简单分析
本文主要介绍一下基因芯片的简单分析过程。
文档参考
Using Bioconductor for Microarray Analysis
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 等的分析。
以上就是基因芯片分析的简单流程。是非常简化的，还有很多细节和相关理论基础有待加深。
