专栏首页生物信息云基因芯片数据分析(一):芯片数据初探

基因芯片数据分析(一):芯片数据初探

关于芯片数据分析,我们之前的文章

基因芯片概述

简单地讲,基因芯片就是一系列微小特征序列的(通常是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(MedBioInfoCloud),作者:DoubleHelix

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-10-18

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 分子对接教程 | (9) VMD可视化对接结果

    能够实现蛋白质三维结构可视化的软件非常多。比专业级的PyMOL(https://pymol.org/2/)。这个软件已经被世界上著名的生物医药软件公司“薛定谔公...

    DoubleHelix
  • R语言数据分析与挖掘(第七章):因子分析

    因子分析(factor analysis, 简称FC)又称因素分析,基于相关关系而进行的数据分析技术,是一种建立在众多的观测数据的基础上的降维处理方法。其主要目...

    DoubleHelix
  • R语言数据分析与挖掘(第一章):数据预处理(3)——数据整理

    在介绍了缺失值处理的方法之后,我们可以得到完整的数据集,但在进行数据分析之前,还需要对数据进行整理,下面我们将介绍数据整理的相关知识。

    DoubleHelix
  • Redis 基本数据结构三:哈希

    几乎所有的编程语言都提供了哈希(hash)类型,例如 Java 中的 Map,python 中的字典,在Redis中,哈希类型是指键的值本身又是一个键值对结构,...

    CoderJed
  • 关于动态重载Lua脚本的一些思考

    平时工作中自己多使用 Lua 脚本,过程中常常会遇到一个痛点:如何动态重载Lua脚本以加快开发的迭代速度.

    用户2615200
  • Infor ERP LN中采购订单状态说明

    崔文远TroyCui
  • shell 基本语法

    jenkins 上构建项目时,经常需要借助 shell 脚本,最近也经常跟服务器打交道,顺便记录些常用命令,方便查阅

    请叫我大苏
  • Linux下使用python脚本执行BCP导入导出操作

    以上python脚本首先从test002中将数据查询出来,将结果集使用BCP写入tempData.csv文件中,然后再使用BCP将文件中的数据写入表test00...

    我是李超人
  • 用Single-spa 创建基于 React 和 Vue 的微型前端

    它使你可以在单页应用中使用多个框架,这样就可以按功能拆分代码,并 能使 Angular、React、Vue.js 程序一起运行。

    疯狂的技术宅
  • SpringBoot2 整合 ClickHouse数据库,实现数据高性能查询分析

    知了一笑

扫码关注云+社区

领取腾讯云代金券