基因芯片分析是生物信息学的重要应用之一，今天总结其中一种简单的使用方法。

一. 获取数据

样本数据选取了论文Molecular subsets of mantle cell lymphoma defined by the IGHV mutational status and SOX11 expression have distinct biologic and clinical features.的实验数据，该套数据的下载网址：
https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE36000&format=file
从NCBI的GEO数据库获取数据的全部信息
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36000

这张芯片主要是研究来自Mantle细胞淋巴瘤（MCL）患者的B细胞淋巴细胞的免疫球蛋白重链可变（IGHV）基因的突变（MUT）与非突变(UNMUT)
芯片共38张，其中24张为突变型（MUT），14张为非突变型（UNMUT）

二. 数据预处理

过程梳理
质量检测——背景处理、归一化、汇总——探针注释——合并重复探针

1.质量检验

数据读入

#初次使用，下载 CLL package
#library(BiocInstaller)
#biocLite("CLL")
#biocLite("simpleaffy")
#install.packages("scales")

library(affy)
library(tcltk)
library(gcrma)
library(RColorBrewer)
library(simpleaffy)
library(CLL)

setwd("E:/生信/大三/转录组学")
geneCELs=list.celfiles('E:/生信/大三/转录组学/GSE36000_RAW/',full.name=T)
data.raw<-ReadAffy(filenames=geneCELs)
#basename(geneCELs)     #查看文件名
#由于数据分布本身极不规则，这里通过sublime text的优质工具便捷导入名称
treats <- strsplit("UNMUT UNMUT UNMUT UNMUT UNMUT MUT MUT UNMUT UNMUT UNMUT MUT UNMUT MUT MUT MUT MUT MUT MUT MUT MUT MUT MUT UNMUT MUT MUT MUT MUT MUT MUT UNMUT UNMUT MUT MUT MUT MUT MUT UNMUT UNMUT", " ")[[1]]
snames <- strsplit("UNMUT1 UNMUT2 UNMUT3 UNMUT4 UNMUT5 MUT1 MUT2 UNMUT6 UNMUT7 UNMUT8 MUT3 UNMUT9 MUT4 MUT5 MUT6 MUT7 MUT8 MUT9 MUT10 MUT11 MUT12 MUT13 UNMUT10 MUT14 MUT15 MUT16 MUT17 MUT18 MUT19 UNMUT11 UNMUT12 MUT20 MUT21 MUT22 MUT23 MUT24 UNMUT13 UNMUT14", " ")[[1]]
sampleNames(data.raw) <- snames

step0.首先对芯片进行质量检测(是否存在spatial artifact)

for (i in 1:length(sampleNames(data.raw))){
  name = paste("image",sampleNames(data.raw)[i],".jpg",sep="")
  jpeg(name)
  image(data.raw[,i])
  dev.off()
}

下面为MUT1芯片，所有芯片均完好，无spatial artifact

2.标准化

step1.背景处理、归一化：MAS5/rma/gcrma

# mas5
data.raw.mas5 <- mas5(data.raw)
data.mas5.Matrix <- exprs(data.raw.mas5)    #
# rma
data.raw.rma <- rma(data.raw)
data.rma.Matrix <- exprs(data.raw.rma)
# gcrma
data.raw.gcrma <- gcrma(data.raw)
data.gcrma.Matrix <- exprs(data.raw.gcrma)
## boxplot
boxplot(data.mas5.Matrix, las=3, ylim=c(0,600), main="MAS 5.0")
boxplot(data.rma.Matrix, las=3, ylim=c(0,8), main="RMA")
boxplot(data.gcrma.Matrix, las=3, ylim=c(0,8), main="gcRMA")

RMA拖尾情况最少，相对而言最好

3.过滤探针

calls <- mas5calls(data.raw) # 获得探针的PMA信息
calls <- exprs(calls) # 提取包含探针的PMA信息的表达矩阵
head(calls)
present <- rowSums(calls == "P") # 统计每个探针在所有样本中分类为"P"的个数
present <- which(present != 0) # 提取至少在一个样本中表达的探针
rmaFiltered <- data.raw.rma[present,] # 留下至少在一个样本中表达的探针
#nrow(rmaFiltered)
#Features 
#   38652 

4.探针注释

library(GEOquery)
gse36<-getGEO('GSE36000',destdir =".") #根据GSE号来下载数据，下载_series_matrix.txt.gz
gpl570<-getGEO('GPL570',destdir =".")  #根据GPL号下载的是芯片设计的信息, soft文件

colnames(Table(gpl570))
 #提取探针注释信息，保留第1列和第4列
genename <- Table(gpl515)[,c(1,4)]
exprSet <- as.data.frame(exprs(gse382[[1]]))# 得到表达矩阵，行名为ID，需要转换
#转换ID为gene name
exprSet$ID <- rownames(exprSet)
express <- merge(x = genename, y = exprSet, by = "ID")
express$ID =NULL

5. 合并重复探针

#安装hgu133plus2.db包
#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("hgu133plus2.db", version = "3.8")

library(hgu133plus2.db) #加载芯片注释包
columns(hgu133plus2.db) #查看芯片注释包提供的注释信息
probe2symbol <- toTable(hgu133plus2SYMBOL) #提取注释信息，这里选择GENESYMBOL进行注释。也可以选择其他ID进行注释，详情请见注释包的技术支持文档。
exprSet <- as.data.frame(exprs(rmaFiltered)) #提取表达矩阵，转换为data.frame对象
head(exprSet) #查看表达矩阵

exprSet$probe_id <- rownames(exprSet) #将行名转化为probe_id

probeset <- rownames(exprSet) #提取行名为探针名
Symbol <- as.character(as.list(hgu133plus2SYMBOL[probeset])) #提取探针对应的GENESYMBOL
exprSet_symbol <- cbind(Symbol,exprSet)

#采取数据中位数代替重复值
exprSet_symbol <- aggregate(x = exprSet_symbol[,-1],by = list(exprSet_symbol$Symbol), FUN = median)
exprSet_symbol$probe_id <- NULL
dim(exprSet_symbol)
#[1] 15672    39

注意，此处的生成数据是经过RMA修正的，因此是求log后的值

三. 基因的差异表达分析

将数据处理成纯数据

exprSet_symbol_1<-exprSet_symbol[,-1]
rownames(exprSet_symbol_1)<-exprSet_symbol[,1]

1. t.test

unmut<-exprSet_symbol_1[,c(1,2,3,4,5,8,9,10,12,23,30,31,37,38)]
mut<-exprSet_symbol_1[,c(6,7,11,13,14,15,16,17,18,19,20,21,22,24,25,26,27,28,29,32,33,34,35,36)]



for(i in 1:nrow(exprSet_symbol_1)){
    p_value[i]<-t.test(unmut[i,],mut[i,])$p.value
}
exprSet_symbol_1$p_value<-p_value
exprSet_symbol_p_005 <- exprSet_symbol_1[p_value<0.05,]
#dim(exprSet_symbol_p)
#[1] 2503   39
exprSet_symbol_p_001 <- exprSet_symbol_1[p_value<0.01,]
#dim(exprSet_symbol_p)
#[1] 969  39


exprSet_symbol_unlog<-2^exprSet_symbol_1
unmut<-exprSet_symbol_unlog[,c(1,2,3,4,5,8,9,10,12,23,30,31,37,38)]
mut<-exprSet_symbol_unlog[,c(6,7,11,13,14,15,16,17,18,19,20,21,22,24,25,26,27,28,29,32,33,34,35,36)]

2.倍数法

unmut_avg<-apply(unmut,1,mean)
mut_avg<-apply(mut,1,mean)
exprSet_symbol_1$unmut_avg<-unmut_avg
exprSet_symbol_1$mut_avg<-mut_avg
exprSet_symbol_fold_2 <- exprSet_symbol_1[2^(log2(unmut_avg/mut_avg))>2,]
# dim(exprSet_symbol_fold)
# [1] 407  41
exprSet_symbol_fold_3 <- exprSet_symbol_1[2^(log2(unmut_avg/mut_avg))>3,]
# dim(exprSet_symbol_fold_3)
# [1] 145  43

3.z值

exprSet_symbol_1$avg_cha<-abs(unmut_avg-mut_avg)
exprSet_symbol_z<-vector(length=nrow(exprSet_symbol_1))
for(i in 1:nrow(exprSet_symbol_1)){
    exprSet_symbol_z[i]<-(exprSet_symbol_1$avg_cha[i]-mean(exprSet_symbol_1$avg_cha))/sd(exprSet_symbol_1$avg_cha)
}
exprSet_symbol_1$exprSet_symbol_z<-abs(exprSet_symbol_z)
exprSet_symbol_zscore_196<-exprSet_symbol_1[exprSet_symbol_z>1.96,]
# dim(exprSet_symbol_zscore)
# [1] 305  43
exprSet_symbol_zscore_258<-exprSet_symbol_1[exprSet_symbol_z>2.58,]
# dim(exprSet_symbol_zscore_258)
# [1] 224  43

4.筛选

#t 0.05  倍数 2
exprSet_symbol_t005_fold2<-exprSet_symbol_fold_2[exprSet_symbol_fold_2$p_value<0.05,]
# dim(exprSet_symbol_t005_fold2)
# [1] 283  43

#t 0.05  倍数 3
exprSet_symbol_t005_fold3<-exprSet_symbol_fold_3[exprSet_symbol_fold_3$p_value<0.05,]
# dim(exprSet_symbol_t005_fold3)
# [1] 116  43

#t 0.05  z 0.05
exprSet_symbol_t005_z005<-exprSet_symbol_p_005[exprSet_symbol_p_005$exprSet_symbol_z>1.96,]
# dim(exprSet_symbol_t005_z005)
# [1] 174  43

#t 0.05  z 0.01
exprSet_symbol_t005_z001<-exprSet_symbol_p_005[exprSet_symbol_p_005$exprSet_symbol_z>2.58,]
# dim(exprSet_symbol_t005_z001)
# [1] 141  43

（P=0.01&&倍数法3）、（P=0.01&&Z值法0.05）

exprSet_symbol_t001_fold3<-exprSet_symbol_fold_3[exprSet_symbol_fold_3$p_value<0.01,]
exprSet_symbol_t001_fold3_z005<-exprSet_symbol_t001_fold3[exprSet_symbol_t001_fold3$exprSet_symbol_z>1.96,]
# dim(exprSet_symbol_t001_fold3_z005)
# [1] 13 43