一、 illumina磁珠芯片制造原理

芯片由玻璃基片和微珠组成，光蚀刻出许多微米级小孔，用于容纳微珠。每个微珠表面偶联一种序列的DNA片段（一个珠子上片段序列相同），每个微珠上有几十万个片段。

5'端address序列是标识微珠的标签序列，每种序列就是微珠的身份证号ID。

3'端probe为探针序列

illumina芯片的制造过程：把要做芯片的几十万种微珠，按设定的比例混合好，撒到玻璃基片上，微珠随机的落入基片的小孔中。然后通过检测芯片上每个小孔中微珠的address序列，就可以知道小孔当中是哪种微珠。因为address序列和probe序列一一对应，因此就知道每个小孔中有那种probe。

因此illumina公司出厂的每张芯片都跟一个.dmap文件，dmap文件标注了每一张芯片的每一个微孔中分别是哪种微珠。用户扫描完芯片后，要从Illumina网站上下载这张芯片对应的dmap文件然后才能解读芯片。在一张芯片的一个反应中，每种珠子平均有15颗或更多（表达量芯片中有30个左右）。

❝

详见：bilibili陈巍学基因，这段引自 "Illumina的SNP芯片原理" 。

❞

二、 lumi: a pipeline for processing Illumina microarray

「摘要：」 illumina使用的BeadArray（磁珠芯片）技术需要不同的预处理和质控，相较于其他芯片技术。不幸的是，大多数分析软件并没有利用BeadArray 系统的独特特性，而只是结合了最初为 Affymetrix 芯片设计的预处理方法。lumi是专门为处理illumina芯片数据设计的R包，可以从Bioconductor下载获得。它包括芯片读入，质控，固定方差，标准化和基因注释部分。具体来说，lumi包 包括「方差固定变换 (variance-stabilizing transformation，VST) 算法」，该算法利用每个 Illumina 芯片上可用的技术性重复（同一序列对应约30个磁珠）。该包提供了不同的标准化方法和众多的质控图。为了更好的注释illumina数据，供应商自主的核苷酸通用标识符（nuID）用于识别illumina芯片的探针。nuID注释包和lumi处理输出的结果可以轻松地与其他Bioconductor包集成，以构建Illumina数据的统计分析的工作流程。

1 介绍

illumina磁珠芯片有约30个随机定位的「重复磁珠」（具有同样的探针序列）。与其他类型的芯片相比，这种额外的设计可产生更高的置信度和更稳健的估计。然而，Illumina 微阵列设计的独特性使得预处理和质量控制步骤与其他类型的微芯片显著不同。不幸的是，到目前为止（2007），大多数分析软件都没有利用 Illumina BeadArray 系统的这种独特特性，而只是结合了最初为 Affymetrix 芯片设计的预处理方法。lumi包提供专门为illumina磁珠芯片设计的算法，并且使用现有的算法和基因注释的框架。除了支持芯片数据的现有算法外，lumi 包还包括几个独特的部分：(1) 利用 Illumina 芯片上可用的技术重复的固定方差变换 (VST)；(2) 为 Illumina 微阵列数据设计的标准化算法 ；（3) 核苷酸通用标识符 (nuID) 注释包。nuID 注释包允许对每个探针进行不依赖版本和供应商的注释。nuID 还通过包括错误检查的过程对原始探针序列进行唯一且准确的编码。

2 应用

2.1 预处理和质控

lumi包包含一个主要类：「LumiBatch」，它继承自Bioconductor中的 「ExpressionSet类」，以实现与其他Bioconductor包的互操作性。LumiBatch类的结构如下所示，它在ExpressionSet类的基础上延伸出了三个元素：在assayData槽下的：se.exprs, beadNum 和 detection，用于存储illumina磁珠芯片的额外信息。

controlData槽保存对照探针的信息，QC槽保存质控总结，history槽追踪所有在LumiBatch对象上进行的操作，可以解释数据来源。用户可以选择将BeadStudio输出的Illumina注释信息保留在LumiBach对象的featureData中。

lumi包中有几种主要的处理方式。lumiR 通过智能读取所有版本的 Illumina BeadStudio 软件的原始数据来初始化 LumiBatch 对象，并且 lumiR.batch 方法旨在读取一批数据文件。
lumiB调整芯片背景；lumiT 对数据进行方差固定；lumiN 对方差固定后的数据进行标准化，lumiQ 评估数据质量。所有这些方法构成了一个预处理流程。为了方便起见，「lumiExpresso 函数」将所有四种方法封装为一个。

这些方法还包括调用先前为 Affymetrix 数据设计的其他处理方法的选项。出于为更好的可视化和质控的目的，lumi 包还提供了不同类型的可视化功能。这些绘图函数可以处理表达和对照探针数据。更多详细信息请参阅教程和函数帮助文件。

2.2 注释包

Illumina 注释包是使用 Bioconductor 注释工具构建的，并使用每个探针的 nuID 作为标识符。比对 TargetID 或 ProbeId 到 nuID 也包含在注释包中。由于Illumina的芯片均采用50mers，通过使用nuID通用标识符，我们可以为同一物种的不同版本芯片建立一个注释数据库。此外，nuID可以直接转换为探针序列，并用于获取最新的refSeq匹配和注释。目前所有Illumina表达芯片的注释包（包名称以“lumi”为前缀，后跟物种名称和版本号，例如lumiHumanAll.db）可以从Bioconductor下载。

3 使用案例

图2 显示数据处理流程图。用于预处理的R源代码如图3所示。由于lumi包中的类是从类ExpressionSet扩展而来的，因此Bioconductor中的许多数据分析包可以直接应用于lumi产生的结果。图 2 显示了使用 lumi 包以及 limma、GOstats 和 MLInterfaces 的场景。更多的实现细节可以在lumi包的教程中找到。

总之，lumi 包提供了LumiBatch类的基础框架和相关方法来构建 Illumina 从原始数据开始到功能分析的工作流程。

三、 以GSE67936为例的原始数据处理

可以看到是illumina磁珠的表达量芯片，一共有168个样品。

Supplementary files 中有一个RAW原始数据的压缩包和一个non-normalized为标准化数据的压缩包。这个示例数据中的RAW.tar不可用，存储的平台的注释信息。我们读取"GSE67936_non-normalized.txt" 进行分析。

rm(list = ls())  #一键清空
options(stringsAsFactors = F)

GSE_number<-"GSE67936"

library(GEOquery)
library(stringr) 

#从GEO下好non-normalized.txt.gz
#解压
rawdata <- paste0(GSE_number,"_non-normalized.txt")
series_matrix<-paste0(GSE_number,"_series_matrix.txt.gz")

#下载好series_matrix，用现成文件导入获取临床信息
getGEO(file= series_matrix , destdir = "./",getGPL = F)->geoObject

#提取表达矩阵和临床信息
pData(geoObject[1])->pd

#读取non-normalized_data.txt
library(lumi)

a=read.table(rawdata,header = T,sep = '\t')
colnames(a);ncol(a)

读进来的a中第一列是探针id，第二列是symbol，从第三列起每两列对应一个样本的信号值和pvalue。

整理矩阵行名使之适用于lumiR的输入：

#168个样本这里需要根据上面代码返回结果自己改动，前两行是探针id和对应的symbol
#后面每两列对应一个样本，是荧光信号强度和pvalue?
colnames(a)=c('PROBE_ID',"Gene",paste(rep(rownames(pd), each =2),
                               rep(c('AVG_Signal','Detection Pval'), 168)
                               ,sep='.'))
#如果有probe_id 和gene_symbol的话，就提取出来做ids，没有去GEO下载GPL表格做探针转化
ids<-data.frame(a$PROBE_ID,a$Gene)
head(colnames(a))
write.table(a,file = 'tmp.txt',sep = '\t',quote = F)

之后按照lumi文献提供的方法lumiR读入数据。

lumiB调整芯片背景；lumiT 对数据进行方差固定；lumiN 对方差固定后的数据进行标准化，lumiQ 评估数据质量。这些步骤可以使用封装好的lumiExpresso函数一步完成

x.lumi <- lumiR('tmp.txt',lib.mapping = 'lumiHumanIDMapping')
lumi.N.Q <- lumiExpresso(x.lumi)
class(lumi.N.Q)

因为这里我们处理过的lumi.N.Q已经是ExpressionSet对象了，所以可以用exprs函数提取其表达矩阵

可以看到行名既不是探针id也不是symbol，还有补救方法，从lumi.N.Q@featureData@data中提取探针id，在转换成symbol

## retrieve normalized data
dat <- exprs(lumi.N.Q)
probeid<-lumi.N.Q@featureData@data
probeid<-probeid[match(rownames(dat),rownames(probeid)),]
identical(rownames(dat),rownames(probeid))
probeid$PROBE_ID->rownames(dat)
dat[1:4,1:4]

#有些探针对应不到基因，去掉symbol为空的
colnames(ids)<-c("probeid","symbol")
ids=ids[ids$symbol!= '',]
dat=dat[rownames(dat) %in% ids$probeid,]
ids=ids[match(rownames(dat),ids$probeid),]
head(ids)  
#在这个示例数据集中probeid的组成有字母有数字，有些只有数字的就要转成character
ids$probeid=as.character(ids$probeid)
identical(rownames(dat),ids$probeid)

ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果
dat=dat[ids$probeid,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4]  #保留每个基因ID第一次出现的信息

箱线图检查一下单个样本表达量分布和样本间方差齐性

png("check_boxplot.png")
boxplot(dat,las=2)
dev.off()
limma::normalizeBetweenArrays(dat)->dat
png("check_normalized_boxplot.png")
boxplot(dat,las=2)
dev.off()
save(dat,pd,file="step0_output.R")

之后就可以进行后续分析了