需求
和這一篇的需求一致，不同公司的芯片原始數(shù)據(jù)處理方式也不同，最終都是為了得到表達矩陣，今天介紹的是illumina beadchip 芯片，用神奇的limma搞定它。
找不到某些GEO數(shù)據(jù)的表達矩陣腫么辦
illumina beadchip芯片原始數(shù)據(jù)處理,參考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.
聽起來很厲害，實際上一個函數(shù)搞定。
探針過濾
原始文件里的信息，被拆分成了一個表達矩陣和一個Detection P值矩陣，P值可以用于過濾樣本。
示例數(shù)據(jù)總共有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
這時的表達矩陣，也就能對接常規(guī)芯片數(shù)據(jù)的其他分析了，比如差異分析。
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數(shù)據(jù)實例
去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
