illumina beadchip 芯片原始数据处理

需求

和这一篇的需求一致，不同公司的芯片原始数据处理方式也不同，最终都是为了得到表达矩阵，今天介绍的是illumina beadchip 芯片，用神奇的limma搞定它。
找不到某些GEO数据的表达矩阵肿么办
illumina beadchip芯片原始数据处理,参考limma的userguide。

1.帮助文档里的示例

rm(list=ls())
library(limma)
x <- read.ilmn(files="probe profile.txt",
               other.columns="Detection")

## Reading file probe profile.txt ... ...

x$E[1:4,1:4]
##                     1        2        3        4
## ILMN_1762337 52.34406 46.10429 54.01238 47.65873
## ILMN_2055271 69.91481 73.86729 58.64280 72.36581
## ILMN_1736007 57.47208 53.70050 53.39117 49.43674
## ILMN_2383229 53.60817 57.50813 48.22960 48.18757

读取原始文件，这么简单就拿到表达矩阵了，原示例中有一个control probe profile.txt，文档说了它不是必须的，我直接去掉咯。画个箱线图

boxplot(log2(x$E),range=0,ylab="log2 intensity")

image.png

做背景校正、标准化

y <- neqc(x) 
## Note: inferring mean and variance of negative control probe intensities from the detection p-values.

听起来很厉害，实际上一个函数搞定。

探针过滤

原始文件里的信息，被拆分成了一个表达矩阵和一个Detection P值矩阵，P值可以用于过滤样本。

示例数据总共有12个样本，过滤标准是至少在3个样本里P值<0.05。

x$other$Detection[1:4,1:4]

##                       1         2          3           4
## ILMN_1762337 0.55849580 0.6754875 0.13698630 0.601388900
## ILMN_2055271 0.03064067 0.0000000 0.04931507 0.002777778
## ILMN_1736007 0.27715880 0.2924791 0.15342470 0.486111100
## ILMN_2383229 0.47353760 0.1866295 0.36575340 0.563888900

dim(y)

## [1] 48803    12

expressed <- rowSums(y$other$Detection < 0.05) >= 3 ;table(expressed)

## expressed
## FALSE  TRUE 
## 24112 24691

y <- y[expressed,]
dim(y)

## [1] 24691    12

这时的表达矩阵，也就能对接常规芯片数据的其他分析了，比如差异分析。

exp = y$E
exp[1:4,1:4]

##                      1        2        3        4
## ILMN_2055271  5.085517 5.294789 5.047373 5.274919
## ILMN_1653355  5.803398 5.901535 5.378863 5.567170
## ILMN_1787689  5.164395 4.929249 5.401202 5.006137
## ILMN_1745607 10.508010 9.922797 5.689654 5.317240

2.GEO数据实例

去GSE16997的页面下载它的补充文件“GSE16997_raw.txt”。

读取、背景校正和标准化

rm(list = ls())
x <- read.ilmn(files="GSE16997_raw.txt",
               expr="Sample",
               probeid="ID_REF",
               other.columns="Detection Pval")

## Reading file GSE16997_raw.txt ... ...

y <- neqc(x,detection.p="Detection Pval")

## Note: inferring mean and variance of negative control probe intensities from the detection p-values.

探针过滤

x$other$Detection[1:4,1:4]

##                       1         2          3           4
## ILMN_1762337 0.55849580 0.6754875 0.13698630 0.601388900
## ILMN_2055271 0.03064067 0.0000000 0.04931507 0.002777778
## ILMN_1736007 0.27715880 0.2924791 0.15342470 0.486111100
## ILMN_2383229 0.47353760 0.1866295 0.36575340 0.563888900

dim(y)

## [1] 48803    12

expressed <- rowSums(y$other$`Detection Pval` < 0.05) >= 3 ;table(expressed)

## expressed
## FALSE  TRUE 
## 24112 24691

y <- y[expressed,]
dim(y)

## [1] 24691    12

y$E[1:4,1:4]

##                      1        2        3        4
## ILMN_2055271  5.085517 5.294789 5.047373 5.274919
## ILMN_1653355  5.803398 5.901535 5.378863 5.567170
## ILMN_1787689  5.164395 4.929249 5.401202 5.006137
## ILMN_1745607 10.508010 9.922797 5.689654 5.317240