搜索
查看: 3014|回复: 0

[software] oligo 阵列数据读取和归一化

[复制链接]

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-7-9 09:22:54 | 显示全部楼层 |阅读模式
## Overview
  
本文档介绍了用于处理Affymetrix和Nimblegen微阵列的`r Biocpkg(“oligo”)包,特别是基因表达,外显子表达和SNP阵列数据。

## Dependencies

This document has the following dependencies:

[AppleScript] 纯文本查看 复制代码
```{r dependencies, warning=FALSE, message=FALSE}
library(oligo)
library(GEOquery)
```


Use the following commands to install these packages in R.

[AppleScript] 纯文本查看 复制代码
```{r biocLite, eval=FALSE}
source("http://www.bioconductor.org/biocLite.R")
biocLite(c("oligo", "GEOquery"))
```

## Getting the data

我们将使用保存为GEO登录号“GSE38792”的数据集。 在这个数据集中,实验者从两个不同的条件分析脂肪活检:10例阻塞性睡眠呼吸暂停(OSA)和8例健康对照。

使用Affymetrix人类基因ST 1.0阵列进行分析。

首先我们需要获取原始数据; 这将是一组CEL格式的二进制文件。 每个样本将有一个文件。 CEL文件可作为GEO的补充信息访问; 我们使用`r Biocpkg(“GEOquery”)`获取文件。

[AppleScript] 纯文本查看 复制代码
```{r getData}
library(GEOquery)
getGEOSuppFiles("GSE38792")
list.files("GSE38792")
untar("GSE38792/GSE38792_RAW.tar", exdir = "GSE38792/CEL")#utils包中的命令
list.files("GSE38792/CEL")
```


`r Biocpkg(“oligo”)`和许多其他这样的软件包可以方便地读取许多文件。 在这种情况下,我们构造一个文件名向量,并将其提供给`read.celfiles()`。

[AppleScript] 纯文本查看 复制代码
```{r readData, message=FALSE}
library(oligo)
celfiles <- list.files("GSE38792/CEL", full = TRUE)
rawData <- read.celfiles(celfiles)
```


[AppleScript] 纯文本查看 复制代码
```{r show}
rawData
```


这是一个“GeneFeatureSet”的形式; 这是一个类似“ExpressionSet”container。 根据S4,我们可以通过类定义看到这一点
[AppleScript] 纯文本查看 复制代码
```{r getClass}
getClass("GeneFeatureSet")
```


我们看到这是一个特殊情况的“FeatureSet”,这是一个特殊情况的“NChannelSet”,它是一个“eSet”。 我们可以看到强度被衡量:
[AppleScript] 纯文本查看 复制代码
```{r rawPeak}
exprs(rawData)[1:4,1:3]
```


我们看到这是原始强度数据; 测量单位是16位扫描仪,所以我们得到的值在0和$ 2 ^ 16 = 65,536 $之间。 这很容易验证:
[AppleScript] 纯文本查看 复制代码
```{r maxExpr}
max(exprs(rawData))
```


注意这个数据集中的大量特征,超过100万。 由于技术限制,Affymetrix只能制造非常短的寡核苷酸(约25bp),但其廉价,高品质。 短寡核苷酸意味着寡核苷酸的结合特异性不是很好。 为了弥补这一点,Affymetrix使用一种设计,其中同时通过许多不同的探针测量基因; 这被称为探针组。 作为Affymetrix阵列的预处理步骤的一部分,探针组中所有探针的测量需要合并成一个表达式测量。

让我们去除`rawData`的表型信息。
[AppleScript] 纯文本查看 复制代码
```{r pData}
filename <- sampleNames(rawData)
pData(rawData)$filename <- filename
sampleNames <- sub(".*_", "", filename)
sampleNames <- sub(".CEL.gz$", "", sampleNames)
sampleNames(rawData) <- sampleNames
pData(rawData)$group <- ifelse(grepl("^OSA", sampleNames(rawData)),
                               "OSA", "Control")
pData(rawData)
```


## Normalization

让我们来看一下样品的探针强度,使用`boxplot()'函数。
[AppleScript] 纯文本查看 复制代码
```{r rawBox, plot=TRUE}
boxplot(rawData)
```


Boxplots非常适合用于比较许多样品,因为它很容易并排显示许多盒子图。我们看到样本之间的位置和分布差异很大。有三个样品,强度很低;几乎所有探针在log2尺度上的强度都小于7。从Affymetrix微阵列的经验来看,可知道这是一个非常低的强度。也许这些阵列的阵列杂乱失败。要确定这将需要更多的调查。

用于预处理Affymetrix基因表达阵列的经典和强大的方法是RMA方法。经验告诉我们,RMA基本上总是表现良好,所以很多人喜欢这种方法;人们可以认为最好使用一种总是做得很好的方法,而不是在某些数据集上做得非常好,而对其他数据集不好的方法。

RMA方法最初在`r Biocpkg(“affy”)`package中实现,后来被`r Biocpkg(“oligo”)`包取代。我们正在分析的数据来自于基于随机启动的“新”样式Affymetrix阵列; `r Biocpkg(“affy”)`package不支持这些类型的数组。运行RMA非常简单:

[AppleScript] 纯文本查看 复制代码
```{r rma}
normData <- rma(rawData)
normData
```


请注意,“normData”对33k特征的数量是否更接近于人类基因组中的基因数量。

我们可以再次查看箱形图来检查RMA的性能。

[AppleScript] 纯文本查看 复制代码
```{r normBox, plot=TRUE}
boxplot(normData)
```



在这里,重要的是要记住,第一组框图位于探测级(〜1M探针),而第二组框图位于探针组级(〜33k探针组),因此它们以不同的摘要级别显示数据。 然而,重要的分析是,探针分布在样本之间进行归一化,乍一看它看起来很好。 可以看出,3个可疑样本之前仍然略有不同,但至少再有两个样本与之相似。

对于归一化感兴趣的人,注意尽管分布是相似的,但它们并不相同,尽管RMA包括分位数归一化。 这是因为分位数归一化在探测摘要之前完成; 如果您将分位数归一化为不同的分布,则它们将保证以后具有相同的分布。

归一化后,数据 可以用于差异表达分析。

回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-11-20 13:26 , Processed in 0.028424 second(s), 23 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.