导读

个人笔记，初次分析芯片数据，纯属练习，不涉及生物学意义。

一、数据获取

1 GEO官网中输入GSE编号，检索

2 获取分组信息，下载丰度表矩阵、平台探针信息等数据

3 解压数据如下：

说明：
GPL96-tbl-1.txt：含探针<=>基因对应关系
GSE14359_family.soft：含探针<=>基因对应关系
GSE14359_series_matrix.txt：探针、样品信号值矩阵
GSM*****-tbl-1.txt：芯片/样品探针信号值

二、制作分组文件

cat group.txt 
FileName    Target
GSM359137   nnpoc
GSM359138   nnpoc
GSM359139   cot
GSM359140   cot
GSM359141   olmt
GSM359142   olmt
GSM359143   cot
GSM359144   cot
GSM359145   olmt
GSM359146   olmt
GSM359147   cot
GSM359148   cot
GSM359149   cot
GSM359150   cot
GSM359151   olmt
GSM359152   olmt
GSM359153   olmt
GSM359154   olmt
GSM359155   cot
GSM359156   cot

三、limma包做探针的差异值分析

1 数据准备

library(limma)
options(digits=3)  # 全局，保留三位小数

# 读取文件：矩阵、分组
data = read.table("GSE14359_series_matrix.txt", comment.char="!", header=T, sep="\t", row.names=1)
group = read.table("group.txt", header=T, sep="\t")

# 分组文件处理：排序，矩阵化，因子化
group = group[order(group$Target, decreasing=T),]
category = model.matrix(~0 + factor(group$Target, levels=unique(group$Target)))
colnames(category) = unique(group$Target)

category

2 选两组，lmFit线性模型拟合

# 测试：选两组，lmFit线性模型拟合
category_1 = makeContrasts("olmt-nnpoc", levels = category)
fit = lmFit(data, category)
fit2 = contrasts.fit(fit, category_1)
fit3 = eBayes(fit2)

3 topTable求差异基因

result = topTable(fit3, adjust.method="BH", sort.by="logFC", n=20, p.value=0.001, lfc=2, confint=FALSE)

参数：
lfc: min fold change
adjust.method: "none", "BH", "BY" and "holm"
sort.by: "logFC", "AveExpr", "t", "P", "p", "B" or "none"
n=Inf: output for all probes in original (unsorted) order
confint: should confidence 95% intervals be output for logFC
number: maximum number of genes to list [与n=Inf冲突]

4 结果

result_select = subset(result, select=c("adj.P.Val","P.Value","logFC"))
colnames(result_select) = c("FDR", "P.Value", "logFC")
head(result_select)

说明
第一列：探针
第二列：FDR矫正P
第三列：log(fold chang value)

待续：探针-基因的匹配用merge可以搞定，一对多的关系可能需要取平均值。

参考：
R语言limma包差异基因分析