需求
和这一篇的需求一致，不同公司的芯片原始数据处理方式也不同，最终都是为了得到表达矩阵，今天介绍的是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
