基因芯片数据分析（一）：芯片数据初探
基因芯片概述
简单地讲，基因芯片就是一系列微小特征序列的（通常是DNA探针，也可能是蛋白质）的集合，它们可以被用于定性或者定量检查样品内特异分子的成份。比如说，基因芯片可以检测几十个gene marker在细胞样品中的表达量。现在最常见的是用于整个基因组的表达量分析。它的雏形来自于同位素杂交技术，又如Southern blots或者dot blots。在上世纪九十年代，2维的具有现代意义的基因芯片才在实验室里诞生。基因芯片自问世以来，已经有超过23年（至2014年）了。现在，世界上主流的芯片制造商有4家，分别是Affymetrix，Agilent，Nimblegen以及Illumina。下图为历年来提交至Gene Expression Omnibus数据库的主流芯片厂商的芯片数据统计分布图（数据截止日期为2014年3月1日）。从下图中可以看出，Affymetrix制造的基因芯片在2008年以前占据了市场的主流，在2008年，因为illumina BeadArray的推广，它的市场份额有较大的攀升，但是2年以后就下降至与Affymetrix公司类似的份额。而Agilent却在2010年以后成为芯片市场份额最大的一家。市场份额的变化有价格的因素，质量的因素，使用习惯的因素，也有受到第二代测序技术冲击的因素。
基因芯片的原理大同小异，都是将生物分子固定在某种介质上，比如平面介质，微管，微槽以及微粒等，用于检测目标样品中一批生物分子的含量。使用玻璃介质的好处是经济，稳定，有利于杂交，并且背景荧光小。我们能够最直观想象的基因芯片制造方法就是象打印机印刷一样，通过将来自于不同墨盒里的DNA探针喷印或者吸印至玻璃平板上。在检测的时候，因为每组探针在介质上的位置是已知的，通过不同位置荧光强度与生物分子浓度的对应关系，就不难了解目的片断的表达强度了。这种方法其实就是基因芯片最初的制造方法。这种方法制造的芯片，密度不可能太高。
Affymetrix公司在基因芯片的高密度生产上做出了里程碑式的贡献，它们首创了在石英介质上原位合成oligo短序的工艺，这一工艺极大地提高了基因芯片探针的排布密度。那么oligonucleotide是在原位如何合成的呢？首先，合成位点被对UV光敏感 的化学物质保护起来，当受到UV光照的时候，合成位点被打开，这时候可以把A, C, G或者T碱基连接在合成位点上。而后再去掉光源，位点关闭，洗脱。之后再UV光照打开。如此反复，oligonucleotide的原位合成就完成了。
参考文献:[PMID: 19822891]
芯片数据初探
不同的芯片设计，使得后期数据处理也需要按不同的芯片区别对待。这里将分别针对四家主流的芯片制造商进行简单的代码示例。这里不解释代码的含义，后面我们的教程会讲解，这里只是让大家对基因芯片数据分析有一个初步的认识。
处理Affymetrix的微阵列
## 安装分析所需要的软件包。
BiocManager::install("affy")
BiocManager::install("limma")
BiocManager::install("arrays")
## 载入分析所需要的软件包
library(affy) ## Affymetrix 预处理
library(limma) ## 双色芯片预处理; 差异表达分析
## 从之前安装的arrays软件包中import "phenotype" 数据
## 以及实验设计。
phenoData <-
read.AnnotatedDataFrame(system.file("extdata", "pdata.txt",
package="arrays"))
## RMA normalization
celfiles <- system.file("extdata", package="arrays")# 读取cel文件
eset <- justRMA(phenoData=phenoData,
celfile.path=celfiles)
## differential expression
combn <- factor(paste(pData(phenoData)[,1],
pData(phenoData)[,2], sep = "_"))
design <- model.matrix(~combn) # describe model to be fit
fit <- lmFit(eset, design) # fit each probeset to model
efit <- eBayes(fit) # empirical Bayes adjustment
topTable(efit, coef=2) # table of differentially expressed probesets处理Nimblegen的双色芯片
## 安装分析所需要的软件包，limma也要装，前面装了就不装了
BiocManager::install(c("oligo", "maqcExpression4plex", "genefilter"))
library(oligo)
library(genefilter)
library(limma)
xysPath <- system.file("extdata", package="maqcExpression4plex")
xys.files <- list.xysfiles(xysPath, full.name=TRUE)
basename(xys.files)
theData <- data.frame(Key=rep(c("brain", "universal reference"), each=3))
rownames(theData) <- basename(xys.files)
lvls <- c("exprs", "_ALL_")
vMtData <- data.frame(channel=factor('exprs', levels=lvls),
labelDescription="Sample type")
pd <- new("AnnotatedDataFrame", data=theData, varMetadata=vMtData)
maqc <- read.xysfiles(xys.files, phenoData=pd)
class(maqc)
eset <- rma(maqc)
design <- model.matrix(~factor(eset[["Key"]]))
fit <- lmFit(eset, design)
efit <- eBayes(fit)
topTable(efit, coef=2)
##pair双色文件
BiocManager::install("Ringo")
library(Ringo)
exDir <- system.file("exData",package="Ringo")
list.files(exDir, pattern="pair.txt")
head(read.delim(file.path(exDir,"MOD_20551_PMT1_pair.txt"),
skip=1))[,c(1,4:7,9)]
read.delim(file.path(exDir,"example_targets.txt"), header=TRUE)
read.delim(file.path(exDir,"spottypes.txt"), header=TRUE)
exRG <- readNimblegen("example_targets.txt","spottypes.txt",path=exDir)
exampleX <- preprocess(exRG)
sampleNames(exampleX) <-
with(exRG$targets, paste(Cy5,"vs",Cy3,sep="_"))
print(exampleX)
load(file.path(exDir,"exampleProbeAnno.rda")) ## see details in retrieveGenomicFeatureAnnotation.R
smoothX <- computeRunningMedians(exampleX, probeAnno=exProbeAnno,
modColumn = "Cy5", allChr = "9",
winHalfSize = 400)
sampleNames(smoothX) <- paste(sampleNames(exampleX),"smoothed")
(y0 <- apply(exprs(smoothX),2,upperBoundNull))
chersX <- findChersOnSmoothed(smoothX, probeAnno=exProbeAnno, thresholds=y0,
allChr="9", distCutOff=600, cellType="human")
chersX <- relateChers(chersX, exGFF)
chersXD <- as.data.frame.cherList(chersX)
chersXD[order(chersXD$maxLevel, decreasing=TRUE),]处理Agilent的txt文件
##下载示例数据
library(GEOquery)
GSE_NUM <- "GSE29769"
gset <- getGEO(GSE_NUM)[[1]]
pd <- pData(gset)
pd.c <- apply(pd, 2, function(.ele) length(unique(.ele)))
pd.s <- pd[, pd.c>1] ## 下载得到数据描述
getGEOSuppFiles(GSE_NUM)
setwd(GSE_NUM)
dir()
untar(paste0(GSE_NUM, "_RAW.tar"))
files <- dir(pattern="gz$")
sapply(files, gunzip)
library(limma)
BiocManager::install("RGList")
library(RGList)
BiocManager::install("backgroundCorrect")
library(backgroundCorrect)
filelist <- sub(".gz", "", basename(as.character(pd.s$supplementary_file)))
rownames(pd.s) <- filelist
data <- read.maimages(files=filelist, source = "genepix")
RG <- backgroundCorrect(data, method = "normexp") # 代码有问题
MA <- normalizeWithinArrays(RG)
f <- sub("_rep.", "", pd.s$title)
f <- factor(f)
design <- model.matrix(~-1+f)
colnames(design) <- sub("^f", "", colnames(design))
constr <- c("Adult_mouse_heart_infarcted-Adult_mouse_heart_control")
contrast.matrix <- makeContrasts(contrasts=constr, levels=design)
fit <- lmFit(MA, design)
fit1 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit1)
topTable(fit2)处理Illumina的BeadArray
BiocManager::install(c("beadarray", "beadarrayExampleData", "illuminaHumanv3.db"))
library(beadarray)
require(beadarrayExampleData)
data(exampleSummaryData)
exampleSummaryData.unlogged <- channel(exampleSummaryData, "G.ul")
BSData <-
addFeatureData(exampleSummaryData.unlogged,
toAdd = c("SYMBOL", "PROBEQUALITY",
"CODINGZONE", "PROBESEQUENCE", "GENOMICLOCATION"))
BSData.norm <- normaliseIllumina(BSData,
method="quantile",
transform="log2")
qual <- fData(BSData.norm)$PROBEQUALITY
rem <- qual == "No match" | qual == "Bad" | is.na(qual)
BSData.norm <- BSData.norm[!rem, ]
dim(BSData.norm)
rna <- factor(pData(BSData.norm)[,"SampleFac"])
design <- model.matrix(~-1+rna)
colnames(design) <- levels(rna)
design
library(limma)
aw <- arrayWeights(exprs(BSData.norm), design)
fit <- lmFit(exprs(BSData.norm), design, weights=aw)
contrasts <- makeContrasts(UHRR-Brain, levels=design)
contrasts
contr.fit <- eBayes(contrasts.fit(fit, contrasts))
colnames(contr.fit$coefficients)
res <- topTable(contr.fit, coef="UHRR - Brain", number=nrow(contr.fit))
features <- fData(BSData.norm)
res <- cbind(res, features[rownames(res), ])
head(res)作者水平有限，欢迎批评指正，互相进步！
本文分享自 MedBioInfoCloud 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体分享计划 ，欢迎热爱写作的你一起参与！